text:
# function transpose_kernel!(
#         state, At, A::AbstractArray{T}, width, height, ::Val{BLOCK}, ::Val{LMem}
#     ) where {BLOCK, LMem, T}
#
#     bidx_x = blockidx_x(state) - 1
#     bidx_y = blockidx_y(state) - 1
#     tidx_x = threadidx_x(state) - 1
#     tidx_y = threadidx_y(state) - 1
#
#     A_local = @LocalMemory(state, T, LMem)
#
#     base_idx_a = bidx_x * BLOCK + bidx_y * (BLOCK * width)
#     base_idx_a_t = bidx_y * BLOCK + bidx_x * (BLOCK * height)
#
#     glob_idx_a = base_idx_a + tidx_x + width * tidx_y
#     glob_idx_a_t = base_idx_a_t + tidx_x + height * tidx_y
#     glob_idx_a >= length(A) && return
#     A_local[tidx_y * BLOCK + tidx_x + 1] = A[glob_idx_a + 1]
#     synchronize_threads(state)
#     At[glob_idx_a_t + 1] = A_local[tidx_x * BLOCK + tidx_y + 1]
#     return
# end

function transpose_blocks!(
        state, odata::AbstractArray{T}, idata, ::Val{SHMEM}, ::Val{TDIM}, ::Val{BLOCK_ROWS}, ::Val{NROW}
    ) where {T, SHMEM, TDIM, BLOCK_ROWS, NROW}

    tile = @LocalMemory(state, T, SHMEM)
    bidx_x = blockidx_x(state) - 1
    bidx_y = blockidx_y(state) - 1
    tidx_x = threadidx_x(state) - 1
    tidx_y = threadidx_y(state) - 1

    x = bidx_x * TDIM + tidx_x + 1
    y = bidx_y * TDIM + tidx_y + 1
    dims = size(idata)

    (x <= dims[2] && (y + (BLOCK_ROWS * 3)) <= dims[1]) || return

    for j = 0:3
        j0 = j * BLOCK_ROWS
        @inbounds tile[tidx_x + 1, tidx_y + j0 + 1] = idata[y + j0, x]
    end

    synchronize_threads(state)
    for j = 0:3
        j0 = j * BLOCK_ROWS
        @inbounds odata[x, y + j0] = tile[tidx_x + 1, tidx_y + j0 + 1]
    end

    return
end

function LinearAlgebra.transpose!(At::GPUArray{T, 2}, A::GPUArray{T, 2}) where T
    if size(A, 1) == size(A, 2) && all(x-> x % 32 == 0, size(A))
        outsize = size(At)
        TDIM = 32; BLOCK_ROWS = 8
        nrows = TDIM ÷ BLOCK_ROWS
        shmemdim = (TDIM, (TDIM + 1))
        static_params = map(x-> Val(x), (shmemdim, TDIM, BLOCK_ROWS, nrows))
        args = (At, A, static_params...)

        griddim = ceil.(Int, size(A) ./ (TDIM, TDIM))
        blockdim = (TDIM, BLOCK_ROWS)
        # optimized version for 32x & square dimensions
        gpu_call(transpose_blocks!, At, args, (griddim, blockdim))
    else
        # simple fallback
        gpu_call(At, (At, A)) do state, At, A
            idx = @cartesianidx A state
            @inbounds At[idx[2], idx[1]] = A[idx[1], idx[2]]
            return
        end
    end
    At
end

function genperm(I::NTuple{N}, perm::NTuple{N}) where N
    ntuple(d-> (@inbounds return I[perm[d]]), Val(N))
end

function LinearAlgebra.permutedims!(dest::GPUArray, src::GPUArray, perm) where N
    perm isa Tuple || (perm = Tuple(perm))
    gpu_call(dest, (dest, src, perm)) do state, dest, src, perm
        I = @cartesianidx src state
        @inbounds dest[genperm(I, perm)...] = src[I...]
        return
    end
    return dest
end

function Base.copyto!(A::AbstractArray, B::Adjoint{<:Any, <:GPUArray})
    copyto!(A, Adjoint(Array(parent(B))))
end
function Base.copyto!(A::AbstractArray, B::Transpose{<:Any, <:GPUArray})
    copyto!(A, Transpose(Array(parent(B))))
end
function Base.copyto!(A::AbstractArray, B::UpperTriangular{<:Any, <:GPUArray})
    copyto!(A, UpperTriangular(Array(parent(B))))
end
function Base.copyto!(A::AbstractArray, B::LowerTriangular{<:Any, <:GPUArray})
    copyto!(A, LowerTriangular(Array(parent(B))))
end

function Base.copyto!(A::GPUArray, B::Adjoint{T, <: GPUArray}) where T
    transpose!(A, B.parent)
end


text:
function random_spatial_network(n::Integer; # number of nodes
    d::Integer = 2, # dimension
    Pk = Poisson(10), # Expected degree distribution
    Pf = Normal(0.0, 0.03) # Distance kernel
    )
    ks = rand(Pk, n)
    xs = rand(n,d)
    dists = Array{Float64,2}(undef, n, n)
    for i in 1:n, j in i:n
        dists[i,j] = norm(xs[i,:]-xs[j,:])
        dists[j,i] = dists[i,j]
    end
    f = x -> pdf(Pf, x)
    return random_spatial_network(dists, ks, f, 1), xs
end

function random_spatial_network(distances::AbstractMatrix, ks, f, ρ)
    k = mean(ks)
    ps = (ks * ks') .* f.(distances) / (ρ * k)
    return edge_probability_graph(ps)
end

function edge_probability_graph(ps::AbstractMatrix)
    n = size(ps, 1)
    edges = [(i,j) for i in 1:n-1, j in (i+1):n if rand() < ps[i,j]]
    return SimpleGraphFromIterator(edges)
end

function geometric_inhomogeneous_random_graph(n, τ, α)

end

ext: jl


text:
using Knet
using ReinforcementLearning
using ArgParse
using JLD
import Base.push!
import Base.length

const SCREEN_WIDTH = 600
const BATCH_SIZE = 128
const GAMMA = 0.95
const EPS_INIT = 0.9
const EPS_FINAL = 0.005
const EPS_DECAY = 200
const NUM_EPISODES = 10
const LR = 0.001
const GCLIP = 5.0
const WINIT = 0.1
const CAPACITY = 10000
const PERIOD = 200

function main(args)
    s = ArgParseSettings()
    s.description = ""

    @add_arg_table s begin
        ("--hidden"; nargs='*'; arg_type=Int64)
        ("--lr"; default=LR; arg_type=Float64)
        ("--gclip"; default=GCLIP; arg_type=Float64)
        ("--batchsize"; default=BATCH_SIZE; arg_type=Int64)
        ("--gamma"; default=GAMMA; arg_type=Float64)
        ("--epsinit"; default=EPS_INIT; arg_type=Float64)
        ("--epsfinal"; default=EPS_FINAL; arg_type=Float64)
        ("--epsdecay"; default=EPS_DECAY; arg_type=Int64)
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}":"Array{Float32}"))
        ("--winit"; default=WINIT; arg_type=Float64)
        ("--capacity"; default=CAPACITY; arg_type=Int64)
        ("--numepisodes"; default=NUM_EPISODES; arg_type=Int64)
        ("--monitor"; default=nothing)
        ("--seed"; default=-1; arg_type=Int64)
        ("--generate"; default=0; arg_type=Int64)
        ("--loadfile"; default=nothing)
        ("--savefile"; default=nothing)
        ("--period"; default=PERIOD; arg_type=Int64)
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:atype] = eval(parse(o[:atype]))
    sr = o[:seed] > 0 ? srand(o[:seed]) : srand()

    env = GymEnv("CartPole-v0")
    w = opts = nothing
    if o[:loadfile] == nothing
        w = initweights(o[:atype], o[:hidden], o[:winit])
    else
        w = load(o[:loadfile], "w")
        w = map(wi->convert(o[:atype], wi), w)
    end
    opts = map(wi->Adam(;gclip=o[:gclip]), w)

    memory = ReplayMemory(o[:capacity])
    if o[:monitor] != nothing
        monitor_start(env, o[:monitor])
    end

    steps_done = 0
    currerr = 0
    preverr = 0
    iter = 0
    for k = 1:o[:numepisodes]
        # reset environment
        s0 = getInitialState(env); state = s0

        while true
            s = convert(o[:atype], state.data)
            action = select_action(w,s,steps_done; o=o)
            steps_done += 1

            next_state, reward = transfer(env, state, env.actions[action])
            t = Transition(state, action, next_state, reward)
            push!(memory, t)
            state = next_state

            if length(memory) >= o[:batchsize]
                val = train!(w,memory,opts; o=o)

                preverr = currerr
                if iter < 100
                    currerr = (iter * currerr) + val
                    currerr = currerr / (iter+1)
                else
                    currerr = 0.01*val+0.99*currerr
                end

                iter += 1
                if iter % o[:period] == 0
                    println("($iter,$currerr)")
                end
            end

            if state.done
                break
            end
        end

        if o[:savefile] != nothing
            save(o[:savefile], "w", map(wi->Array(wi),w))
        end
    end

    for k = 1:o[:generate]
        done = false
        state = getInitialState(env)
        while !done
            s = convert(o[:atype], state.data)
            action = select_action(w,s,steps_done; o=o)
            steps_done += 1
            render_env(env)
            next_state, reward = transfer(env, state, env.actions[action])
            state = next_state
            done = next_state.done
        end
    end

    if o[:monitor] != nothing
        monitor_close(env)
    end
end

function initweights(atype, hidden, winit=0.1)
    w = []
    x = 4
    for y in [hidden..., 2]
        push!(w, winit*randn(y,x))
        push!(w, zeros(y,1))
        x = y
    end
    return map(wi->convert(atype,wi), w)
end

function predict(w,x)
    for i=1:2:length(w)
        x = w[i]*x .+ w[i+1]
        if i < length(w)-1
            x = relu(x)
        end
    end
    return x
end

type ReplayMemory
    capacity
    memory

    function ReplayMemory(capacity)
        memory = Transition[]
        new(capacity, memory)
    end
end

type Transition
    state
    action
    next_state
    reward
end

function push!(obj::ReplayMemory, t::Transition)
    push!(obj.memory,t)
    length(obj.memory) > obj.capacity && shift!(obj.memory)
end

function length(obj::ReplayMemory)
    return length(obj.memory)
end

function sample(obj::ReplayMemory, batchsize)
    batchsize = min(length(obj), batchsize)
    indices = randperm(length(obj))[1:batchsize]
    return obj.memory[indices]
end

function select_action(w, s, steps_done; o=Dict())
    ei = get(o, :epsinit,  EPS_INIT)
    ef = get(o, :epsfinal, EPS_FINAL)
    ed = get(o, :epsdecay, EPS_DECAY)
    et = ef + (ei-ef) * exp(-steps_done/ed)

    if rand() > et
        # @show w,s
        s = reshape(s,length(s),1)
        y = predict(w,s)
        y = Array(y)
        return indmax(y)
    else
        return rand(1:2)
    end
end

function train!(w, memory, opts; o=Dict())
    # get params
    gamma = get(o, :gamma, GAMMA)
    atype = get(o, :atype, Array{Float32})
    batchsize = get(o, :batchsize, BATCH_SIZE)

    # sample batch
    batch = sample(memory, batchsize)

    # prepare batch
    states  = mapreduce(i->i.state.data, hcat, batch)
    actions = map(i->i.action, batch)
    rewards = map(i->i.reward, batch)
    mask = map(i->!i.next_state.done, batch)
    mask = find(mask)
    next_states = filter(i->!i.next_state.done, batch)
    next_states = mapreduce(i->i.next_state.data, hcat, next_states)

    # convert batch
    states = convert(atype, states)
    rewards = convert(atype, rewards)
    next_states = convert(atype, next_states)

    # expected return
    y0 = predict(w,next_states)
    y1 = maximum(y0,1)
    y2 = sum(y0 .* (y0.==y1), 1)
    y3 = reshape(y2, 1, length(y2))
    expret = reshape(rewards, 1, length(rewards))
    expret[:,mask] += gamma * y3

    # training
    values = []
    g = gradient(w, states, actions, expret, mask; values=values)
    update!(w,g,opts)
    return values[1]
end

function objective(w, states, actions, expret, mask; values=[])
    qsa = predict(w, states)
    nrows,ncols = size(qsa)
    index = actions + nrows*(0:(length(actions)-1))
    qs  = qsa[index]
    qs  = reshape(qs, 1, length(qs))
    val = sumabs2(expret-qs) / length(qs)
    push!(values, val)
    return val
end

gradient = grad(objective)

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/rl

max_issues_repo_name: ilkerkesen/rl

max_forks_repo_name: ilkerkesen/rl


text:
module Utils
# Module containing the necessary structs and functions for LIIC execution

using LinearAlgebra
using Statistics
using Formatting

export import_molecule, write_file

struct ZMatrix
    Atoms::Vector{String}
    IntVars::Vector{Float64}
    VarNames::Vector{String}
    Number::Int64
end

struct Molecule
    Atoms::Vector{String}
    Coord::Array{Float64,2}
    Number::Int64
end

function translate_to_centroid(coord_matrix)
    # Normalises the molecular coordinates by centering them.
    center = mean(coord_matrix,dims=1)
    translated_geom = coord_matrix .- center
    return translated_geom
end

function optimal_rotation_matrix(CCmatrix)
    # Returns 3x3 matrix that can be applied to P to get Q
    # Previous implementation, not general (for cases of non-invertible matrices)
    # ORmatrix = sqrt(transpose(CCmatrix)*CCmatrix)*inv(CCmatrix)
    # Using SVD instead
    u,sing,vt = svd(CCmatrix)
    if sign(det(transpose(vt)*transpose(u))) == 1.0
        ORmatrix = transpose(vt)*transpose(u)
    elseif sign(det(transpose(vt)*transpose(u))) == -1.0
        mat = [1. 0. 0.;0. 1. 0.; 0. 0. -1.]
        ORmatrix = transpose(vt)*mat*transpose(u)
    else
        println("Error: Issue with SVD routine in finding optimal rotation matrix")
    end
    return ORmatrix
end

function kabsch_rotate(Pgeom,Qgeom)
    # Returns the two geometries in xyz format, one of which has been translated, the
    # other translated and rotated
    # Pgeom = import_molecule(Pxyz)
    # Qgeom = import_molecule(Qxyz)

    normalisedP = (translate_to_centroid(Pgeom.Coord))
    normalisedQ = (translate_to_centroid(Qgeom.Coord))

    # Calculate cross covariance matrix
    xcov = transpose(normalisedP)*normalisedQ
    orot = optimal_rotation_matrix(xcov)

    num_atoms = size(normalisedP)[1]
    rotated = zeros(Float64,num_atoms,3)

    for i=1:num_atoms
        rotated[i,:] = orot*normalisedP[i,:]
    end

    RMSD_value = norm(rotated-normalisedQ)

    println("The RMSD between these two structures is ",RMSD_value)

    RotatedP = Molecule(Pgeom.Atoms,rotated,num_atoms)
    NormalisedQ = Molecule(Qgeom.Atoms,normalisedQ,num_atoms)

    return RotatedP, NormalisedQ
end

function import_molecule(xyz_file::String)
    io = open(xyz_file)
    natoms = parse(Int,readline(io))
    readline(io);
    reg=r"([a-zA-Z]+)\s+(-?\d+.?\d*)\s+(-?\d+.?\d*)\s+(-?\d+.?\d*)"
    atoms=Vector{String}()
    atom_coords=Array{Float64,2}(undef,natoms,3)
    for num in 1:natoms
        m=match(reg,readline(io))
        push!(atoms,string(m.captures[1]::SubString))
        atom_coords[num,:]=[parse(Float64,m.captures[2]),parse(Float64,m.captures[3]),parse(Float64,m.captures[4])]
    end
    close(io)
    return Molecule(atoms,atom_coords,natoms)
end

function bond(a1::Vector{Float64},a2::Vector{Float64})
    return norm(a1-a2)::Float64
end

function angle(a1::Vector{Float64},a2::Vector{Float64},a3::Vector{Float64})
    bond1 = a3 - a2
    bond2 = a2 - a1
    sintheta = norm(cross(bond1,bond2))
    costheta = dot(bond1,bond2)
    theta = atan(sintheta,costheta)
    angle = 180.0 * theta/pi
    # angle = acosd(dot(bond1,bond2)/(norm(bond1)*norm(bond2)))
    return angle::Float64
end

function dihedral(a1::Vector{Float64},a2::Vector{Float64},a3::Vector{Float64},a4::Vector{Float64})
    bond1 = a4-a3
    bond2 = a3-a2
    bond3 = a2-a1
    plane1 = cross(bond1,bond2)
    plane2 = cross(bond2,bond3)
    dihedral = acosd(dot(plane1,plane2)/(norm(plane1)*norm(plane2)))
    if dot(cross(plane1,plane2),bond2) < 0
        dihedral = -dihedral
    end
    return dihedral::Float64
end

function write_zmat(zmat::ZMatrix,withvars::Bool=true)
    if withvars==true
        natoms = zmat.Number
        for atom in 1:natoms
            if atom < 2
                printfmtln("{:3s}",zmat.Atoms[atom])
            elseif atom < 3
                printfmtln("{:3s} {:>4d}  {:11s}",zmat.Atoms[atom],(atom-1),zmat.VarNames[atom-1])
            elseif atom < 4
                printfmtln("{:3s} {:>4d}  {:11s} {:>4d}  {:11s}",zmat.Atoms[atom],(atom-1),
                zmat.VarNames[atom-1],(atom-2),zmat.VarNames[atom])
            else
                printfmtln("{:3s} {:>4d}  {:11s} {:>4d}  {:11s} {:>4d}  {:11s}",zmat.Atoms[atom],
                (atom-1),zmat.VarNames[3*atom-8],(atom-2),zmat.VarNames[3*atom-7],(atom-3),zmat.VarNames[3*atom-6])
            end
        end
        println("  Variables:")
        for var in 1:length(zmat.IntVars)
            printfmtln("{:3s} {:>11.8f}",zmat.VarNames[var],zmat.IntVars[var])
        end
    else
        natoms = length(zmat.Atoms)
        for atom in 1:natoms
            if atom < 2
                printfmtln("{:3s}",zmat.Atoms[atom])
            elseif atom < 3
                printfmtln("{:3s} {:>4d}  {:>11.8f}",zmat.Atoms[atom],(atom-1),zmat.IntVars[atom-1])
            elseif atom < 4
                printfmtln("{:3s} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f}",zmat.Atoms[atom],(atom-1),
                zmat.IntVars[atom-1],(atom-2),zmat.IntVars[atom])
            else
                printfmtln("{:3s} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f}",zmat.Atoms[atom],
                (atom-1),zmat.IntVars[3*atom-8],(atom-2),zmat.IntVars[3*atom-7],(atom-3),zmat.IntVars[3*atom-6])
            end
        end
    end
end

function write_xyz(mol::Molecule,io=nothing)
    natoms = mol.Number
    if io !== nothing
        printfmtln(io,"{:>4d}",natoms)
        println(io)
        for atom in 1:natoms
            printfmtln(io,"{:3s}   {:>11.8f}   {:>11.8f}    {:>11.8f}",mol.Atoms[atom],mol.Coord[atom,1],mol.Coord[atom,2],mol.Coord[atom,3])
        end
    elseif io === nothing
        printfmtln("{:>4d}",natoms)
        println()
        for atom in 1:natoms
            printfmtln("{:3s}   {:>11.8f}   {:>11.8f}    {:>11.8f}",mol.Atoms[atom],mol.Coord[atom,1],mol.Coord[atom,2],mol.Coord[atom,3])
        end
    else
        println("Error in write_xyz")
    end
end

function write_liic(outArray,filename::String)
    io = open(filename,"w")
    for i in 1:length(outArray)
        write_xyz(outArray[i],io)
    end
    close(io)
end

function xyz_to_zmat(molecule)
    natoms     = molecule.Number
    atoms      = molecule.Atoms
    intvars      = Vector{Float64}()
    varnames   = Vector{String}()
    for atom in 1:natoms
        if atom > 1
            bl     = bond(molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,bl)
            push!(varnames,"R$atom")
        end
        if atom > 2
            agl  = angle(molecule.Coord[atom-2,:],molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,agl)
            push!(varnames,"A$atom")
        end
        if atom > 3
            dhl  = dihedral(molecule.Coord[atom-3,:],molecule.Coord[atom-2,:],molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,dhl)
            push!(varnames,"D$atom")
        end
    end
    zmatrix    = ZMatrix(atoms,intvars,varnames,natoms)
    return zmatrix
end

function rotation_matrix(axis,angle)
    # Rotation about the origin in 3 dimensions using the Euler-Rodrigues formula
    axis = -axis / norm(axis)
    a = cosd(angle*0.5)
    (b,c,d) = axis * sind(angle*0.5)
    rotM = [a*a+b*b-c*c-d*d 2*(b*c-a*d) 2*(b*d+a*c);
    2*(b*c+a*d) a*a+c*c-b*b-d*d 2*(c*d-a*b);
    2*(b*d-a*c) 2*(c*d+a*b) a*a+d*d-b*b-c*c]
    return rotM::Matrix{Float64}
end

function zmat_to_xyz(zmat::ZMatrix)
    natoms=zmat.Number
    atom_coords = Array{Float64,2}(undef,natoms,3)

    if natoms > 0
        atom_coords[1,:]=[0.0,0.0,0.0]

        if natoms > 1
            r2 = zmat.IntVars[1]
            atom_coords[2,:]=[r2,0.0,0.0]

            if natoms > 2
                vec = atom_coords[2,:] - atom_coords[1,:]
                r3 = zmat.IntVars[2]
                a3 = zmat.IntVars[3]
                vec = r3 * vec / norm(vec)
                atom_coords[3,:] = atom_coords[2,:] + (rotation_matrix([0, 0, 1], a3) * vec)

                if natoms > 3
                    for atom in 4:natoms
                        rn = zmat.IntVars[3*atom-8]
                        an = zmat.IntVars[3*atom-7]
                        dn = zmat.IntVars[3*atom-6]

                        # ox = rn * cosd(an)
                        # oy = rn * cosd(dn) * sind(an)
                        # oz = rn * sind(an) * sind(dn)

                        # Coordinates of the last three atoms
                        am1 = atom_coords[atom-1,:]
                        am2 = atom_coords[atom-2,:]
                        am3 = atom_coords[atom-3,:]

                        # Vectors between them
                        vec1 = am1 - am2
                        vec2 = am2 - am3
                        # Plane normal defined by these vectors
                        nm = cross(vec1,vec2)

                        # Vector of length rn pointing along vec1
                        vecrn = rn * vec1 / norm(vec1)

                        # Rotate vecrn by angle an around the normal
                        vecan = rotation_matrix(nm,an)*vecrn

                        # Rotate vecan around dihedral angle
                        vecdn = rotation_matrix(vec1,dn)*vecan

                        # Add this vector to previous coordinate
                        atom_coords[atom,:] = am1 + vecdn
                    end
                end
            end
        end
    end

    mol=Molecule(zmat.Atoms,atom_coords,natoms)

    return mol::Molecule
end

function cartesian(mol1::Molecule,mol2::Molecule,steps::Int64)
	"""
    This function returns an "steps" length array of Molecule objects
	"""
    # difference = (mol2.Coord - mol1.Coord) / (steps - 1)
    natoms = mol1.Number
	liic = Array{Float64,3}(undef,natoms,3,steps)
    arrayOfMolecules = Array{Molecule, 1}(undef, steps)
    for step in 1:steps
		fracA = (steps-step)/(steps-1)
		fracB = (step-1)/(steps-1)
		liic[:,:,step] = mol1.Coord*fracA + mol2.Coord*fracB
        arrayOfMolecules[step] = Molecule(mol1.Atoms,liic[:,:,step],natoms)
    end
    return arrayOfMolecules
end

function internal_babel(int_arr_1,int_arr_2,steps,header)
	"""
	Takes array of Z-matrix internal coordinates and interpolates between them
	to return an array of xyzs.
	"""
	numInternalCrd = size(int_arr_1)[1]
	natoms = Int((numInternalCrd + 6)/3)
	diffVec = Vector(undef,numInternalCrd)
	for i in 1:numInternalCrd
		if ( -180. <= (int_arr_2[i,2]-int_arr_1[i,2]) <= 180. )
			diffVec[i] = (int_arr_2[i,2]-int_arr_1[i,2])/steps
		elseif (int_arr_2[i,2] - int_arr_1[i,2]) < -180.
			diffVec[i] = (int_arr_2[i,2] - int_arr_1[i,2] + 360.)/steps
		elseif (int_arr_2[i,2]-int_arr_1[i,2]) > 180.
			diffVec[i] = (int_arr_2[i,2] - int_arr_1[i,2] - 360.)/steps
		else
			println("Something has gone terribly wrong")
		end
	end
	zmatArray = Vector(undef,steps)
	interXYZ = Array{Float32,3}(undef,natoms,3,steps)
	atom_names = Array{String,1}(undef,natoms)
	for j in 1:steps
		zmatArray[j] = header * "Variables:\n" * intlVec(int_arr_1,diffVec,j)
		open("tmp","w") do io
			write(io,zmatArray[j])
		end
		convXyz = read(`obabel -igzmat tmp -oxyz`,String)
		run(`rm tmp`)
		just_coords = Array{Float32}(readdlm(IOBuffer(convXyz),skipstart=2)[:,2:4])
		atom_names = Array{String}(readdlm(IOBuffer(convXyz),skipstart=2)[:,1])
		interXYZ[:,:,j] = just_coords
	end
	return interXYZ, atom_names
end

function internal(zmat_f::ZMatrix,zmat_l::ZMatrix,steps::Int64)
	int_f = zmat_f.IntVars
	int_l = zmat_l.IntVars
	natoms = zmat_l.Number
	nvars = length(zmat_l.IntVars)
	difference = Vector(undef,nvars)
	for i in 1:nvars
		toAngle = (zmat_l.IntVars[i]+180.)
		fromAngle = (zmat_f.IntVars[i]+180.)
		diff = abs(toAngle - fromAngle)
		if diff < 180.
			difference[i] = toAngle - fromAngle
		else
			if fromAngle > toAngle
				fromAngle = fromAngle - 360.
				difference[i] = toAngle - fromAngle
			elseif toAngle > fromAngle
				toAngle = toAngle - 360.
				difference[i] = toAngle - fromAngle
			end
		end
	end
	difference = difference / (steps-1)
	liic = Array{Float64,2}(undef,nvars,steps)
	arrayOfMolecules = Array{Molecule, 1}(undef, steps)
	for step in 1:steps
		liic[:,step] = int_f + (step-1)*difference
		arrayOfMolecules[step] = zmat_to_xyz(ZMatrix(zmat_f.Atoms,liic[:,step],zmat_f.VarNames,natoms))
	end
	# @assert isapprox(int_l,liic[:,steps],atol=1e-5)
	return arrayOfMolecules
end

end

max_stars_repo_name: dasha-shchep/LIIC

max_issues_repo_name: dasha-shchep/LIIC

max_forks_repo_name: dasha-shchep/LIIC


text:
# Basic acquisition example
using AcqSynth, PlotlyJS, Statistics
const acq = AcqSynth

# Let's see if we have a board plugged in
acq.get_num_boards()
# If the number you get is bigger than 0, congrats your board works!

# Select board 0 and check its serial number
boardnum = 0
acq.get_serial(boardnum)

# OK we can start working. First, we need to calibrate/setup the board.
acq.init_board(boardnum)

# Make sure we have an ADC12D2000!
acq.is_AD12(boardnum)

# We can see have many channels the board has
acq.get_all_channels(boardnum)

# Now, we can configure the default acquisition session and digitize data
# Let's start by checking what number of channels are in use right now
acq.get_num_channels(boardnum)

# Check the frequency
acq.get_frequency(boardnum)

# At this point, without any other configuration we can start acquiring by
# choosing how many 1MiB blocks we want and sending a soft trigger.
# This is done with the setup_acquire() function.
numblocks = 4 # let's grab a bunch of blocks
acq.setup_acquire(boardnum,numblocks) # acquisition starts here!

# Grab the data from the card with get_samples_12 utility function
samples = acq.get_samples_12(boardnum,numblocks) # transfer data as 12-bit samples
data = permutedims(reshape(samples,2,:))

# Plot the stuff
t = (1:size(data,1))./2E9 # time vector, we are sampling at 2GSPS
plot([scatter(x=t[1:2000],y=data[1:2000,1],name="AIN0"),
      scatter(x=t[1:2000],y=data[1:2000,2],name="AIN1")])

println(" AIN0: $(mean(data[:,1])) ± $(std(data[:,1]))")
println(" AIN1: $(mean(data[:,2])) ± $(std(data[:,2]))")


text:
#///////////////////////////////////////
#// File Name: type_definition.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: Type definitions for Risk Sensitive Stochastic SAC
#///////////////////////////////////////

import Base: ==, isequal, isapprox
import RobotOS: Time

using DataStructures
using LinearAlgebra
using RobotOS

abstract type Parameter end;

abstract type State end;

# WayPoint
struct WayPoint2D
    pos::Vector{Float64}
    t::Time
    function WayPoint2D(pos::Vector{Float64}, t::Time)
        @assert length(pos) == 2 "Invalid waypoint dimension!"
        return new(pos, t)
    end
end

# Trajectory
const Trajectory2D = OrderedDict{Time, Vector{Float64}};

function Trajectory2D(points::Vector{WayPoint2D})
    trajectory = Trajectory2D();
    for waypoint in points
        trajectory[waypoint.t] = waypoint.pos;
    end
    sort!(trajectory)
    return trajectory
end

function get_position(trajectory::Trajectory2D, t::Time)
    if haskey(trajectory, t)
        return trajectory[t]
    else
        traj_array = collect(trajectory);
        index2 = findfirst(x -> x[1] > t, traj_array);
        @assert !(index2 == 1 || index2 == nothing) "Trajectory Interpolation failed for t = $(to_sec(t)) [s]!"
        t2 = traj_array[index2][1];
        index1 = index2 - 1;
        t1 = traj_array[index1][1];
        # linear interpolation
        pos = traj_array[index1][2] + to_nsec(t - t1)/to_nsec(t2 - t1)*
              (traj_array[index2][2] - traj_array[index1][2])
        return pos
    end
end

# Ego Robot State
struct RobotState <: State
    x::Vector{Float64} # [pos_x, pos_y, vel_x, vel_y]
    t::Time
    function RobotState(x::Vector{Float64}, t::Time)
        @assert length(x) == 4 "Invalid ego state dimension!"
        return new(x, t)
    end
end

RobotState(x::Vector{Float64}) = RobotState(x, Time());

function isequal(e1::RobotState, e2::RobotState)
    return e1.x == e2.x && e1.t == e2.t
end

==(e1::RobotState, e2::RobotState) = isequal(e1, e2)

function isapprox(e1::RobotState, e2::RobotState; kwargs...)
    return isapprox(e1.x, e2.x; kwargs...) && e1.t == e2.t
end

get_position(e::RobotState) = e.x[1:2];
get_velocity(e::RobotState) = e.x[3:4];

# Joint State
mutable struct WorldState <: State
    e_state::RobotState
    ap_dict::Dict{String, Vector{Float64}}
    t::Time
    t_last_m::Union{Time, Nothing}
    function WorldState(e_state::RobotState,
                        ap_dict::Dict{String, Vector{Float64}},
                        t::Time, t_last_m::Union{Time, Nothing})
        @assert all(length.(values(ap_dict)) .== 2) "Invalid ado state dimension!"
        @assert e_state.t == t "Inconsistent time stamps!";
        if t_last_m != nothing
            @assert t_last_m <= t "Last measurement time cannot be in future!"
        end
        return new(e_state, ap_dict, t, t_last_m)
    end
end

function WorldState(e_state::RobotState,
                    ap_dict::Dict{String, Vector{Float64}},
                    t_last_m=nothing)
    return WorldState(e_state, ap_dict, e_state.t, t_last_m)
end

function WorldState(x::Vector{Float64},
                    ap_dict::Dict{String, Vector{Float64}},
                    t::Time=Time(), t_last_m=nothing)
    return WorldState(RobotState(x, t), ap_dict, t, t_last_m)
end
function isequal(w1::WorldState, w2::WorldState)
    return begin w1.e_state == w2.e_state && w1.ap_dict == w2.ap_dict &&
                 w1.t == w2.t && w1.t_last_m == w2.t_last_m
           end
end

==(w1::WorldState, w2::WorldState) = isequal(w1, w2);

function isapprox(w1::WorldState, w2::WorldState; kwargs...)
    e_approx = isapprox(w1.e_state, w2.e_state; kwargs...);
    ap_approx = keys(w1.ap_dict) == keys(w2.ap_dict) &&
                all([isapprox(w1.ap_dict[key], w2.ap_dict[key]; kwargs...)
                     for key in keys(w1.ap_dict)]);
    return e_approx && ap_approx && w1.t == w2.t && w1.t_last_m == w2.t_last_m
end


text:
type AnnealingEpsilonGreedy <: BanditAlgorithm
  counts::Vector{Int64}
  values::Vector{Float64}
end

function AnnealingEpsilonGreedy(n_arms::Int64)
  AnnealingEpsilonGreedy(zeros(Int64, n_arms), zeros(n_arms))
end

function initialize(algo::AnnealingEpsilonGreedy, n_arms::Int64)
  algo.counts = zeros(Int64, n_arms)
  algo.values = zeros(n_arms)
end

function initialize(algo::AnnealingEpsilonGreedy, arms::Vector)
  initialize(algo, length(arms))
end

function select_arm(algo::AnnealingEpsilonGreedy)
  t = sum(algo.counts) + 1
  epsilon = 1 / log(t + eps(1.0))

  if rand() > epsilon
    return ind_max(algo.values)
  else
    return randi(length(algo.values))
  end
end

function update(algo::AnnealingEpsilonGreedy, chosen_arm::Int64, reward::Real)
  algo.counts[chosen_arm] += 1
  n = algo.counts[chosen_arm]

  value = algo.values[chosen_arm]
  algo.values[chosen_arm] = ((n - 1) / n) * value + (1 / n) * reward
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
using JsonGrinder, JSON, Test, SparseArrays

@testset "Generating HTML Schema" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a":1}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")
	j6 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a":1}], "b":1}""")
	j7 = JSON.parse("""{"a": [{"a":4,"b":5},{"b":6,"a":7}], "b":2}""")
	j8 = JSON.parse("""{"a": [{"a":9,"b":10},{"b":11,"a":12}], "b":2}""")
	j9 = JSON.parse("""{"a": [{"a":4,"b":3},{"b":2,"a":2}], "b":2}""")
	j10 = JSON.parse("""{"a": [{"a":11,"b":12},{"b":13,"a":14}], "b":2}""")
	j11 = JSON.parse("""{"a": [{"a":7,"b":5},{"b":6,"a":6}], "b":2}""")

	JsonGrinder.updatemaxkeys!(6)

	sch = JsonGrinder.schema([j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11])

	JsonGrinder.generate_html(sch, "schema_test_maxvals=5.html", max_vals=5)
	@test isfile("schema_test_maxvals=5.html")
end


text:
using ChunkedArrays
using Base.Test

tic()
u = [2 2
     2 2]

chunkRand = ChunkedArray(randn,float(u),10)

println(next(chunkRand))

chunkRand2 = ChunkedArray(randn,(3,3),100)

println(next(chunkRand2))

chunkRand3 = ChunkedArray(randn,100)

println(next(chunkRand3))

u = [2;2;2;2]
chunkRand4 = ChunkedArray(randn,float(u),10)
println(next(chunkRand4))

println("Random Generating Benchmark")
const loopSize = 1000
const buffSize = 100
const numRuns = 40

function test1()
  j=[0;0;0;0]
  for i = 1:loopSize
    j += randn(4)
  end
end

function test2()
  j=[0;0;0;0]
  chunks = 1000
  rands = randn(4,chunks)
  for k = 1:(loopSize÷chunks)
    rands[:] = randn(4,chunks)
    for i = 1:chunks
      j += rands[:,i]
    end
  end
end

function test3()
  j=[0;0;0;0]
  rands = randn(4,loopSize)
  for i = 1:loopSize
    j += rands[:,i]
  end
end

function test4()
  j=[0;0;0;0]
  chunks = 100
  rands = randn(4,chunks)
  for k = 1:(loopSize÷chunks)
    rands[:] = randn(4,chunks)
    for i = 1:chunks
      j += rands[:,i]
    end
  end
end

function test5()
  rands = randn(4,loopSize)
  j=[0;0;0;0]
  for i = 1:loopSize
    j += rands[:,i]
  end
end

const savedRands = randn(4,loopSize)
function test6()
  j=[0;0;0;0]
  for i = 1:loopSize
    j += savedRands[:,i]
  end
end

const chunkRands = ChunkedArray(randn,(4,),loopSize,parallel=false)
function test7()
  j=[0;0;0;0]
  for i = 1:loopSize
    j += next(chunkRands)
  end
end

function test8()
  rands2 = ChunkedArray(randn,(4,),buffSize,parallel=false)
  j=[0;0;0;0]
  for i = 1:loopSize
    j += next(rands2)
  end
end

function test9()
  rands3 = ChunkedArray(randn,(4,),loopSize,parallel=false)
  j=[0;0;0;0]
  for i = 1:loopSize
    j += rands3.randBuffer[i]
  end
end

function test10()
  rands4 = ChunkedArray(randn,(4,),loopSize,parallel=false)
  j=[0;0;0;0]
  for i = 1:loopSize
    j += next(rands4)
  end
end

function test11()
  rands5 = ChunkedArray(randn,(4,),buffSize,parallel=true)
  j=[0;0;0;0]
  for i = 1:loopSize
    j += next(rands5)
  end
end

t1 = @elapsed(for i=1:numRuns test1() end)/numRuns
t2 = @elapsed(for i=1:numRuns test2() end)/numRuns
t3 = @elapsed(for i=1:numRuns test3() end)/numRuns
t4 = @elapsed(for i=1:numRuns test4() end)/numRuns
t5 = @elapsed(for i=1:numRuns test5() end)/numRuns
t6 = @elapsed(for i=1:numRuns test6() end)/numRuns
t7 = @elapsed(for i=1:numRuns test7() end)/numRuns
t8 = @elapsed(for i=1:numRuns test8() end)/numRuns
t9 = @elapsed(for i=1:numRuns test9() end)/numRuns
t10= @elapsed(for i=1:numRuns test10() end)/numRuns
t11= @elapsed(for i=1:numRuns test11() end)/numRuns

println("""Test Results For Average Time:
One-by-one:                             $t1
Thousand-by-Thousand:                   $t2
Altogether:                             $t3
Hundred-by-hundred:                     $t4
Take at Beginning:                      $t5
Pre-made Rands:                         $t6
Chunked Rands Premade:                  $t7
Chunked Rands $buffSize buffer:              $t8
Chunked Rands Direct:                   $t9
Chunked Rands Max buffer:               $t10
Parallel Chunked Rands $buffSize buffer:     $t11
""")
toc()


text:
#!/usr/bin/env julia

export pe7

# need Primes.jl
# https://github.com/JuliaMath/Primes.jl
using Primes

pe7()=prime(10001)

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe7()|>print
end

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_name: miRoox/ProjectEuler


text:
# initialize.jl

using CUDA
using OffsetArrays
using TimerOutputs

""" Initialize the state vector of the model."""
function initModelState(grid, h0, u0, v0, w0, useCUDA, myFloat)
    # Allocate model state
    # with 4 parameters
    # for now cell-centered,
    # i.e. in all direction the same number of degrees of freedom
    nx = grid.nx
    ny = grid.ny
    nz = grid.nz

    # Generate initial arrays with start index 0
    useOffset = true

    h    = initField(h0, nx, ny, nz, useCUDA, useOffset, myFloat)
    u    = initField(u0, nx, ny, nz, useCUDA, useOffset, myFloat)
    v    = initField(v0, nx, ny, nz, useCUDA, useOffset, myFloat)
    w    = initField(w0, nx, ny, nz, useCUDA, useOffset, myFloat)

    hⁿ⁺¹ = initField(h0, nx, ny, nz, useCUDA, useOffset, myFloat)
    uⁿ⁺¹ = initField(u0, nx, ny, nz, useCUDA, useOffset, myFloat)
    vⁿ⁺¹ = initField(v0, nx, ny, nz, useCUDA, useOffset, myFloat)
    wⁿ⁺¹ = initField(w0, nx, ny, nz, useCUDA, useOffset, myFloat)

    state = State(h, u, v, w, hⁿ⁺¹, uⁿ⁺¹, vⁿ⁺¹, wⁿ⁺¹)
    return state
end

""" Instantiate a simulation object (without timer output)."""
function initSimulation(modelInput, useCUDA, myFloat)
    initSimulation(modelInput, useCUDA, myFloat, TimerOutput())
end

""" Instantiate a simulation object (with timer output)."""
function initSimulation(modelInput, useCUDA, myFloat, to)

    @synctimeit to "initialization" begin
        # Initialize the correct sizes/dimensions of the model
        nx        = modelInput.nx
        ny        = modelInput.ny
        nz        = modelInput.nz
        Δx        = modelInput.Δx
        Δy        = modelInput.Δy
        Δz        = modelInput.Δz
        Δt        = modelInput.Δt
        tend      = modelInput.tend
        K0        = modelInput.K0
        S0        = modelInput.S0
        h0        = modelInput.h0
        u0        = modelInput.u0
        v0        = modelInput.v0
        w0        = modelInput.w0
        source    = modelInput.source
        i_src     = modelInput.i_src
        j_src     = modelInput.j_src
        k_src     = modelInput.k_src
        duration  = modelInput.duration
        ΔhConv    = modelInput.ΔhConv
        constRecharge = modelInput.constRecharge
        rechargeFactor = modelInput.rechargeFactor
        boundaryPressure = modelInput.boundaryPressure

        useOffset = true
        noOffset  = false

        # The grid
        x, y, z    = gridCoords(nx, ny, nz, Δx, Δy, Δz, useCUDA, useOffset, myFloat)
        AT         = typeof(x)
        grid       = Grid{AT,myFloat}(nx, ny, nz, Δx, Δy, Δz, x, y, z, useCUDA ? GPUArch() : CPUArch())

        # State vector
        state      = initModelState(grid, h0, u0, v0, w0, useCUDA, myFloat)

        # External forcing
        externals  = initField(0.0, nx, ny, nz, useCUDA, noOffset, myFloat)

        # Model parameters
        K          = initField(K0, nx, ny, nz, useCUDA, useOffset, myFloat)
        specificStorage = S0

        # Sources/sinks
        AT = typeof(externals)
        source     = Source{AT,myFloat}(i_src, j_src, k_src, duration, source, externals)
        setSources!(0.0, source)

        # Recharge
        recharge = Recharge{myFloat}(constRecharge, 0.0, rechargeFactor)

        # Store the boundary conditions
        if (useCUDA)
            boundaryPressure = CuArray(boundaryPressure)
        end
        boundaryConditions = BoundaryConditions(boundaryPressure)

        # Group some parameters in the model.
        # For now sources and boundary conditions
        model      = Model(source, recharge, boundaryConditions)

        # Initialize the set of parameters (for now only K)
        AT = typeof(K)
        parameters = Parameters{AT,myFloat}(K, specificStorage)

        # Input object is no longer needed
        modelInput = nothing

        # Time related data
        maxsteps  = round(Int64, tend / Δt)
        time      = myFloat(0.0)
        timeData  = TimeData(Δt, tend, time, maxsteps)

        # Solver data
        Δh = initField(0.0, nx, ny, nz, useCUDA, useOffset, myFloat)
        AT = typeof(Δh)
        solverData = SolverData{AT,myFloat}(ΔhConv, Δh)

        return Simulation(grid, model, state, parameters, timeData, solverData)
    end # end timer
end


text:
export ParticleCollisions

struct ParticleCollisions

    dt

    function ParticleCollisions( p :: Particles)

        dt = zeros(p.n, p.n)
        fill!(dt, Inf)

        for k in 1:p.n
            for l in (k+1):p.n
                dt[k, l] = compute_dt(p, l, k)
            end
        end
        new( dt )
    end

end

function reset!(pc :: ParticleCollisions, i)
    pc.dt[i, :] .= Inf
    pc.dt[:, i] .= Inf
end

export BoxCollisions

struct BoxCollisions

    dt

    function BoxCollisions( p :: HardSpheres)

        dt = zeros(p.n, 2)
        fill!(dt, Inf)

        for k in 1:p.n
            t, m = compute_dt(p, k)
            dt[k, m] = t
        end

        new(dt)

    end

    function BoxCollisions( p :: SquareParticles )

        dt = zeros(p.n, 4)
        fill!(dt, Inf)

        for k in 1:p.n
            t, m = compute_dt(p, k)
            dt[k, m] = t
        end

        new(dt)

    end

end

function reset!(bc :: BoxCollisions, i)
    bc.dt[i, :] .= Inf
end

function compute_dt!(pc :: ParticleCollisions, i, particles :: Particles)

    for j in 1:particles.n
        if i != j
            pc.dt[i, j] = compute_dt(particles, j, i)
        end
    end
end

function dt_min_position( collisions )
    p = argmin(collisions.dt)
    dt_min = collisions.dt[p]
    return dt_min, p[1], p[2]
end


text:
include("../imports.jl")

@testset ExtendedTestSet "projective" begin
    image = Image(rand(RGB, 32, 32))
    tfm = Project(Translation(20, 10))

    @test_nowarn apply(tfm, image)
    timage = apply(tfm, image)
    @test timage.bounds.rs == (21:52, 11:42)

    @test_nowarn apply!(timage, tfm, image)

end


text:
"""
    Frames

Stores the actual computed frames and the user input
which can be `:same` or `Rel(10)`.
The `frames` are computed in `javis`.
"""
mutable struct Frames{T}
    frames::Union{Nothing,UnitRange}
    user::T
end

Base.convert(::Type{Frames}, x::Union{Symbol,Rel}) = Frames(nothing, x)
Base.convert(::Type{Frames}, x::UnitRange) = Frames(x, x)

"""
    set_frames!(a::AbstractAction, last_frames::UnitRange)

Compute the frames based on a.frames and `last_frames`.
Save the result in `a.frames.frames` which can be accessed via [`get_frames`](@ref).
"""
function set_frames!(a::AbstractAction, last_frames::UnitRange)
    frames = a.frames.user
    a.frames.frames = get_frames(frames, last_frames)
end

"""
    get_frames(a::AbstractAction)

Return `a.frames.frames` which holds the computed frames for the AbstractAction `a`.
"""
get_frames(a::AbstractAction) = a.frames.frames

"""
    get_frames(frames::Symbol, last_frames::UnitRange)

Get the frames based on a symbol (currently only `same`) and the `last_frames`.
Throw `ArgumentError` if symbol is unknown
"""
function get_frames(frames::Symbol, last_frames::UnitRange)
    if frames === :same
        return last_frames
    else
        throw(ArgumentError("Currently the only symbol supported for defining frames is `:same`"))
    end
end

"""
    get_frames(frames::Rel, last_frames::UnitRange)

Return the frames based on a relative frames [`Rel`](@ref) object and the `last_frames`.
"""
function get_frames(frames::Rel, last_frames::UnitRange)
    start_frame = last(last_frames) + first(frames.rel)
    last_frame = last(last_frames) + last(frames.rel)
    return start_frame:last_frame
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule RadeonProRender_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("RadeonProRender")
JLLWrappers.@generate_main_file("RadeonProRender", UUID("e92fa02b-1b7a-5e4e-a5cb-14fe26dfc45b"))
end  # module RadeonProRender_jll


text:
using Base.Test

rqa = Rqa.rqadata(collect("rqa"))
@test_approx_eq Rqa.RR(rqa) 0
@test Rqa.DET(rqa) |> isnan
@test Rqa.LMEAN(rqa) |> isnan
@test Rqa.LENT(rqa) |> isnan

alabama = Rqa.rqadata(collect("alabama"))
@test_approx_eq Rqa.RR(alabama) 6/21
@test_approx_eq Rqa.DET(alabama) 0
@test Rqa.LMEAN(alabama) |> isnan
@test Rqa.LENT(alabama) |> isnan

zzz = Rqa.rqadata(collect("zzz"))
@test_approx_eq Rqa.RR(zzz) 1
@test_approx_eq Rqa.DET(zzz) 2/3
@test_approx_eq Rqa.LMEAN(zzz) 2
@test_approx_eq Rqa.LENT(zzz) 0

w7 = Rqa.crqadata(collect("www"), collect("wwww"))
@test_approx_eq Rqa.RR(w7) 1
@test_approx_eq Rqa.DET(w7) 10/12
@test_approx_eq Rqa.LMEAN(w7) 2.5
@test_approx_eq Rqa.LENT(w7) 1

taran = Rqa.crqadata(collect("ratatar"), collect("taran"))
@test_approx_eq Rqa.RR(taran) 10/35
@test_approx_eq Rqa.DET(taran) 7/10
@test_approx_eq Rqa.LMEAN(taran) 7/3
@test_approx_eq Rqa.LENT(taran) (-(1/3) * log2(1/3) - (2/3) * log2(2/3))
@test_approx_eq Rqa.DET(taran, 3) 3/10
@test_approx_eq Rqa.LMEAN(taran, 3) 3
@test_approx_eq Rqa.LENT(taran, 3) 0
@test_approx_eq Rqa.DET(taran, 4) 0
@test Rqa.LMEAN(taran, 4) |> isnan
@test Rqa.LENT(taran, 4) |> isnan


text:
using LinearAlgebra
using Plots

temps = [ 71, 71, 68, 69, 68, 69, 68, 74, 77, 82, 85, 86,
 88, 86, 85, 86, 84, 79, 77, 75, 73, 71, 70, 70, 69, 69, 69,
 69, 67, 68, 68, 73, 76, 77, 82, 84, 84, 81, 80, 78, 79, 78,
 73, 72, 70, 70, 68, 67 ];

plot(temps, marker = :circle, legend = false, grid = false)
#savefig(temps,"temperature.pdf")

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
module InvertedIndices

export InvertedIndex, Not

using Base: tail

"""
    InvertedIndex(idx)
    Not(idx)

Construct an inverted index, selecting all indices not in the passed `idx`.

Upon indexing into an array, the `InvertedIndex` behaves like a 1-dimensional
collection of the indices of the array that are not in `idx`. Bounds are
checked to ensure that all indices in `idx` are within the bounds of the array
— even though they are skipped. The `InvertedIndex` behaves like a
1-dimensional collection of its inverted indices. If `idx` spans multiple
dimensions (like a multidimensional logical mask or `CartesianIndex`), then the
inverted index will similarly span multiple dimensions.
"""
struct InvertedIndex{S}
    skip::S
end
const Not = InvertedIndex

# Like Base.LogicalIndex, the InvertedIndexIterator is a pseudo-vector that is
# just used as an iterator and does not support getindex.
struct InvertedIndexIterator{T,S,P} <: AbstractVector{T}
    skips::S
    picks::P
end
InvertedIndexIterator(skips, picks) = InvertedIndexIterator{eltype(picks), typeof(skips), typeof(picks)}(skips, picks)
Base.size(III::InvertedIndexIterator) = (length(III.picks) - length(III.skips),)

@inline Base.iterate(I::InvertedIndexIterator) = iterate(I, (iterate(I.skips), iterate(I.picks)))
Base.iterate(I::InvertedIndexIterator, ::Tuple{Any, Nothing}) = nothing
@inline function Base.iterate(I::InvertedIndexIterator, (skipitr, pickitr))
    while should_skip(skipitr, pickitr)
        skipitr = iterate(I.skips, tail(skipitr)...)
        pickitr = iterate(I.picks, tail(pickitr)...)
        pickitr === nothing && return nothing
    end
    return (pickitr[1], (skipitr, iterate(I.picks, tail(pickitr)...)))
end
Base.collect(III::InvertedIndexIterator) = [i for i in III]

should_skip(::Nothing, ::Any) = false
should_skip(s::Tuple, p::Tuple) = _should_skip(s[1], p[1])
_should_skip(s, p) = s == p
_should_skip(s::Integer, p::CartesianIndex{1}) = s == p.I[1]
_should_skip(s::CartesianIndex{1}, p::Integer) = s.I[1] == p

@inline Base.checkbounds(::Type{Bool}, A::AbstractArray, I::InvertedIndexIterator{<:Integer}) =
    checkbounds(Bool, A, I.skips) && eachindex(IndexLinear(), A) == eachindex(IndexLinear(), I.picks)
@inline Base.checkbounds(::Type{Bool}, A::AbstractArray, I::InvertedIndexIterator) = checkbounds(Bool, A, I.skips) && axes(A) == axes(I.picks)
@inline Base.checkindex(::Type{Bool}, indx::AbstractUnitRange, I::InvertedIndexIterator) = checkindex(Bool, indx, I.skips) && (indx,) == axes(I.picks)
@inline Base.checkindex(::Type{Bool}, inds::Tuple, I::InvertedIndexIterator) = checkindex(Bool, indx, I.skips) && inds == axes(I.picks)

@inline Base.ensure_indexable(I::Tuple{InvertedIndexIterator, Vararg{Any}}) = (collect(I[1]), Base.ensure_indexable(tail(I))...)

Base.show(io::IO, I::InvertedIndexIterator) = show(io, collect(I))

# We use a custom sort method that supports Base.LogicalIndex and CartesianIndex
sort(A::AbstractArray) = Base.sort(A)
sort(A::Union{Number,CartesianIndex,Base.LogicalIndex}) = A

@inline function Base.to_indices(A, inds, I::Tuple{InvertedIndex, Vararg{Any}})
    new_indices = to_indices(A, inds, (I[1].skip, tail(I)...))
    skips = sort(new_indices[1])
    picks = spanned_indices(inds, skips)[1]
    return (InvertedIndexIterator(skips, picks), tail(new_indices)...)
end

struct ZeroDArray{T} <: AbstractArray{T,0}
    x::T
end
Base.size(::ZeroDArray) = ()
Base.getindex(Z::ZeroDArray) = Z.x

# Be careful with CartesianIndex as they splat out to a variable number of new indices and do not iterate
function Base.to_indices(A, inds, I::Tuple{InvertedIndex{<:CartesianIndex}, Vararg{Any}})
    new_indices = to_indices(A, inds, (I[1].skip, tail(I)...))
    skips = ZeroDArray(I[1].skip)
    picks, tails = spanned_indices(inds, skips)
    return (InvertedIndexIterator(skips, picks), to_indices(A, tails, tail(I))...)
end

# Either return a CartesianRange or an axis vector
@inline spanned_indices(inds, ::Any) = inds[1], tail(inds)
const NIdx{N} = Union{CartesianIndex{N}, AbstractArray{CartesianIndex{N}}, AbstractArray{Bool,N}}
@inline spanned_indices(inds, ::NIdx{0}) = CartesianIndices(()), inds
@inline spanned_indices(inds, ::NIdx{1}) = inds[1], tail(inds)
@inline function spanned_indices(inds, ::NIdx{N}) where N
    heads, tails = Base.IteratorsMD.split(inds, Val(N))
    return CartesianIndices(heads), tails
end

# It's possible more indices were specified than there were axes
@inline spanned_indices(inds::Tuple{}, ::Any) = Base.OneTo(1), ()
@inline spanned_indices(inds::Tuple{}, ::NIdx{0}) = CartesianIndices(()), ()
@inline spanned_indices(inds::Tuple{}, ::NIdx{1}) = Base.OneTo(1), ()
@inline spanned_indices(inds::Tuple{}, ::NIdx{N}) where N = CartesianIndices(ntuple(i->Base.OneTo(1), N)), ()

# This is an interesting need — we need this because otherwise indexing with a
# single multidimensional boolean array ends up comparing a multidimensional cartesian
# index to a linear index. Does this need addressing in Base, too?
@inline Base.to_indices(A, I::Tuple{Not{<:NIdx{1}}}) = to_indices(A, (eachindex(IndexLinear(), A),), I)
@inline Base.to_indices(A, I::Tuple{Not{<:NIdx}}) = to_indices(A, axes(A), I)
# Arrays of Bool are even more confusing as they're sometimes linear and sometimes not
@inline Base.to_indices(A, I::Tuple{Not{<:AbstractArray{Bool, 1}}}) = to_indices(A, (eachindex(IndexLinear(), A),), I)
@inline Base.to_indices(A, I::Tuple{Not{<:Union{Array{Bool}, BitArray}}}) = to_indices(A, (eachindex(A),), I)

end # module


text:
using Weave
using Base.Test

include("chunk_options.jl")
#include("pyplot_formats.jl") cause segfaults, but OK if run by itself

#Test Gadfly and markdown reader, Gadfly only works with 0.3
if VERSION.minor == 3
    weave("documents/gadfly_markdown_test.jmd", doctype="github",plotlib="gadfly", informat="markdown")
    result = readall(open("documents/gadfly_markdown_test.md"))
    ref = readall(open("documents/gadfly_markdown_test_ref.md"))
    @test result == ref

    include("winston_formats.jl")
    include("gadfly_formats.jl")
end


text:
function get_fields
end

function get_records
end

function build_result
end

# interface for `DataFrame`s
get_fields(df :: DataFrame) = names(df)
get_records(df :: DataFrame) = zip(DataFrames.columns(df)...)
function build_result(::Type{DataFrame}, fields, typs, source :: Base.Generator)
    res = Tuple(typ[] for typ in typs)
    for each in source
        push!.(res, each)
    end
    DataFrame(collect(res), fields)
end

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: cherichy/MLStyle-Playground

max_forks_repo_name: cherichy/MLStyle-Playground


text:
## ===== Interface for an AbstractDFG =====

"""
    $(SIGNATURES)

De-serialization of IncrementalInference objects require discovery of foreign types.

Example:

Template to tunnel types from a user module:
```julia
# or more generic solution -- will always try Main if available
IIF.setSerializationNamespace!(Main)

# or a specific package such as RoME if you import all variable and factor types into a specific module.
using RoME
IIF.setSerializationNamespace!(RoME)
```
"""
function setSerializationModule!(dfg::G, mod::Module)::Nothing where G <: AbstractDFG
    @warn "Setting serialization module from AbstractDFG - override this in the '$(typeof(dfg)) structure! This is being ignored."
end

function getSerializationModule(dfg::G)::Module where G <: AbstractDFG
    @warn "Retrieving serialization module from AbstractDFG - override this in the '$(typeof(dfg)) structure! This is returning Main"
    return Main
end

# Accessors
function getLabelDict(dfg::G) where G <: AbstractDFG
	error("getLabelDict not implemented for $(typeof(dfg))")
end
function getDescription(dfg::G) where G <: AbstractDFG
	error("getDescription not implemented for $(typeof(dfg))")
end
function setDescription(dfg::G, description::String) where G <: AbstractDFG
	error("setDescription not implemented for $(typeof(dfg))")
end
function getAddHistory(dfg::G) where G <: AbstractDFG
	error("getAddHistory not implemented for $(typeof(dfg))")
end
function getSolverParams(dfg::G) where G <: AbstractDFG
	error("getSolverParams not implemented for $(typeof(dfg))")
end
function setSolverParams(dfg::G, solverParams::T) where {G <: AbstractDFG, T <: AbstractParams}
	error("setSolverParams not implemented for $(typeof(dfg))")
end

# Get user, robot, and session "small" data.
# function getUserData(dfg::G)::Dict{Symbol, String} where {G <: AbstractDFG}
# 	error("getUserData not implemented for $(typeof(dfg))")
# end
# function setUserData(dfg::G, data::Dict{Symbol, String})::Bool where {G <: AbstractDFG}
# 	error("setUserData not implemented for $(typeof(dfg))")
# end
# function getRobotData(dfg::G)::Dict{Symbol, String} where {G <: AbstractDFG}
# 	error("getRobotData not implemented for $(typeof(dfg))")
# end
# function setRobotData(dfg::G, data::Dict{Symbol, String})::Bool where {G <: AbstractDFG}
# 	error("setRobotData not implemented for $(typeof(dfg))")
# end
# function getSessionData(dfg::G)::Dict{Symbol, String} where {G <: AbstractDFG}
# 	error("getSessionData not implemented for $(typeof(dfg))")
# end
# function setSessionData(dfg::G, data::Dict{Symbol, String})::Bool where {G <: AbstractDFG}
# 	error("setSessionData not implemented for $(typeof(dfg))")
# end

getUserData(dfg::AbstractDFG)::Dict{Symbol, String} = return dfg.userData
function setUserData(dfg::AbstractDFG, data::Dict{Symbol, String})::Bool
	dfg.userData = data
	return true
end
getRobotData(dfg::AbstractDFG)::Dict{Symbol, String} = return dfg.robotData
function setRobotData(dfg::AbstractDFG, data::Dict{Symbol, String})::Bool
	dfg.robotData = data
	return true
end
getSessionData(dfg::AbstractDFG)::Dict{Symbol, String} = return dfg.sessionData
function setSessionData(dfg::AbstractDFG, data::Dict{Symbol, String})::Bool
	dfg.sessionData = data
	return true
end

pushUserData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.userData, pair)
pushRobotData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.userData, pair)
pushSessionData!(dfg::AbstractDFG, pair::Pair{Symbol,String}) = push!(dfg.userData, pair)

popUserData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.userData, key)
popRobotData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.userData, key)
popSessionData!(dfg::AbstractDFG, key::Symbol) = pop!(dfg.userData, key)

"""
    $(SIGNATURES)
True if the variable or factor exists in the graph.
"""
function exists(dfg::G, node::N) where {G <: AbstractDFG, N <: DFGNode}
	error("exists not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Add a DFGVariable to a DFG.
"""
function addVariable!(dfg::G, variable::V)::Bool where {G <: AbstractDFG, V <: AbstractDFGVariable}
	error("addVariable! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::G, variables::Vector{V}, factor::F)::Bool where {G <: AbstractDFG, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
	error("addFactor! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::G, variableIds::Vector{Symbol}, factor::F)::Bool where {G <: AbstractDFG, F <: AbstractDFGFactor}
	error("addFactor! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its underlying integer ID.
"""
function getVariable(dfg::G, variableId::Int64)::AbstractDFGVariable where G <: AbstractDFG
	error("getVariable not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its label.
"""
function getVariable(dfg::G, label::Union{Symbol, String})::AbstractDFGVariable where G <: AbstractDFG
	return getVariable(dfg, Symbol(label))
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its underlying integer ID.
"""
function getFactor(dfg::G, factorId::Int64)::AbstractDFGFactor where G <: AbstractDFG
	error("getFactor not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its label.
"""
function getFactor(dfg::G, label::Union{Symbol, String})::AbstractDFGFactor where G <: AbstractDFG
	return getFactor(dfg, Symbol(label))
end

"""
    $(SIGNATURES)
Update a complete DFGVariable in the DFG.
"""
function updateVariable!(dfg::G, variable::V)::AbstractDFGVariable where {G <: AbstractDFG, V <: AbstractDFGVariable}
	error("updateVariable! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG.
"""
function updateFactor!(dfg::G, factor::F)::AbstractDFGFactor where {G <: AbstractDFG, F <: AbstractDFGFactor}
	error("updateFactor! not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Delete a DFGVariable from the DFG using its label.
"""
function deleteVariable!(dfg::G, label::Symbol)::AbstractDFGVariable where G <: AbstractDFG
	error("deleteVariable! not implemented for $(typeof(dfg))")
end

#Alias
"""
    $(SIGNATURES)
Delete a referenced DFGVariable from the DFG.
"""
function deleteVariable!(dfg::G, variable::V)::AbstractDFGVariable where {G <: AbstractDFG, V <: AbstractDFGVariable}
	return deleteVariable!(dfg, variable.label)
end

"""
    $(SIGNATURES)
Delete a DFGFactor from the DFG using its label.
"""
function deleteFactor!(dfg::G, label::Symbol)::AbstractDFGFactor where G <: AbstractDFG
	error("deleteFactors not implemented for $(typeof(dfg))")
end

# Alias
"""
    $(SIGNATURES)
Delete the referened DFGFactor from the DFG.
"""
function deleteFactor!(dfg::G, factor::F)::AbstractDFGFactor where {G <: AbstractDFG, F <: AbstractDFGFactor}
	return deleteFactor!(dfg, factor.label)
end

"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function getVariables(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{AbstractDFGVariable} where G <: AbstractDFG
	error("getVariables not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a list of IDs of the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.

Example
```julia
getVariableIds(dfg, r"l", tags=[:APRILTAG;])
```

Related

ls
"""
function getVariableIds(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{Symbol} where G <: AbstractDFG
  vars = getVariables(dfg, regexFilter, tags=tags)
  # mask = map(v -> length(intersect(v.tags, tags)) > 0, vars )
  return map(v -> v.label, vars)
end

# Alias
"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function ls(dfg::G, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{Symbol} where G <: AbstractDFG
	return getVariableIds(dfg, regexFilter, tags=tags)
end

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactors(dfg::G, regexFilter::Union{Nothing, Regex}=nothing)::Vector{AbstractDFGFactor} where G <: AbstractDFG
	error("getFactors not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a list of the IDs of the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactorIds(dfg::G, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} where G <: AbstractDFG
	return map(f -> f.label, getFactors(dfg, regexFilter))
end

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
# Alias
function lsf(dfg::G, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol} where G <: AbstractDFG
	return getFactorIds(dfg, regexFilter)
end

"""
	$(SIGNATURES)
Alias for getNeighbors - returns neighbors around a given node label.
"""
function lsf(dfg::G, label::Symbol)::Vector{Symbol} where G <: AbstractDFG
  return getNeighbors(dfg, label)
end

"""
    $(SIGNATURES)
Checks if the graph is fully connected, returns true if so.
"""
function isFullyConnected(dfg::G)::Bool where G <: AbstractDFG
	error("isFullyConnected not implemented for $(typeof(dfg))")
end

#Alias
"""
    $(SIGNATURES)
Checks if the graph is not fully connected, returns true if it is not contiguous.
"""
function hasOrphans(dfg::G)::Bool where G <: AbstractDFG
	return !isFullyConnected(dfg)
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function getNeighbors(dfg::G, node::T; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}  where {G <: AbstractDFG, T <: DFGNode}
	error("getNeighbors not implemented for $(typeof(dfg))")
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function getNeighbors(dfg::G, label::Symbol; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol} where G <: AbstractDFG
	error("getNeighbors not implemented for $(typeof(dfg))")
end

# Aliases
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function ls(dfg::G, node::T)::Vector{Symbol} where {G <: AbstractDFG, T <: DFGNode}
    return getNeighbors(dfg, node)
end
"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function ls(dfg::G, label::Symbol)::Vector{Symbol} where G <: AbstractDFG
    return getNeighbors(dfg, label)
end

"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::G)::G where G <: AbstractDFG
	error("_getDuplicatedEmptyDFG not implemented for $(typeof(dfg))")
end

# TODO: NEED TO FIGURE OUT SIGNATURE FOR DEFAULT ARGS

"""
    $(SIGNATURES)
Retrieve a deep subgraph copy around a given variable or factor.
Optionally provide a distance to specify the number of edges should be followed.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraphAroundNode(dfg::G, node::T, distance::Int64=1, includeOrphanFactors::Bool=false, addToDFG::H=_getDuplicatedEmptyDFG(dfg))::G where {G <: AbstractDFG, H <: AbstractDFG, P <: AbstractParams, T <: DFGNode}
	error("getSubgraphAroundNode not implemented for $(typeof(dfg))")
end

"""
    $(SIGNATURES)
Get a deep subgraph copy from the DFG given a list of variables and factors.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraph(dfg::G,
					 variableFactorLabels::Vector{Symbol},
					 includeOrphanFactors::Bool=false,
					 addToDFG::H=_getDuplicatedEmptyDFG(dfg))::H where {G <: AbstractDFG, H <: AbstractDFG}
    for label in variableFactorLabels
        if !exists(dfg, label)
            error("Variable/factor with label '$(label)' does not exist in the factor graph")
        end
    end

    _copyIntoGraph!(dfg, addToDFG, variableFactorLabels, includeOrphanFactors)
    return addToDFG
end

"""
    $(SIGNATURES)
Common function for copying nodes from one graph into another graph.
This is overridden in specialized implementations for performance.
"""
function _copyIntoGraph!(sourceDFG::G, destDFG::H, variableFactorLabels::Vector{Symbol}, includeOrphanFactors::Bool=false)::Nothing where {G <: AbstractDFG, H <: AbstractDFG}
    # Split into variables and factors
    sourceVariables = map(vId->getVariable(sourceDFG, vId), intersect(getVariableIds(sourceDFG), variableFactorLabels))
    sourceFactors = map(fId->getFactor(sourceDFG, fId), intersect(getFactorIds(sourceDFG), variableFactorLabels))
    if length(sourceVariables) + length(sourceFactors) != length(variableFactorLabels)
        rem = symdiff(map(v->v.label, sourceVariables), variableFactorLabels)
        rem = symdiff(map(f->f.label, sourceFactors), variableFactorLabels)
        error("Cannot copy because cannot find the following nodes in the source graph: $rem")
    end

    # Now we have to add all variables first,
    for variable in sourceVariables
        if !exists(destDFG, variable)
            addVariable!(destDFG, deepcopy(variable))
        end
    end
    # And then all factors to the destDFG.
    for factor in sourceFactors
        # Get the original factor variables (we need them to create it)
        sourceFactorVariableIds = getNeighbors(sourceDFG, factor)
        # Find the labels and associated variables in our new subgraph
        factVariableIds = Symbol[]
        for variable in sourceFactorVariableIds
            if exists(destDFG, variable)
                push!(factVariableIds, variable)
            end
        end
        # Only if we have all of them should we add it (otherwise strange things may happen on evaluation)
        if includeOrphanFactors || length(factVariableIds) == length(sourceFactorVariableIds)
            if !exists(destDFG, factor)
                addFactor!(destDFG, factVariableIds, deepcopy(factor))
            end
        end
    end
    return nothing
end

"""
    $(SIGNATURES)
Merges and updates solver and estimate data for a variable (variable can be from another graph).
Note: Makes a copy of the estimates and solver data so that there is no coupling
between graphs.
"""
function mergeUpdateVariableSolverData!(dfg::AbstractDFG, sourceVariable::AbstractDFGVariable)::AbstractDFGVariable
    if !exists(dfg, sourceVariable)
        error("Source variable '$(sourceVariable.label)' doesn't exist in the graph.")
    end
    var = getVariable(dfg, sourceVariable.label)
    # We don't know which graph this came from, must be copied!
    merge!(var.estimateDict, deepcopy(sourceVariable.estimateDict))
	# If this variable has solverDataDict (summaries do not)
	:solverDataDict in fieldnames(typeof(var)) && merge!(var.solverDataDict, deepcopy(sourceVariable.solverDataDict))
    return sourceVariable
end

"""
    $(SIGNATURES)
Common function to update all solver data and estimates from one graph to another.
This should be used to push local solve data back into a cloud graph, for example.
"""
function mergeUpdateGraphSolverData!(sourceDFG::G, destDFG::H, varSyms::Vector{Symbol})::Nothing where {G <: AbstractDFG, H <: AbstractDFG}
    # Update all variables in the destination
    # (For now... we may change this soon)
    for variableId in varSyms
        mergeUpdateVariableSolverData!(destDFG, getVariable(sourceDFG, variableId))
    end
end

"""
    $(SIGNATURES)
Get an adjacency matrix for the DFG, returned as a Matrix{Union{Nothing, Symbol}}.
Rows are all factors, columns are all variables, and each cell contains either nothing or the symbol of the relating factor.
The first row and first column are factor and variable headings respectively.
"""
function getAdjacencyMatrix(dfg::G)::Matrix{Union{Nothing, Symbol}} where G <: AbstractDFG
	varLabels = sort(map(v->v.label, getVariables(dfg)))
    factLabels = sort(map(f->f.label, getFactors(dfg)))
    vDict = Dict(varLabels .=> [1:length(varLabels)...].+1)

    adjMat = Matrix{Union{Nothing, Symbol}}(nothing, length(factLabels)+1, length(varLabels)+1)
    # Set row/col headings
    adjMat[2:end, 1] = factLabels
    adjMat[1, 2:end] = varLabels
    for (fIndex, factLabel) in enumerate(factLabels)
        factVars = getNeighbors(dfg, getFactor(dfg, factLabel))
        map(vLabel -> adjMat[fIndex+1,vDict[vLabel]] = factLabel, factVars)
    end
    return adjMat
end

"""
    $(SIGNATURES)
Get an adjacency matrix for the DFG, returned as a tuple: adjmat::SparseMatrixCSC{Int}, var_labels::Vector{Symbol) fac_labels::Vector{Symbol).
Rows are the factors, columns are the variables, with the corresponding labels in fac_labels,var_labels.
"""
function getAdjacencyMatrixSparse(dfg::G)::Tuple{LightGraphs.SparseMatrixCSC, Vector{Symbol}, Vector{Symbol}} where G <: AbstractDFG
	varLabels = map(v->v.label, getVariables(dfg))
	factLabels = map(f->f.label, getFactors(dfg))

	vDict = Dict(varLabels .=> [1:length(varLabels)...])

	adjMat = spzeros(Int, length(factLabels), length(varLabels))

	for (fIndex, factLabel) in enumerate(factLabels)
		factVars = getNeighbors(dfg, getFactor(dfg, factLabel))
	    map(vLabel -> adjMat[fIndex,vDict[vLabel]] = 1, factVars)
	end
	return adjMat, varLabels, factLabels
end

"""
    $SIGNATURES

Return boolean whether a factor `label` is present in `<:AbstractDFG`.
"""
function hasFactor(dfg::G, label::Symbol)::Bool where {G <: AbstractDFG}
	@warn "hasFactor() deprecated, please use exists()"
	return exists(dfg, label)
end

"""
    $(SIGNATURES)

Return `::Bool` on whether `dfg` contains the variable `lbl::Symbol`.
"""
function hasVariable(dfg::G, label::Symbol)::Bool where {G <: AbstractDFG}
	@warn "hasVariable() deprecated, please use exists()"
	return exists(dfg, label) # haskey(vertices(dfg.g), label)
end

"""
    $SIGNATURES

Returns state of vertex data `.initialized` flag.

Notes:
- used by both factor graph variable and Bayes tree clique logic.
"""
function isInitialized(var::DFGVariable; key::Symbol=:default)::Bool
  	data = solverData(var, key)
  	if data == nothing
		@error "Variable does not have solver data $(key)"
		return false
  	else
  		return data.initialized
	end
end

function isInitialized(dfg::AbstractDFG, label::Symbol; key::Symbol=:default)::Bool
  return isInitialized(getVariable(dfg, label), key=key)
end

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a variable vertex in the graph.
"""
isVariable(dfg::G, sym::Symbol) where G <: AbstractDFG = hasVariable(dfg, sym)

"""
    $SIGNATURES

Return whether `sym::Symbol` represents a factor vertex in the graph.
"""
isFactor(dfg::G, sym::Symbol) where G <: AbstractDFG = hasFactor(dfg, sym)

"""
    $SIGNATURES

Return reference to the user factor in `<:AbstractDFG` identified by `::Symbol`.
"""
getFactorFunction(fcd::GenericFunctionNodeData) = fcd.fnc.usrfnc!
getFactorFunction(fc::DFGFactor) = getFactorFunction(solverData(fc))
function getFactorFunction(dfg::G, fsym::Symbol) where G <: AbstractDFG
  getFactorFunction(getFactor(dfg, fsym))
end

"""
    $SIGNATURES

Display and return to console the user factor identified by tag name.
"""
showFactor(fgl::G, fsym::Symbol) where G <: AbstractDFG = @show getFactor(fgl,fsym)

"""
    $(SIGNATURES)
Produces a dot-format of the graph for visualization.
"""
function toDot(dfg::AbstractDFG)::String
    #TODO implement convert
    graphsdfg = GraphsDFG{AbstractParams}()
    DistributedFactorGraphs._copyIntoGraph!(dfg, graphsdfg, union(getVariableIds(dfg), getFactorIds(dfg)), true)

	# Calls down to GraphsDFG.toDot
    return toDot(graphsdfg)
end

"""
    $(SIGNATURES)
Produces a dot file of the graph for visualization.
Download XDot to see the data

Note
- Default location "/tmp/dfg.dot" -- MIGHT BE REMOVED
- Can be viewed with the `xdot` system application.
- Based on graphviz.org
"""
function toDotFile(dfg::AbstractDFG, fileName::String="/tmp/dfg.dot")::Nothing
    #TODO implement convert
    graphsdfg = GraphsDFG{AbstractParams}()
    DistributedFactorGraphs._copyIntoGraph!(dfg, graphsdfg, union(getVariableIds(dfg), getFactorIds(dfg)), true)

    open(fileName, "w") do fid
        write(fid,Graphs.to_dot(graphsdfg.g))
    end
    return nothing
end

"""
    $(SIGNATURES)
Get a summary of the graph (first-class citizens of variables and factors).
Returns a AbstractDFGSummary.
"""
function getSummary(dfg::G)::AbstractDFGSummary where {G <: AbstractDFG}
	vars = map(v -> convert(DFGVariableSummary, v), getVariables(dfg))
	facts = map(f -> convert(DFGFactorSummary, f), getFactors(dfg))
	return AbstractDFGSummary(
		Dict(map(v->v.label, vars) .=> vars),
		Dict(map(f->f.label, facts) .=> facts),
		dfg.userId,
		dfg.robotId,
		dfg.sessionId)
end

"""
$(SIGNATURES)
Get a summary graph (first-class citizens of variables and factors) with the same structure as the original graph.
Note this is a copy of the original.
Returns a LightDFG{NoSolverParams, DFGVariableSummary, DFGFactorSummary}.
"""
function getSummaryGraph(dfg::G)::LightDFG{NoSolverParams, DFGVariableSummary, DFGFactorSummary} where {G <: AbstractDFG}
	summaryDfg = LightDFG{NoSolverParams, DFGVariableSummary, DFGFactorSummary}(
		description="Summary of $(dfg.description)",
		userId=dfg.userId,
		robotId=dfg.robotId,
		sessionId=dfg.sessionId)
	for v in getVariables(dfg)
		newV = addVariable!(summaryDfg, convert(DFGVariableSummary, v))
	end
	for f in getFactors(dfg)
		addFactor!(summaryDfg, getNeighbors(dfg, f), convert(DFGFactorSummary, f))
	end
	return summaryDfg
end


text:
@theme XcodeTheme Dict(
    :name => "Xcode",
    :description => "Style similar to the Xcode default colouring theme.",
    :comments => "Based on the Xcode theme from Pygments.",
    :tokens => Dict(
        COMMENT             => S"fg: 177500",
        COMMENT_PREPROC     => S"fg: 633820",

        ERROR               => S"fg: 000000",

        KEYWORD             => S"fg: A90D91",

        LITERAL             => S"fg: 1C01CE",

        NAME                => S"fg: 000000",
        NAME_ATTRIBUTE      => S"fg: 836C28",
        NAME_BUILTIN        => S"fg: A90D91",
        NAME_BUILTIN_PSEUDO => S"fg: 5B269A",
        NAME_CLASS          => S"fg: 3F6E75",
        NAME_DECORATOR      => S"fg: 000000",
        NAME_FUNCTION       => S"fg: 000000",
        NAME_LABEL          => S"fg: 000000",
        NAME_TAG            => S"fg: 000000",
        NAME_VARIABLE       => S"fg: 000000",

        NUMBER              => S"fg: 1C01CE",

        OPERATOR            => S"fg: 000000",

        STRING              => S"fg: C41A16",
        STRING_CHAR         => S"fg: 2300CE",

        TEXT                => S"",
    )
)


text:
# Author: Lukasz Bratos

module Diff
using Plots
gr()

export ilorazyRoznicowe, warNewton, naturalna, rysujNnfx

"""
Funkcja oblicza iloraz różnicowy

Dane:
x – wektor długości n + 1 zawierający węzły x0, ..., xn
    x[1] = x0, ..., x[n+1] = xn
f – wektor długości n + 1 zawierający wartości interpolowanej funkcji w węzłach
    f(x0), ..., f(xn)

Wyniki:
fx– wektor długości n + 1 zawierający obliczone ilorazy różnicowe
    fx[1] = f[x0],
    fx[2] = f[x0, x1], ..., fx[n] = f[x0, ..., xn−1], fx[n+1] = f[x0, ..., xn]
"""
function ilorazyRoznicowe(x::Vector{Float64}, f::Vector{Float64})
    len = length(f)
    fx = Vector{Float64}(undef, len)

    for i in 1:len
        fx[i] = f[i]
    end

    for i in 2:len
        for j in len:-1:i
            fx[j] = (fx[j] - fx[j - 1]) / (x[j] - x[j - i + 1])
        end
    end

    return fx
end

"""
Funkcja wyliczajaca wartosc wielomianu interpolacyjnego stopnia n w postaci Newtona Nn(x)
za pomoca uogolnionego algorytmu Hornera

Dane:
x – wektor długości n + 1 zawierający węzły x0, ..., xn
    x[1] = x0, ..., x[n+1] = xn
fx – wektor długości n+ 1 zawierający ilorazy różnicowe
    fx[1] = f[x0],
    fx[2] = f[x0, x1], ..., fx[n] = f[x0, ..., xn−1], fx[n+1] =f [x0, ..., xn]
t – punkt, w którym należy obliczyć wartość wielomianu

Wyniki:
nt – wartość wielomianu w punkcie t
"""
function warNewton(x::Vector{Float64}, fx::Vector{Float64}, t::Float64)
    len = length(x)
    nt = fx[len]

    for i in (len - 1):-1:1
        nt = fx[i] + (t - x[i]) * nt
    end

    return nt
end

"""
Funkcja wyliczajaca wspolczynniki wielomianu w postaci normanej z wielomianu
interpolacyjnego w postaci Newtona

Dane:
x – wektor długości n+ 1 zawierający węzły x0, ..., xn
    x[1] = x0, ..., x[n+1] = xn
fx – wektor długości n + 1 zawierający ilorazy różnicowe
    fx[1] = f[x0],
    fx[2] = f[x0, x1], ..., fx[n] = f[x0, ..., xn−1], fx[n+1] = f[x0, ..., xn]

Wyniki:
a – wektor długości n+ 1 zawierający obliczone współczynniki postaci naturalnej
    a[1] = a0,
    a[2] = a1, ..., a[n] = an−1, a[n+1] = an.
"""
function naturalna(x::Vector{Float64}, fx::Vector{Float64})
    len = length(x)
    a = Vector{Float64}(undef, len)
    a[len] = fx[len]

    for i in (len - 1):-1:1
        a[i] = fx[i] - a[i + 1] * x[i]

        for j in (i + 1):(len - 1)
            a[j] = a[j] - a[j + 1] * x[i]
        end
    end

    return a
end

"""
Funkcja rysująca wielomian interpolacyjny i interpolowaną funkcję na danym przedziale

Dane:
f – funkcja f(x) zadana jako anonimowa funkcja,
a, b – przedział interpolacji
n – stopień wielomianu interpolacyjnego
filename - nazwa pliku wynikowego

Wyniki:
– funkcja rysuje wielomian interpolacyjny i interpolowaną funkcję w przedziale [a, b]
"""
function rysujNnfx(f, a::Float64, b::Float64, n::Int, filename::String)
    dist = (b - a) / n

    x = Vector{Float64}(undef, n + 1)
    y = Vector{Float64}(undef, n + 1)

    for i in 1:(n + 1)
        x[i] = a + (i - 1) * dist
        y[i] = f(x[i]);
    end

    fx = ilorazyRoznicowe(x, y)

    numOfPoints = 100
    plot_dist = (b - a) / numOfPoints

    fx_val = Vector{Float64}(undef, numOfPoints + 1)
    wx_val = Vector{Float64}(undef, numOfPoints + 1)

    for i in 1:(numOfPoints + 1)
        t = a + (i - 1) * plot_dist
        fx_val[i] = f(t)
        wx_val[i] = warNewton(x, fx, t)
    end

    plot(range(a, stop=b, length=(numOfPoints + 1)), fx_val, color = "red", label = "f(x)")
    plot!(range(a, stop=b, length=(numOfPoints + 1)), wx_val, color = "blue", label = "w(x)")
    savefig("./plots/$filename.png")
end

end


text:
"""
    parse_robinx(filename::String)

Parses an instance given in the RobinX format.

# Arguments
- `filename::String` : The project path to the file to be parsed

# Returns
- `Instance` : the problem instance
"""
function parse_robinx(filename::String)
    xdoc = XML.parse_file(filename)
    root = XML.root(xdoc)

    formatNode = root["Structure"][1]["Format"][1]
    resourceNode = root["Resources"][1]
    constraintNode = root["Constraints"][1]

    nteams = length(resourceNode["Teams"][1]["team"])
    nslots = 2nteams - 2
    isPhased = XML.content(only(formatNode["gameMode"])) == "P"

    rr2 = RR2(nteams, nslots, isPhased)

    # capacity constraints
    constrCA = Vector{CAcon}(undef, 0)
    parse_ca!(rr2, constrCA, constraintNode["CapacityConstraints"][1]["CA1"])
    parse_ca!(rr2, constrCA, constraintNode["CapacityConstraints"][1]["CA2"])
    parse_ca!(rr2, constrCA, constraintNode["CapacityConstraints"][1]["CA3"])
    parse_ca!(rr2, constrCA, constraintNode["CapacityConstraints"][1]["CA4"])

    # break constraints
    constrBR = Vector{BRcon}(undef, 0)
    parse_br!(constrBR, constraintNode["BreakConstraints"][1]["BR1"])
    parse_br!(constrBR, constraintNode["BreakConstraints"][1]["BR2"])

    # game constraints
    constrGA = Vector{GAcon}(undef, 0)
    parse_ga!(constrGA, constraintNode["GameConstraints"][1]["GA1"])

    # fairness constraints
    constrFA = Vector{FAcon}(undef, 0)
    parse_fa!(constrFA, constraintNode["FairnessConstraints"][1]["FA2"])

    # separation constraints
    constrSE = Vector{SEcon}(undef, 0)
    parse_se!(constrSE, constraintNode["SeparationConstraints"][1]["SE1"])

    rr2 = RR2(nteams, 2*nteams-2, isPhased)

    return Instance(rr2, constrCA, constrBR, constrGA, constrFA, constrSE)
end

"""
    parse_ca!(rr2::RR2, constrArray::Vector{CAcon}, node)

Parses all capacity constraints given by `node` and adds them to `constrArray`.
"""
function parse_ca!(rr2::RR2, constrArray::Vector{CAcon}, node)
    for constr in node
        attr = XML.attributes_dict(constr)
        ub = parse(Int64, attr["max"])
        pen = attr["type"] == "SOFT" ? parse(Int64, attr["penalty"]) : 0
        modestr = XML.name(constr) == "CA1" ? "mode" : "mode1"
        if attr[modestr] == "H"
            mode = MODE_H
        elseif attr[modestr] == "A"
            mode = MODE_A
        else
            mode = MODE_HA
        end

        if XML.name(constr) == "CA1"
            slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
            teams1 = IdxSet(2^parse(Int64, attr["teams"]))
            teams2 = IdxSet(2^rr2.nteams - 1)
            push!(constrArray, CAcon(teams1, teams2, slots, ub, mode, pen))
        elseif XML.name(constr) == "CA2"
            slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
            teams1 = IdxSet(parse.(Int64, split(attr["teams1"], ";")) .+ 1)
            teams2 = IdxSet(parse.(Int64, split(attr["teams2"], ";")) .+ 1)
            push!(constrArray, CAcon(teams1, teams2, slots, ub, mode, pen))
        elseif XML.name(constr) == "CA3"
            teams1 = parse.(Int64, split(attr["teams1"], ";"))
            teams2 = IdxSet(parse.(Int64, split(attr["teams2"], ";")) .+ 1)
            intp = parse(Int64, attr["intp"])
            binrepr = 2^intp - 1
            for team in teams1
                for start = 0:rr2.nslots-intp
                    push!(
                        constrArray,
                        CAcon(
                            IdxSet(2^team), # teams1
                            teams2,
                            IdxSet(binrepr << start), # slots
                            ub,
                            mode,
                            pen,
                        ),
                    )
                end
            end
        elseif XML.name(constr) == "CA4"
            teams1 = IdxSet(parse.(Int64, split(attr["teams1"], ";")) .+ 1)
            teams2 = IdxSet(parse.(Int64, split(attr["teams2"], ";")) .+ 1)
            if attr["mode2"] == "GLOBAL"
                slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
                push!(constrArray, CAcon(teams1, teams2, slots, ub, mode, pen))
            else
                for s in parse.(Int64, split(attr["slots"], ";"))
                    slots = IdxSet(2^s)
                    push!(
                        constrArray,
                        CAcon(teams1, teams2, slots, ub, mode, pen),
                    )
                end
            end
        else
            throw(ArgumentError("invalid node"))
        end
    end
end

"""
    parse_br!(rr2::RR2, constrArray::Vector{BRcon}, node)

Parses all break constraints given by `node` and adds them to `constrArray`.
"""
function parse_br!(constrArray, node)
    for constr in node
        attr = XML.attributes_dict(constr)
        slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
        pen = attr["type"] == "SOFT" ? parse(Int64, attr["penalty"]) : 0
        ub = parse(Int64, attr["intp"])
        modestr = XML.name(constr) == "BR1" ? "mode2" : "homeMode"
        if attr[modestr] == "H"
            mode = MODE_H
        elseif attr[modestr] == "A"
            mode = MODE_A
        else
            mode = MODE_HA
        end
        if XML.name(constr) == "BR1"
            for team in parse.(Int64, split(attr["teams"], ";"))
                teams = IdxSet(2^team)
                push!(constrArray, BRcon(teams, slots, ub, mode, pen))
            end
        else
            teams = IdxSet(parse.(Int64, split(attr["teams"], ";")) .+ 1)
            push!(constrArray, BRcon(teams, slots, ub, mode, pen))
        end
    end
end

"""
    parse_ga!(rr2::RR2, constrArray::Vector{GAcon}, node)

Parses all game constraints given by `node` and adds them to `constrArray`.
"""
function parse_ga!(constrArray, node)
    for constr in node
        attr = XML.attributes_dict(constr)
        pen = attr["type"] == "SOFT" ? parse(Int64, attr["penalty"]) : 0
        slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
        ub = parse(Int64, attr["max"])
        lb = parse(Int64, attr["min"])
        splittotuple = z -> tuple(parse.(Int64, split(z, ",")) .+ 1...)
        meetings = splittotuple.(split(strip(attr["meetings"], [';']), ";"))
        push!(constrArray, GAcon(slots, lb, ub, meetings, pen))
    end
end

"""
    parse_fa!(rr2::RR2, constrArray::Vector{FAcon}, node)

Parses all fairness constraints given by `node` and adds them to `constrArray`.
"""
function parse_fa!(constrArray, node)
    for constr in node
        attr = XML.attributes_dict(constr)
        teams = IdxSet(parse.(Int64, split(attr["teams"], ";")) .+ 1)
        slots = IdxSet(parse.(Int64, split(attr["slots"], ";")) .+ 1)
        ub = parse(Int64, attr["intp"])
        pen = attr["type"] == "SOFT" ? parse(Int64, attr["penalty"]) : 0
        push!(constrArray, FAcon(teams, slots, ub, pen))
    end
end

"""
    parse_se!(rr2::RR2, constrArray::Vector{SEcon}, node)

Parses all separation constraints given by `node` and adds them to `constrArray`.
"""
function parse_se!(constrArray, node)
    for constr in node
        attr = XML.attributes_dict(constr)
        teams = IdxSet(parse.(Int64, split(attr["teams"], ";")) .+ 1)
        pen = attr["type"] == "SOFT" ? parse(Int64, attr["penalty"]) : 0
        lb = parse(Int64, attr["min"])
        push!(constrArray, SEcon(teams, lb, pen))
    end
end

"""
    toxml(sched::Matrix{Int}, filename::String)

Writes a xml representation of the schedule `sched` to a file specified by
`filename`.
"""
function toxml(sched::Matrix{Int}, filename::String)
    open(filename, "w") do io
        write(
            io,
            "<Solution>
            <MetaData>
            <SolutionName>mySolution.xml</SolutionName>
            <InstanceName>myInstance.xml</InstanceName>
            <ObjectiveValue infeasibility=\"0\" objective=\"0\"/>
            </MetaData>
            <Games>",
        )
        for i = 1:size(sched, 1)
            for s = 1:size(sched, 2)
                if sign(sched[i, s]) == 1
                    write(
                        io,
                        "<ScheduledMatch home=\"$(i-1)\" away=\"$(sched[i,s]-1)\" slot=\"$(s-1)\"/>\n",
                    )
                end
            end
        end
        write(io, "</Games>\n</Solution>")
    end
end

function tofile(sched::Matrix{Int}, filename::String)
    open(filename, "w") do io
        for i = 1:size(sched, 1)
            write(io, join(sched[i,:], " "))
            write(io, "\n")
        end
    end
end


text:
using Kyulacs: Observable, QuantumCircuit, QuantumState
using Kyulacs.Gate: CNOT, Y, merge
using PyCall

using Test

@pyinclude("readme_example.py")

@testset "readme_example" begin
    state = QuantumState(3)
    seed = 0 # set random seed
    state.set_Haar_random_state(seed)

    circuit = QuantumCircuit(3)
    circuit.add_X_gate(0)
    merged_gate = merge(CNOT(0, 1), Y(1))
    circuit.add_gate(merged_gate)
    circuit.add_RX_gate(1, 0.5)
    circuit.update_quantum_state(state)

    observable = Observable(3)
    observable.add_operator(2.0, "X 2 Y 1 Z 0")
    observable.add_operator(-3.0, "Z 2")
    value = observable.get_expectation_value(state)
    @test value == py"readme_example"()
end


text:
@testset "IntervalsFold" begin
    i01 = Interval{Open,Closed}(0.0, 1.0)
    i23 = Interval{Open,Closed}(2.0, 3.0)
    i45 = Interval{Open,Closed}(4.0, 5.0)
    i67 = Interval{Open,Closed}(6.0, 7.0)
    i89 = Interval{Open,Closed}(8.0, 9.0)
    isf_dict = Dict([
        IntervalsFold(i01, 2.0, 5) => Dict(
            :pattern => i01,
            :gap => 2.0,
            :folds => 5,
            :length => 5,
            :size => 5.0,
            :unfold => [i01, i23, i45, i67, i89],
            :reverse => reverse([i01, i23, i45, i67, i89]),
        ),
    ])

    for (isf, results) in isf_dict
        @test pattern(isf) == results[:pattern]
        @test gap(isf) == results[:gap]
        @test folds(isf) == results[:folds]
        @test length(isf) == results[:length]
        @test size(isf) == results[:size]
        @test unfold(isf) == results[:unfold]
        @test ndims(isf) == 1
        for i in 1:1000
            @test rand(isf) ∈ isf
        end
        @test collect(isf) == [i for i in isf] == unfold(isf)
        @test collect(Iterators.reverse(isf)) == reverse(collect(isf)) == results[:reverse]
    end
    @test isempty(IntervalsFold(Interval{Open, Closed}(1.0, 1.0), 1.0, 1))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule boostOpenMPI_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("boostOpenMPI")
JLLWrappers.@generate_main_file("boostOpenMPI", UUID("a6d0cbf1-fcdd-5116-9b1a-1a7a207d6a38"))
end  # module boostOpenMPI_jll


text:
using Documenter
using ConcreteStructs

makedocs(
    modules = [ConcreteStructs],
    sitename = "ConcreteStructs.jl",
    pages =[
        "Home" => "index.md",
        "Walkthrough" => "walkthrough.md",
        "API" => "api.md",
    ],
    format = Documenter.HTML(
        canonical = "https://jonniedie.github.io/ConcreteStructs.jl/stable",
    ),
    repo="https://github.com/jonniedie/ConcreteStructs.jl/blob/{commit}{path}#L{line}",
    authors = "Jonnie Diegelman",
    assets = String[],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/jonniedie/ConcreteStructs.jl.git"
)