text:
x_dim = 28^2
z_dim = 2
hidden_dim = 128
batch_size = 100

train_dl = HDF5DataLayer(name="train-data",
        source=ENV["DATASETS"]"/mnist/train.txt",
        batch_size=batch_size)
test_dl = HDF5DataLayer(name="test-data",
        source=ENV["DATASETS"]"/mnist/test.txt",
        batch_size=batch_size)

enc_layers = [
    InnerProductLayer(name="enc1",
        output_dim=hidden_dim,
        neuron=Neurons.ReLU(),
        bottoms=[:data], tops=[:enc1]),
    SplitLayer(name="enc1_split",
        bottoms=[:enc1], tops=[:enc1a, :enc1b]),
    InnerProductLayer(name="z_mean",
        output_dim=z_dim,
        neuron=Neurons.Identity(),
        bottoms=[:enc1a], tops=[:z_mean]),
    InnerProductLayer(name="z_sd",
        output_dim=z_dim,
        bias_init=ConstantInitializer(-1.),
        weight_init=ConstantInitializer(0.),
        neuron=Neurons.Exponential(),
        bottoms=[:enc1b], tops=[:z_sd])]

sampling_layers = [
    SplitLayer(name="z_mean_split",
        bottoms=[:z_mean], tops=[:z_mean1, :z_mean2]),
    SplitLayer(name="z_sd_split",
        bottoms=[:z_sd], tops=[:z_sd1, :z_sd2]),

    RandomNormalLayer(output_dims=[z_dim],
        batch_sizes=[batch_size],
        tops=[:noise01]),
    ElementWiseLayer(operation=ElementWiseFunctors.Multiply(),
        bottoms=[:noise01, :z_sd1],
        tops=[:noise0sd]),
    ElementWiseLayer(name="z",
        operation=ElementWiseFunctors.Add(),
        bottoms=[:z_mean1, :noise0sd],
        tops=[:z]),

    GaussianKLLossLayer(name="encoder-loss",
        bottoms=[:z_mean2, :z_sd2])]

dec_layers = [
    InnerProductLayer(name="dec1", output_dim=hidden_dim,
        neuron=Neurons.Sigmoid(),
        bottoms=[:z], tops=[:dec1]),
    InnerProductLayer(name="x_mean", output_dim=x_dim,
        neuron=Neurons.Sigmoid(),
        bottoms=[:dec1], tops=[:x_mean])]

expected_recon_loss = BernoulliReconLossLayer(bottoms=[:x_mean, :data])

non_data_layers = [
    enc_layers...,
    sampling_layers...,
    dec_layers...,
    expected_recon_loss]


text:
# collection of GRAPE algorithms
using Zygote
using Optim

import Base.@kwdef

"""
Simple. Use Zygote to solve all of our problems
Compute the gradient of a functional f(x) wrt. x
"""
#  f_tol = 1e-3, g_tol = 1e-6, iters=1000
function _ADGRAPE(functional, x, optim_options)

    function grad_functional!(G, x)
        G .= real(Zygote.gradient(functional, x)[1])
    end

    res = optimize(functional, grad_functional!, x, LBFGS(), optim_options)

end

"""
Evaluate the figure of merit and the gradient (updated in-place) for a given specification of the problem type problem.
"""
function _fom_and_gradient_GRAPE!(
    A::T,
    B,
    control_array,
    n_timeslices,
    duration,
    n_controls,
    gradient,
    fwd_state_store,
    bwd_costate_store,
    propagators,
    X_init,
    X_target,
    evolve_store,
    system_type,
) where {T}

    dt = duration / n_timeslices

    fwd_state_store[1] .= X_init
    bwd_costate_store[end] .= X_target
    # this seems really stupid since we dont ever use the generators again, we can just keep the propagators instead
    # pw_ham_save!(A, B, control_array, n_controls, n_timeslices, @view generators[:])
    # @views propagators[:] .= exp.(generators[:] .* (-1.0im * dt))
    pw_prop_save!(A, B, control_array, n_controls, n_timeslices, dt, propagators)
    # function pw_prop_save!(H₀::T, Hₓ_array::Array{T,1}, x_arr::Array{Float64,2}, n_pulses, timeslices, timestep, out::Array{T, 1}) where T

    for t = 1:n_timeslices
        evolve_func!(
            system_type,
            t,
            fwd_state_store,
            bwd_costate_store,
            propagators,
            evolve_store,
            forward = true,
        )
    end

    for t in reverse(1:n_timeslices)
        evolve_func!(
            system_type,
            t,
            fwd_state_store,
            bwd_costate_store,
            propagators,
            evolve_store,
            forward = false,
        )
    end

    t = n_timeslices

    for c = 1:n_controls
        for t = 1:n_timeslices
            @views gradient[c, t] = grad_func!(
                system_type,
                t,
                dt,
                B[c],
                fwd_state_store,
                bwd_costate_store,
                propagators,
                evolve_store,
            )
        end
    end

    return fom_func(system_type, t, fwd_state_store, bwd_costate_store, propagators)

end

"""
Static GRAPE

Flexible GRAPE algorithm for use with StaticArrays where the size is always fixed. This works best if there are < 100 elements in the arrays. The result is that you can avoid allocations (this whole function allocates just 4 times in my tests). If your system is too large then try the GRAPE! algorithm above which should work for generic array types!
"""
function _fom_and_gradient_sGRAPE(
    A::T,
    B,
    control_array,
    n_timeslices,
    duration,
    n_controls,
    gradient,
    fwd_state_store,
    bwd_costate_store,
    X_init,
    X_target,
    system_type,
) where {T}

    dt = duration / n_timeslices
    # arrays that hold static arrays?
    fwd_state_store[1] = X_init
    bwd_costate_store[end] = X_target

    # now we want to compute the generators
    propagators = pw_evolve_save(A, B, control_array, n_controls, dt, n_timeslices)

    # forward evolution of states
    for t = 1:n_timeslices
        fwd_state_store[t+1] = evolve_func(
            system_type,
            t,
            fwd_state_store,
            bwd_costate_store,
            propagators,
            forward = true,
        )
    end
    # backward evolution of costates
    for t in reverse(1:n_timeslices)
        bwd_costate_store[t] = evolve_func(
            system_type,
            t,
            fwd_state_store,
            bwd_costate_store,
            propagators,
            forward = false,
        )
    end
    # update the gradient array
    for c = 1:n_controls
        for t = 1:n_timeslices
            gradient[c, t] = grad_func(
                system_type,
                t,
                dt,
                B[c],
                fwd_state_store,
                bwd_costate_store,
                propagators,
            )

        end
    end

    t = n_timeslices
    return fom_func(system_type, t, fwd_state_store, bwd_costate_store, propagators)
end

"""
Evolution functions for various GRAPE tasks
"""

"""
Function to compute the evolution for Unitary synthesis, since here we simply stack propagators
"""
function evolve_func(
    system_type::UnitaryGate,
    t,
    state,
    costate,
    propagator;
    forward = true,
)
    if forward
        propagator[t] * state[t]
    else
        propagator[t]' * costate[t+1]
    end
end

"""
Evolution function for use with density matrices
"""
function evolve_func(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    state,
    costate,
    propagator;
    forward = true,
)
    if forward
        propagator[t] * state[t] * propagator[t]'
    else
        propagator[t]' * costate[t+1] * propagator[t]
    end
end

"""
In-place evolution functions
"""
function evolve_func!(
    system_type::UnitaryGate,
    t,
    state,
    costate,
    propagator,
    store;
    forward = true,
)
    if forward
        @views mul!(state[t+1], propagator[t], state[t])
    else
        @views mul!(costate[t], propagator[t]', costate[t+1])
    end
end

"""
In-place evolution functions, I think these don't allocate at all
"""
function evolve_func!(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    state,
    costate,
    propagator,
    store;
    forward = true,
)
    if forward
        @views mul!(store, state[t], propagator[t]')
        @views mul!(state[t+1], propagator[t], store)
    else
        @views mul!(store, costate[t+1], propagator[t])
        @views mul!(costate[t], propagator[t]', store)
    end
end

"""
Grad functions for GRAPE
"""

"""
First order in dt gradient from Khaneja paper for unitary synthesis
"""
function grad_func!(
    system_type::UnitaryGate,
    t,
    dt,
    B,
    state,
    costate,
    props,
    store,
)::Float64
    @views mul!(store, state[t]', costate[t])
    @views 2.0 * real((1.0im * dt) .* tr(costate[t]' * B * state[t]) * tr(store))
end

function grad_func!(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    dt,
    B,
    state,
    costate,
    props,
    store,
)::Float64
    @views mul!(store, costate[t]', commutator(B, state[t]))
    real(tr((1.0im * dt) .* store))
end

function grad_func(system_type::UnitaryGate, t, dt, B, state, costate, props)
    2.0 * real((-1.0im * dt) * tr(costate[t]' * B * state[t]) * tr(state[t]' * costate[t]))
end

function grad_func(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    dt,
    B,
    state,
    costate,
    props,
)
    real(tr((1.0im * dt) .* (costate[t]' * commutator(B, state[t]))))
end

ext: jl


text:
using Flux
using LinearAlgebra: I # Identity matrix
using Distributions
using Random
using Plots

using Minijyro

Random.seed!(42)
theme(:ggplot2)

# Convenience functions
function get_num_params(nn)
    # For each parameter vector count the number of elements and collect the total
    # sum.
    return mapreduce(p -> prod(size(p)), +, Flux.params(nn))
end

function set_params!(nn, weights)
    new_params = Flux.params([])
    i = 1
    for p in Flux.params(nn)
        num_elements = prod(size(p))
        new_p = weights[i:i+num_elements-1]
        push!(new_params, reshape(new_p, size(p)))
        i += num_elements
    end
    Flux.loadparams!(nn, new_params)
end

@jyro function model(xs, neural_net, sigma)
    num_params = get_num_params(neural_net)
    w ~ MvNormal(num_params, 1.0)

    set_params!(neural_net, w)
    y ~ MvNormal(neural_net(xs), sigma*I)
end

D = 1 # Dimensionality of data
N = 100 # Number of data points
sigma = 0.1 # Noise level
dnn = Chain(Dense(D, 5, tanh), Dense(5, 1), x -> vec(x))

# Create some data
xs = collect(range(-3, stop=3, length=N))
xs = reshape(xs, (D,N))
ys = cos.(xs) + sigma * randn(D,N)
ys = reshape(ys, (N,))

display(scatter(xs[1,:], ys, grid=:none))

# Run NUTS sampler
cond_model = condition(model, Dict(:y => ys))
trace!(cond_model)
samples, stats = nuts(cond_model, (xs, dnn, sigma), 2000; autodiff=:reverse, progress=true)

# Plot posterior samples
p = plot(xlims=(-10,10))
for posterior_sample in samples[end-500:end]
    wide_xs = range(-10, stop=10, length=400)
    wide_xs = reshape(wide_xs, (1,400))
    set_params!(dnn, posterior_sample)
    preds = dnn(wide_xs)
    plot!(p, wide_xs[1,:], preds, color=:red, alpha=0.3)
end
display(scatter!(p, xs[1,:], ys, color=:blue, grid=:none, legend=false))


text:
using TrajectoryOptimization
using RobotDynamics
using BenchmarkTools
using ForwardDiff
using Test
using SparseArrays
using LinearAlgebra
const TO = TrajectoryOptimization

# Test NLPTraj iteration
n,m,N = 3,2,101
NN = RobotDynamics.num_vars(n,m,N)
@test NN == N*n + (N-1)*m
Z0 = Traj(n,m,0.1,N)
Zdata = TO.TrajData(Z0)
Z = rand(NN)
Z_ = TO.NLPTraj(Z,Zdata)
@test Z_[1] isa StaticKnotPoint
@test state(Z_[2]) == Z[(n+m) .+ (1:n)]
@test state(Z_[end]) == Z[end-n+1:end]
@test length(Z_) == N
@test size(Z_) == (N,)
@test [z for z in Z_] isa Vector{<:StaticKnotPoint}
@test eltype(Z_) == StaticKnotPoint{n,m,Float64,n+m}

# Test with problem
prob = DubinsCarProblem(:parallel_park)
prob.constraints
TO.add_dynamics_constraints!(prob)
n,m,N = size(prob)
cons = prob.constraints
NN = N*n + (N-1)*m
P = sum(TO.num_constraints(prob))

# Jacobian structure
jac = TO.JacobianStructure(cons)
TO.get_inds(cons[2],n,m)
TO.widths(cons[end])
@test jac.cinds[1][1] == 1:n
@test jac.cinds[2][1] == n .+ (1:4)
@test jac.cinds[5][1] == (n+4) .+ (1:n)
@test jac.cinds[5][end] == P-2n+1:P-n
@test jac.cinds[4][1] == P-n+1:P
@test jac.zinds[1][1] == 1:n
@test jac.zinds[2][1] == 1:n+m
@test jac.zinds[5][1] == 1:n+m
@test jac.zinds[5][1,2] == (n+m) .+ (1:n)
@test jac.zinds[5][end,1] == NN-2n-m+1:NN-n
@test jac.zinds[5][end,2] == NN-n+1:NN
@test jac.zinds[4][1] == NN-n+1:NN
@test jac.linds[1][1] == 1:n^2
@test jac.linds[2][1] == n^2 .+ (1:4*(n+m))
@test jac.linds[5][1] == (n^2 + 4*(n+m)) .+ (1:n*(n+m))

D = spzeros(P,NN)
TO.jacobian_structure!(D, jac)
@test D[1:n,1:n] == reshape(1:n^2,n,n)
@test Matrix(D[n .+ (1:4), 1:n+m]) == n^2 .+ reshape(1:4*(n+m), 4, n+m)
@test nnz(D) == D[end,end]

Dv = zeros(nnz(D))
d = zeros(P)
C,c = TO.gen_convals(Dv,d,cons)
@test C[1][1].parent.indices == (1:n^2,)

D = spzeros(P,NN)
C,c = TO.gen_convals(D,d,cons)
@test C[1][1].indices == (1:n,1:n)

# Test second-order constraint term
data = TO.NLPData(NN, P, jac.nD)
conSet = TO.NLPConstraintSet(prob.model, prob.constraints, data)
prob.Z[2].z += rand(n+m)
conSet.λ[5][1] .= rand(n)
TO.evaluate!(conSet, prob.Z)
@test conSet.convals[end].vals[1] != zeros(n)
@test conSet.convals[end].vals[3] == zeros(n)
TO.∇jacobian!(conSet.hess, conSet, prob.Z, conSet.λ)
@test conSet.hess[end][1] != zeros(n+m,n+m)
@test conSet.hess[end][2] == zeros(n+m,n+m)

# Build NLP
nlp = TrajOptNLP(prob)
Zdata = nlp.Z.Zdata
@test states(nlp) ≈ states(prob.Z)
@test controls(nlp) ≈ controls(prob.Z)

Z_ = TO.NLPTraj(prob.Z)
Z = Z_.Z
@test Z_[10].z == prob.Z[10].z
z = StaticKnotPoint(Z,Zdata,1)
@test state(z) ≈ state(prob.Z[1])

# Cost functions
@test cost(prob) ≈ TO.eval_f(nlp, Z)

@test TO.grad_f!(nlp, Z) ≈ ForwardDiff.gradient(x->TO.eval_f(nlp,x), Z)
# @btime TO.grad_f!($nlp, $Z)
# @btime ForwardDiff.gradient(x->TO.eval_f($nlp,x), $Z)

G = ForwardDiff.hessian(x->TO.eval_f(nlp,x), Z)
G0 = nlp.data.G
G0 .*= 0
@test TO.hess_f!(nlp, Z) ≈ G
@test nlp.E.cost[1].Q.parent === G0
# @btime TO.hess_f!($nlp, $Z)
# @btime ForwardDiff.hessian(x->TO.eval_f($nlp,x), $Z)

# Constraint Functions
initial_trajectory!(nlp, prob.Z)
# conSet = TO.ALConstraintSet(prob)
# TO.evaluate!(conSet, prob.Z)
# TO.evaluate!(nlp.conSet, prob.Z)
# c_max = max_violation(nlp.conSet)
# @test c_max ≈ max_violation(conSet)

c = TO.eval_c!(nlp, Z)
# @test max_violation(nlp) ≈ c_max
# @btime eval_c!($nlp, $Z)
# @btime evaluate!($nlp.conSet, $(prob.Z))

TO.jac_c!(nlp, Z)
D = nlp.data.D
@test D[1:n,1:n] == I(n)
@test D[end-n+1:end,end-n+1:end] == I(n)
@test D[end-2n+1:end-n,end-n+1:end] == -I(n)

# @btime jac_c!($nlp, $Z)

# Test Hessian lagrangian
nlp.conSet.λ[end][1] .= 0
@test TO.hess_L!(nlp, Z) ≈ G
nlp.conSet.λ[end][1] .= rand(n)
@test !(TO.hess_L!(nlp, Z) ≈ G)

# Test cost hessian structure
@test nlp.obj isa Objective{<:TO.DiagonalCostFunction}
G_ = TO.hess_f_structure(nlp)
@test nnz(G_) == NN
@test diag(G_) == 1:NN

obj_ = TO.QuadraticObjective(n,m,N)
prob_ = Problem(prob, obj=obj_)
nlp_ = TrajOptNLP(prob_)
@test !(nlp_.obj isa Objective{<:TO.DiagonalCostFunction})
G_ = TO.hess_f_structure(nlp_)
@test nnz(G_) == (N-1)*(n+m)^2 + n^2
@test G_[1:n+m, 1:n+m] == reshape(1:(n+m)^2, n+m, n+m)

r,c = TO.get_rc(G_)
@test [G_[r[i], c[i]] for i = 1:nnz(G_)] == 1:nnz(G_)

# Test jacobian structure
D_ = TO.jacobian_structure(nlp)
Matrix(D_)

@test D_[1:n,1:n] == reshape(1:n^2, n,n)
@test D_[n .+ (1:4), 1:n+m] == reshape(9 .+ (1:(n+m)*4), 4, n+m)
nlp.data.D .= 0
TO.jac_c!(nlp, Z)

nlp_ = TrajOptNLP(prob, jac_type=:vector)
D_ = TO.jacobian_structure(nlp_)
@test D_[1:n,1:n] == reshape(1:n^2, n,n)
@test D_[n .+ (1:4), 1:n+m] == reshape(9 .+ (1:(n+m)*4), 4, n+m)

nlp_ = TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
D_ = TO.jacobian_structure(nlp_)
@test D_[1:n,1:2n+m] == reshape(1:(2n+m)n, n, 2n+m)
@test Matrix(D_[n .+ (1:n), (n+m) .+ (1:2n+m)]) ≈ reshape((2n+m)n .+ (1:(2n+m)n), n, 2n+m)
@test nnz(D_) == (N-1)*(2n+m)*n

# Constraint type
IE = TO.constraint_type(nlp)
@test IE[1:n] == fill(:Equality,n)              # initial state constraint
@test IE[n .+ (1:4)] == fill(:Inequality,4)     # bound constraint
@test IE[(n+4) .+ (1:n)] == fill(:Equality,n)   # dynamics constraint

cL,cU = TO.constraint_bounds(nlp)
@test cL[1:n] == zeros(n)
@test cL[n .+ (1:4)] == fill(-Inf,4)     # bound constraint
@test cL[(n+4) .+ (1:n)] == zeros(n)  # dynamics constraint
@test all(cU .== 0 )

# Test view reset
nlp = TO.TrajOptNLP(prob)
conSet = TO.get_constraints(nlp)
TO.eval_c!(nlp)
TO.jac_c!(nlp)
@test conSet.convals[1].vals[1].parent === nlp.data.d
D = spzeros(P,NN)
d = zeros(P)
D0 = copy(nlp.data.D)
d0 = copy(nlp.data.d)
nlp.data.D = D
nlp.data.d = d

TO.reset_views!(conSet, nlp.data)
@test conSet.convals[1].vals[1].parent === d
@test conSet.convals[end].jac[1,2].parent === D
@test d == zeros(P)
@test D == spzeros(P,NN)
TO.eval_c!(nlp)
TO.jac_c!(nlp)
@test d == d0 != zeros(P)
@test D == D0 != spzeros(P,NN)

nlp_ = TO.TrajOptNLP(prob, jac_type=:vector)
conSet = TO.get_constraints(nlp_)
TO.eval_c!(nlp_)
TO.jac_c!(nlp_)
@test typeof(nlp_.conSet.convals[1].jac[1]) <: Base.ReshapedArray{<:Any,2,<:SubArray}
@test nlp_.data.d ≈ d0
@test conSet.convals[1].jac[1].parent.parent === nlp_.data.v
v = zero(nlp_.data.v)
d = zero(nlp_.data.d)
v0 = copy(nlp_.data.v)
d0 = copy(nlp_.data.d)
nlp_.data.v = v
nlp_.data.d = d

TO.reset_views!(conSet, nlp_.data)
@test conSet.convals[1].vals[1].parent === d
@test conSet.convals[end].jac[1,2].parent.parent === v
@test d ≈ zeros(P)
@test v == zero(v0)
TO.eval_c!(nlp_)
TO.jac_c!(nlp_)
@test d == d0 != zeros(P)
@test v == v0 != zero(v)

# Test view reset on cost expansion
nlp = TO.TrajOptNLP(prob)
TO.grad_f!(nlp)
TO.hess_f!(nlp)
G = spzeros(NN,NN)
g = zeros(NN)
G0 = copy(nlp.data.G)
g0 = copy(nlp.data.g)
nlp.data.G = G
nlp.data.g = g

TO.reset_views!(nlp.E, nlp.data)
@test nlp.E[1].Q.parent === G
@test nlp.E[1].R.parent === G
@test nlp.E[end].Q.parent === G
@test nlp.E[1].q.parent === g
@test nlp.E[N].q.parent === g
@test nlp.E[1].r.parent === g
TO.grad_f!(nlp)
TO.hess_f!(nlp)
@test nlp.data.G !== G0
@test nlp.data.G == G0
@test nlp.data.g !== g0
@test nlp.data.g == g0

# Test bounds removal
cons = prob.constraints
zL = fill(-Inf,NN)
zU = fill(+Inf,NN)
TO.primal_bounds!(zL, zU, cons, false)
@test zL[1:n] == zeros(n)
@test zU[1:n] == zeros(n)
@test zL[n+1:n+m] == fill(-2,m)
@test zU[n+m+1:2n+m] == [0.25, 1.501, Inf]
@test zL[n+m+1:2n+m] == [-0.25, -0.001, -Inf]
@test zL[end-n+1:end] == prob.xf
@test zU[end-n+1:end] == prob.xf

@test length(cons) == 5
cons2 = copy(cons)
TO.primal_bounds!(zL, zU, cons2, true)
@test length(cons2) == 1
@test cons2[1] isa TO.DynamicsConstraint
@test length(cons) == 5

prob = DubinsCarProblem(:parallel_park)
TO.add_dynamics_constraints!(prob)
n,m,N = size(prob)
cons = prob.constraints

nlp = TrajOptNLP(prob, remove_bounds=true)
@test length(prob.constraints) == 5  # make sure it doesn't modify the problem
@test length(nlp.conSet.convals) == 1
@test (zL,zU) == TO.primal_bounds!(nlp)
TO.jac_c!(nlp)

# Test vector sparsity
nlp_ = TrajOptNLP(prob, remove_bounds=true, jac_type=:vector)
TO.jac_c!(nlp_)
@test nlp_.conSet.convals[end].jac[1].parent.parent == nlp_.data.v
@test nlp_.data.v[1:n*(n+m)] == vec(nlp.data.D[1:n, 1:n+m])
@test all(TO.primal_bounds!(nlp_) .≈ (zL,zU))


text:
@testset "complete graph" begin

    testvideo = Video(400, 400)
    Background(1:30, ground())
    g = complete_graph(8)
    gd = JGraphData(g, NetworkLayout.Shell(),  width = 100, height = 100)
    jg = JGraph(gd)
    for node in jnodes(jg)
        act!(node, Action(anim_rotate_around(2π, O)))
    end
    render(testvideo, tempdirectory = "images", pathname = "")

    @test_reference "refs/test_complete_graph_5.png" load("images/$(lpad(5, 10, "0")).png")

    for image in readdir("images", join = true)
        endswith(image, "png") && rm(image)
    end
end

@testset "morph" begin
    testvideo = Video(400, 400)
    Background(1:30, ground())
    g = complete_graph(8)
    gd = JGraphData(g, NetworkLayout.Shell(), width = 100, height = 100)
    jg = JGraph(gd)

    jgraph_morph(jg, [NetworkLayout.Spring(), NetworkLayout.Shell()], [(50, 50), (100, 100)])
    render(testvideo, tempdirectory = "images", pathname = "")

    for frame_id in [7, 15, 23, 30]
        @test_reference "refs/test_morph_$(frame_id).png" load(
            "images/$(lpad(frame_id, 10, "0")).png",
        )
    end

    for image in readdir("images", join = true)
        endswith(image, "png") && rm(image)
    end
end

@testset "jgraph_walk" begin
    testvideo = Video(400, 400)
    Background(1:30, ground())
    g = complete_graph(8)
    gd = JGraphData(g, NetworkLayout.Shell(), width = 100, height = 100)
    jg = JGraph(gd)

    jgraph_walk(jg, g -> [vertices(g); 1])
    render(testvideo, tempdirectory = "images", pathname = "")

    for frame_id in [1, 8, 15, 23, 30]
        @test_reference "refs/test_jgraph_walk_$(frame_id).png" load(
            "images/$(lpad(frame_id, 10, "0")).png",
        )
    end

    for image in readdir("images", join = true)
        endswith(image, "png") && rm(image)
    end
end

@testset "jgraph_walk jgraph_morph" begin
    testvideo = Video(400, 400)
    Background(1:100, ground("white", "black"))
    g = barbell_graph(5, 10)
    gd = JGraphData(g, NetworkLayout.Spring(), width = 40, height = 40)
    jg = JGraph(gd)
    jgraph_walk(jg, g -> vertices(g))
    jgraph_morph(jg, NetworkLayout.Shell(), (100, 100))
    render(testvideo, tempdirectory = "images", pathname = "")

    for frame_id in [1, 5, 10, 50, 75, 100]
        @test_reference "refs/test_walk_morph_$(frame_id).png" load(
            "images/$(lpad(frame_id, 10, "0")).png",
        )
    end

    for image in readdir("images", join = true)
        endswith(image, "png") && rm(image)
    end
end


text:
mutable struct StandardSimplexe{T} <: AbstractStandard{T}
    M::Array{T, 2}
    A::Array{T, 2}
    b::Array{T, 1}
    c::Array{T, 1}
    Binv::Array{T, 2}
    b_idx_initial::Array{Int, 1}
    xstar::Array{T, 1}
    vstar::T
    b_idx::Array{Int, 1}
    status::AbstractStatus
    function StandardSimplexe(A::Array{T, 2}, b::Array{T, 1}, c::Array{T, 1}; verbose::Bool = true, b_idx::Array{Int, 1} = [-1]) where T
        m,n = size(A)
        @assert !(typeof(T) <: Integer)  "Type $T cannot be a subtype of Integer"
        @assert length(b) == m "dimension of A and b mismatch, size(A) = ($m, $n), length(b) = $(length(b)) != $m"
        @assert length(c) == n "dimension of A and c mismatch, size(A) = ($m, $n), length(c) = $(length(c)) != $n"
        @assert rank(A) == m "A is not a full rank Matrix"
        ss = new{T}()
        ss.xstar = T[],
        ss.vstar = T(Inf)
        ss.status = status
        ss.A = copy(A)
        ss.b = copy(b)
        ss.c = copy(c)

        if b_idx != [-1]
            #if a first base is given, we use it
            M = [A b; c' T(0)]
            @assert all(areisolate(M[:, b_idx[i]], i) for i in 1:m) "base given contains non isolated index"
            BInv = Array{T, 2}(I, m, m)
            ss.Binv = Binv
            ss.b_idx = copy(b_idx)
            ss.b_idx_initial = copy(b_idx)
            ss.M = M
        else
            Binv = Array{T, 2}(I, m, m)
            for k in 1:m
                if b[k] <  0
                    Binv[i, i] = -one(T)
                    A[k, :] = -A[k, :]
                    b[k] = -b[k]
                end
            end
            M, b_idx = phase1(A,b,c, verbose = verbose)
            B = ss.A[:, b_idx]
            Binv = Binv*inv(B)#could be much better. Lazy code
            ss.Binv = Binv
            ss.b_idx = copy(b_idx)
            ss.b_idx_initial =  copy(b_idx)
            ss.M = M
            status = (b_idx == [-1]) ? Infeasible() : Unknown()
        end
        return ss
    end
end
function areisolate(v::Vector{T}, i::Int) where T
    m = length(v)
    return v[i] == one(T) && all(v[setdiff(1:m, i)] .== zero(T))
end

function xstar(ss::StandardSimplexe{T}) where T
    ss.status != Optimal() && @warn "probleme not optimal"
    return ss.xstar
end
function vstar(ss::StandardSimplexe{T}) where T
    ss.status != Optimal() && @warn "probleme not optimal"
    return ss.vstar
end

function println(ss::StandardSimplexe)
    m, n = size(ss.M)
    prd = "+" * prod("-------+" for _ in 1:n+1)
    println(prd)
    print("| V.B.  |")
    for k in 1:n-1
        s = ""
        xk = " x_$k "
        s*= xk
        s*= prod(" " for _ in 1:7-length(xk))
        print(s)
        print("|")
    end

    print(" T.D.  |")
    println()
    println(prd)
    for i in 1:m
        if i != m
            s = "| "#------"
            xk = "x_$(ss.b_idx[i])"
            s*= xk
            s*= prod(" " for _ in 1:6-length(xk))
            print(s)
        else
            print("| -z    ")
        end
        for j in 1:n
            print("| ")
            str = string(Float64(ss.M[i, j]))
            if length(str) >= 5
                str = str[1:5]
            end
            print(str)
            if length(1:(5 - length(str))) != 0
                print(prod(" "  for _ in 1:(5 - length(str))))
            end
            print(" ")
        end
        println("|")
        println(prd)
    end
end

function phase1(A::Matrix{T},b::Vector,c::Vector; verbose::Bool = false) where  T
    m, n = size(A)
    base = -ones(Int, m)
    for k in 1:n
        if count(t -> t > 0, A[:, k]) == 1#recherche d'une variable isolée
            i = findfirst(t -> t > 0,  A[:, k])
            base[i] = k
        end
    end
    #we need to add an artificial variable for all line with no isolated variables
    added = findall(t -> t == -1, base)

    for (i, j) in enumerate(base)
        if j > 0
            pivot!(A, i, j)
        end
    end

    if length(added) == 0
        M = [A b; c' T(0)]
        for (i,j) in enumerate(base)
            pivot!(M, i, j)
        end
        return M, base
    end

    M = [A hcat([ei(T, k, m) for k in added]...) b; zeros(T, 1, n) ones(T, 1, length(added)) zero(T)]

    base[added] = n+1:n+length(added)
    for k in 1:length(base)
        if base[k] > n
            pivot!(M, k, base[k])
        end
    end

    phase1 = StandardSimplexe(M, base)
    verbose && println("simplexe array of phase 1")
    verbose && println(phase1)
    solve!(phase1)
    verbose && println("solved simplexe array of phase 1")
    verbose && println(phase1)
    #if the simplexe array  is optimal but the objective value is different then 0
    if (phase1.status == Optimal()) && (phase1.vstar != zero(T))
        println("probleme is unfeasible")
        return [phase1.M[1:m, 1:n] phase1.M[1:m, end]; c' zero(T)], [-1]
    else
        @assert issubset(base, 1:n) "artificial variable in base of phase 1"
        M = [phase1.M[1:m, 1:n] phase1.M[1:m, end]; c' zero(T)]
        base = phase1.b_idx
        for (i,j) in enumerate(base)
            pivot!(M, i, j)
        end

        return M, base
    end
end

function isOptimal(ss::StandardSimplexe{T}) where T
    return all(ss.M[end, 1:end-1] .>= -10*eps())
end

function solve!(ss::StandardSimplexe{T}; verbose::Bool = false, kmax = 1000, pivotrule::Function = findpivot) where T
    k = 0
    if ss.status == Infeasible()
        return ss.status
    end
    m2, n2 = size(ss.M)
    m, n = m2 - 1,  n2 -1
    while !isOptimal(ss) && (k+=1) <= kmax
        verbose && println(ss)
        i, j = pivotrule(ss, verbose = verbose)
        if i == -1
            ss.status = Unbounded()
            break
        end
        pivot!(ss.M, i, j)
        ss.b_idx[i] = j
    end
    if isOptimal(ss)
        ss.status = Optimal()
        ss.xstar = zeros(T, n)
        ss.xstar[ss.b_idx] = ss.M[1:end-1, end]
        ss.vstar = -ss.M[end, end]
    end
    status = ss.status
    verbose &&(println("$status Simplexe");  println(ss))
    ss.Binv = ss.M[1:end-1, b_idx_initial] * ss.Binv
    ss.status
end

function findpivot(ss::StandardSimplexe{T}; verbose::Bool = false) where T
    M = ss.M
    m = size(M, 2) - 1 #changed 1 to 2
    c = ss.M[end, 1:end-1]
    b=  ss.M[1:end-1, end]
    entering = argmin(b) #argmin(c) changed to b
    verbose && @show entering
    @assert c[entering] < 0 "simplexe is in an optimal state" #ss was optimal
    val = T(Inf)
    leaving = -1
    for k in 1:m
        if (M[k, entering] > 0) && ((tobeat = c[k]/M[k, entering]) < val) #changed b to c
            leaving = k
            val = tobeat
        end
    end
    verbose && @show leaving
    return leaving, entering
end

function findpivotUsingBlandrules(ss::StandardSimplexe{T}; verbose::Bool = false) where T
    M = ss.M
    m = size(M, 1) - 1
    c = ss.M[end, 1:end-1]
    b=  ss.M[1:end-1, end]

    entering = findfirst(t -> t < 0, c)
    verbose && @show entering
    @assert !isnothing(entering) "simplexe is in an optimal state"
    yls = M[1:end-1, entering]
    val = T(Inf)
    leaving = -1
    for k in 1:m
        if (M[k, entering] > 0) && ((tobeat = b[k]/M[k, entering]) < val)
            leaving = k
            val = tobeat
        end
    end
    verbose && @show leaving
    return leaving, entering
end


text:
module HierarchicalTemporalMemory

using SparseArrays, Parameters
using Lazy: @>, @>>
import Setfield: @set
import ImageFiltering: mapwindow, imfilter!, Fill
import Random: rand!, randsubseq!, randperm
import LinearAlgebra: Adjoint
import Statistics: mean, median, quantile
import StatsBase: countmap
import Chain: @chain
import Graphs: SimpleDiGraph, add_edge!, edges, vertices, inneighbors, outneighbors, a_star

# Type aliases

const Option{T}= Union{T,Nothing}
const Maybe{T}= Union{T,Missing}
const VecInt{T<:Integer}= Union{Vector{T}, T}
"𝕊𝕢 is the type of connection permanences, defining their quantization domain"
const 𝕊𝕢= UInt8
"𝕊𝕢range is the domain of connection permanence quantization"
const 𝕊𝕢range= 𝕊𝕢(0):typemax(𝕊𝕢)

"Type of neuron layer activations"
const CellActivity= BitArray
const DenseSynapses= Matrix{𝕊𝕢}
const SparseSynapses= SparseMatrixCSC{𝕊𝕢}
const AnySynapses= Union{DenseSynapses,SparseSynapses}
const DenseConnection= Matrix{Bool}
const SparseConnection= SparseMatrixCSC{Bool}
const AnyConnection= Union{DenseConnection,SparseConnection}

include("utils/tuple_utils.jl")
include("utils/iter_utils.jl")
include("utils/array_utils.jl")
include("utils/arithmetic_utils.jl")
include("utils/topology.jl")

# Algorithm imports

include("algorithm_parameters.jl")
include("dynamical_systems.jl")
include("SpatialPooler.jl")
include("TemporalMemory.jl")
include("Region.jl")
include("Network.jl")

export Hypercube, Hypersphere
export SpatialPooler, SPParams
export TemporalMemory, TMParams
export Region, Network
export step!, reset!, Nc, Nₙ, proximalSynapses, distalSynapses
export networkPlot, propagation_delay

# Maybe move encoders to another package, as they're many and independent from the core algorithms?
# TODO create extensible encoder interface
# encode(:arithmetic,...), encode(:powerday,...) etc
include("encoder.jl")
export encode_simpleArithmetic, initenc_simpleArithmetic,
  encode_powerDay, initenc_powerDay

# TODO same for decoders
include("decoder.jl")
export SDRClassifier, predict!, reverse_simpleArithmetic

end # module


text:
@defcomp grosseconomy begin
    regions = Index()

    YGROSS = Variable(index=[time, regions]) # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    K = Variable(index=[time, regions]) # Capital stock (trillions 2005 US dollars)

    al = Parameter(index=[time, regions]) # Level of total factor productivity
    l = Parameter(index=[time, regions]) # Level of population and labor
    I = Parameter(index=[time, regions]) # Investment (trillions 2005 USD per year)
    gama = Parameter() # Capital elasticity in production function
    dk = Parameter(index=[regions]) # Depreciation rate on capital (per year)
    k0 = Parameter(index=[regions]) # Initial capital value (trill 2005 USD)

    # TODO remove this, just a temporary output trick
    L = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        #Define function for K
        for r in d.regions
            if is_first(t)
                v.K[t,r] = p.k0[r]
            else
                v.K[t,r] = (1 - p.dk[r])^10 * v.K[t-1,r] + 10 * p.I[t-1,r]
            end
        end

        #Define function for YGROSS
        for r in d.regions
            v.YGROSS[t,r] = (p.al[t,r] * (p.l[t,r]/1000)^(1-p.gama)) * (v.K[t,r]^p.gama)
        end

        # TODO remove this, just a temporary output trick
        for r in d.regions
            v.L[t,r] = p.l[t,r]
        end
    end
end

max_stars_repo_name: fdennig/RICEupdate

max_issues_repo_name: fdennig/RICEupdate

max_forks_repo_name: fdennig/RICEupdate


text:
# Implementation of basis Galois field operations
# based on the public domain C code from
# http://www.samiam.org/galois.html
# by Sam Trenholme

# Log table using 0xe5 (229) as the generator
global const LTABLE = [
0x00, 0xff, 0xc8, 0x08, 0x91, 0x10, 0xd0, 0x36,
0x5a, 0x3e, 0xd8, 0x43, 0x99, 0x77, 0xfe, 0x18,
0x23, 0x20, 0x07, 0x70, 0xa1, 0x6c, 0x0c, 0x7f,
0x62, 0x8b, 0x40, 0x46, 0xc7, 0x4b, 0xe0, 0x0e,
0xeb, 0x16, 0xe8, 0xad, 0xcf, 0xcd, 0x39, 0x53,
0x6a, 0x27, 0x35, 0x93, 0xd4, 0x4e, 0x48, 0xc3,
0x2b, 0x79, 0x54, 0x28, 0x09, 0x78, 0x0f, 0x21,
0x90, 0x87, 0x14, 0x2a, 0xa9, 0x9c, 0xd6, 0x74,
0xb4, 0x7c, 0xde, 0xed, 0xb1, 0x86, 0x76, 0xa4,
0x98, 0xe2, 0x96, 0x8f, 0x02, 0x32, 0x1c, 0xc1,
0x33, 0xee, 0xef, 0x81, 0xfd, 0x30, 0x5c, 0x13,
0x9d, 0x29, 0x17, 0xc4, 0x11, 0x44, 0x8c, 0x80,
0xf3, 0x73, 0x42, 0x1e, 0x1d, 0xb5, 0xf0, 0x12,
0xd1, 0x5b, 0x41, 0xa2, 0xd7, 0x2c, 0xe9, 0xd5,
0x59, 0xcb, 0x50, 0xa8, 0xdc, 0xfc, 0xf2, 0x56,
0x72, 0xa6, 0x65, 0x2f, 0x9f, 0x9b, 0x3d, 0xba,
0x7d, 0xc2, 0x45, 0x82, 0xa7, 0x57, 0xb6, 0xa3,
0x7a, 0x75, 0x4f, 0xae, 0x3f, 0x37, 0x6d, 0x47,
0x61, 0xbe, 0xab, 0xd3, 0x5f, 0xb0, 0x58, 0xaf,
0xca, 0x5e, 0xfa, 0x85, 0xe4, 0x4d, 0x8a, 0x05,
0xfb, 0x60, 0xb7, 0x7b, 0xb8, 0x26, 0x4a, 0x67,
0xc6, 0x1a, 0xf8, 0x69, 0x25, 0xb3, 0xdb, 0xbd,
0x66, 0xdd, 0xf1, 0xd2, 0xdf, 0x03, 0x8d, 0x34,
0xd9, 0x92, 0x0d, 0x63, 0x55, 0xaa, 0x49, 0xec,
0xbc, 0x95, 0x3c, 0x84, 0x0b, 0xf5, 0xe6, 0xe7,
0xe5, 0xac, 0x7e, 0x6e, 0xb9, 0xf9, 0xda, 0x8e,
0x9a, 0xc9, 0x24, 0xe1, 0x0a, 0x15, 0x6b, 0x3a,
0xa0, 0x51, 0xf4, 0xea, 0xb2, 0x97, 0x9e, 0x5d,
0x22, 0x88, 0x94, 0xce, 0x19, 0x01, 0x71, 0x4c,
0xa5, 0xe3, 0xc5, 0x31, 0xbb, 0xcc, 0x1f, 0x2d,
0x3b, 0x52, 0x6f, 0xf6, 0x2e, 0x89, 0xf7, 0xc0,
0x68, 0x1b, 0x64, 0x04, 0x06, 0xbf, 0x83, 0x38
]

# Anti-log table
global const ATABLE = [
0x01, 0xe5, 0x4c, 0xb5, 0xfb, 0x9f, 0xfc, 0x12,
0x03, 0x34, 0xd4, 0xc4, 0x16, 0xba, 0x1f, 0x36,
0x05, 0x5c, 0x67, 0x57, 0x3a, 0xd5, 0x21, 0x5a,
0x0f, 0xe4, 0xa9, 0xf9, 0x4e, 0x64, 0x63, 0xee,
0x11, 0x37, 0xe0, 0x10, 0xd2, 0xac, 0xa5, 0x29,
0x33, 0x59, 0x3b, 0x30, 0x6d, 0xef, 0xf4, 0x7b,
0x55, 0xeb, 0x4d, 0x50, 0xb7, 0x2a, 0x07, 0x8d,
0xff, 0x26, 0xd7, 0xf0, 0xc2, 0x7e, 0x09, 0x8c,
0x1a, 0x6a, 0x62, 0x0b, 0x5d, 0x82, 0x1b, 0x8f,
0x2e, 0xbe, 0xa6, 0x1d, 0xe7, 0x9d, 0x2d, 0x8a,
0x72, 0xd9, 0xf1, 0x27, 0x32, 0xbc, 0x77, 0x85,
0x96, 0x70, 0x08, 0x69, 0x56, 0xdf, 0x99, 0x94,
0xa1, 0x90, 0x18, 0xbb, 0xfa, 0x7a, 0xb0, 0xa7,
0xf8, 0xab, 0x28, 0xd6, 0x15, 0x8e, 0xcb, 0xf2,
0x13, 0xe6, 0x78, 0x61, 0x3f, 0x89, 0x46, 0x0d,
0x35, 0x31, 0x88, 0xa3, 0x41, 0x80, 0xca, 0x17,
0x5f, 0x53, 0x83, 0xfe, 0xc3, 0x9b, 0x45, 0x39,
0xe1, 0xf5, 0x9e, 0x19, 0x5e, 0xb6, 0xcf, 0x4b,
0x38, 0x04, 0xb9, 0x2b, 0xe2, 0xc1, 0x4a, 0xdd,
0x48, 0x0c, 0xd0, 0x7d, 0x3d, 0x58, 0xde, 0x7c,
0xd8, 0x14, 0x6b, 0x87, 0x47, 0xe8, 0x79, 0x84,
0x73, 0x3c, 0xbd, 0x92, 0xc9, 0x23, 0x8b, 0x97,
0x95, 0x44, 0xdc, 0xad, 0x40, 0x65, 0x86, 0xa2,
0xa4, 0xcc, 0x7f, 0xec, 0xc0, 0xaf, 0x91, 0xfd,
0xf7, 0x4f, 0x81, 0x2f, 0x5b, 0xea, 0xa8, 0x1c,
0x02, 0xd1, 0x98, 0x71, 0xed, 0x25, 0xe3, 0x24,
0x06, 0x68, 0xb3, 0x93, 0x2c, 0x6f, 0x3e, 0x6c,
0x0a, 0xb8, 0xce, 0xae, 0x74, 0xb1, 0x42, 0xb4,
0x1e, 0xd3, 0x49, 0xe9, 0x9c, 0xc8, 0xc6, 0xc7,
0x22, 0x6e, 0xdb, 0x20, 0xbf, 0x43, 0x51, 0x52,
0x66, 0xb2, 0x76, 0x60, 0xda, 0xc5, 0xf3, 0xf6,
0xaa, 0xcd, 0x9a, 0xa0, 0x75, 0x54, 0x0e, 0x01
]

# Addition over GF(2^8)
function gadd(a::UInt8, b::UInt8)
  xor.(a , b)
end

function gadd(v::Array{UInt8})
  s = v[1]
  for x in v[2:end]
    s = gadd(s, x)
  end
  return s
end

# Subtraction over GF(2^8)
function gsub(a::UInt8, b::UInt8)
  xor.(a , b)
end

# Fast multiplication using (anti-)logarithm table.
function gmul(a::UInt8, b::UInt8)
	loga = Int(LTABLE[Int(a) + 1])
	logb = Int(LTABLE[Int(b) + 1])
	s = ATABLE[mod(loga + logb, 255) + 1]
	# Attempt to resist timing attacks.
	# Return 0x00 if a is zero or if b is zero.
	# Otherwise, return s.
	z = 0x00
	q = s
	if a == 0x00
		s = z
	else
		s = q
	end
	if b == 0x00
		s = z
	else
		q = s
	end
	return s
end

# Multiplication over GF(2^8) defined by the polynomial
# x^8 + x^4 + x^3 + x + 1 = 0.
function gmulbits(a::UInt8, b::UInt8)
	p = 0x00

	for counter=1:8
		if (b & 0x01) == 0x01
			p $= a
		end
		# detect if x^8 term is about to be generated
		carry = a & 0x80
		a <<= 1
		if carry == 0x80
			# replace x^8 with x^4 + x^3 + x + 1
			a $= 0x1B
		end
		b >>= 1
	end

	return p
end

# Returns the multiplicative inverse of a.
function gmulinv(a::UInt8)
	if a == 0x00
		# as defined by AES
		return 0x00
	else
		return ATABLE[256 - LTABLE[Int(a) + 1]]
	end
end


text:
## Defines functions to be used for active learning

# Fetch upper triangle of adj matrix and create parameter vector
# Ref: https://gist.github.com/tpapp/87e5675b87dbdd7a9a840e546eb20fae
function vec_triu_loop(adj_matrix::Array)
    m, n = size(adj_matrix)
    m == n || throw(error("not square"))
    len_theta = n*(n+1) ÷ 2

    # Array to hold values of adjacency matrix
    theta = ones(len_theta)

    # Create array of tuples to hold edges of corresponding values in theta
    edges_theta = Array{Tuple{Int, Int}}(undef, len_theta)

    k = 0
    @inbounds for i in 1:n
        for j in 1:i
            theta[k + j] = adj_matrix[j, i]
            edges_theta[k+j] = (j,i)
        end
        k += i
    end
    return len_theta, theta, edges_theta
end

function remove_node_from_edge(edge::Tuple, node::Integer)
    #=
    Functionality (1,2)\{2} = 1. Note that (1,1)\{1} = 1 as desired here.
    =#
    if edge[1] == node
        return edge[2]
    elseif edge[2] == node
        return edge[1]
    else
        return nothing
    end
end

# Function to calculate the fisher information matrix for a given query
function fisher_information_query(adj_matrix::Array,theta::Array,edges_theta::Array,sigma::Array,len_theta::Integer,num_spins::Integer)
    #=
    Inputs:
    adjacency matrix of the graphical model (estimated) -> corresponds to your theta
    sigma which is the particular query
    =#

    # Define Fisher Information matrix
    I_x = zeros((len_theta, len_theta))

    # Start filling up the upper triangle of I_x
    for a in 1:len_theta
        for b in a:len_theta
            # Get common node if exists
            set_common_nodes_temp = intersect(edges_theta[a],edges_theta[b])

            # if s is empty do nothing (i.e., I_x[a,b] = 0) else
            if !(isempty(set_common_nodes_temp))
                for s in set_common_nodes_temp
                    ind1 = remove_node_from_edge(edges_theta[a],s)
                    ind2 = remove_node_from_edge(edges_theta[b],s)

                    # set shared spin s to 1 as A_s = \sum_{u in ∂s} J_{su} σ_u + J_{ss}
                    sigma_temp = copy(sigma)
                    sigma_temp[s] = 1
                    A_s = sum(sigma_temp[i]*adj_matrix[s,i] for i=1:num_spins)

                    # We pick from sigma_temp instead of sigma for scenarios such as (1,1)∩(1,2)
                    I_x[a,b] += (1.0/num_spins)*sigma_temp[ind1]*sigma_temp[ind2]*(1.0 - tanh(A_s)^2)
                end
            end
        end
    end

    # Now copy over the upper triangle to lower triangle
    D = Diagonal(I_x)
    I_x_sut = I_x - D
    I_x = I_x_sut + transpose(I_x_sut) + D
    return I_x
end

# Function to calculate the fisher information matrix for a given query distribution
function fisher_information_query_distrn(adj_matrix::Array,theta::Array,edges_theta::Array,X_σ::Array,q_σ::Array,len_theta::Integer,num_spins::Integer)
    #=
    Inputs:
    X_σ is the set of queries we are considering
    q_σ is the query distribution over these queries
    =#

    # Make sure the given query distribution is a valid pdf
    #if sum(q_σ)!=1
    #    error("Given query distribution is not a valid pdf")
    #end

    I_q = zeros((len_theta, len_theta))
    for ind_σ=1:length(X_σ)
        sigma0 = 2*digits(X_σ[ind_σ], base=2, pad=num_spins).-1
        I_q += q_σ[ind_σ]*fisher_information_query(adj_matrix,theta,edges_theta,sigma0,len_theta,num_spins)
    end

    return I_q
end

# Function to calculate the fisher information matrix for a given query distribution where q is a Variable (Convex.jl)
function fisher_information_query_distrn_obj(adj_matrix::Array,theta::Array,edges_theta::Array,X_σ::Array,q_σ::Variable,len_theta::Integer,num_spins::Integer)
    #=
    Inputs:
    X_σ is the set of queries we are considering
    q_σ is the query distribution over these queries
    =#

    I_q = zeros((len_theta, len_theta))
    for ind_σ=1:length(X_σ)
        sigma0 = 2*digits(X_σ[ind_σ], base=2, pad=num_spins).-1
        I_q += q_σ[ind_σ]*fisher_information_query(adj_matrix,theta,edges_theta,sigma0,len_theta,num_spins)
    end

    return I_q
end

# Function to calculate the fisher information matrix for a given query distribution
function active_learning_M(adj_matrix::Array,X_σ::Array,num_spins::Integer)
    #=
    Inputs:
    adj_matrix is the current estimate of the parameters
    X_σ is the set of queries we are considering
    q_σ is the query distribution over these queries

    Ref: https://www.juliaopt.org/Convex.jl/dev/types/
    =#

    # Create the parameter vector: θ = [J_11, J_12, J_13, ..., J1n, ..., J_{i,i+1}, ..., J_nn]
    len_theta, theta, edges_theta = vec_triu_loop(adj_matrix)

    # Create the variables required for the SDP
    len_X = length(X_σ)
    q = Variable(len_X)  # query distribution
    t = Variable(len_theta)  # auxillary variables

    # Constraints for the query distribution
    constraints = [q>=0, q<=1, sum(q)==1]

    # Constraints involving the matrices
    U = 1.0*Matrix(I,len_theta,len_theta)
    for i=1:len_theta
        constraints += [fisher_information_query_distrn_obj(adj_matrix, theta, edges_theta, X_σ, q, len_theta, num_spins) U[i,:]; transpose(U[i,:]) t[i]] ⪰ 0
        #constraints += [t[i] transpose(U[i,:]); U[i,:] fisher_information_query_distrn_obj(adj_matrix, theta, edges_theta, X_σ, q, len_theta, num_spins)] ⪰ 0
    end

    # Solve the SDP
    problem = minimize(sum(t), constraints)
    #solve!(problem, SCS.Optimizer)
    solve!(problem, Mosek.Optimizer)

    return q.value
end

# Function to calculate the entropy of a given query
function entropy_query(adj_matrix::Array,sigma::Array,num_spins::Integer)
    #=
    Inputs:
    adjacency matrix of the graphical model (estimated) -> corresponds to your theta
    sigma which is the particular query
    =#

    # Define entropy
    entropy = 0.0

    # Start filling up the upper triangle of I_x
    for k in 1:num_spins
        sigma_temp = copy(sigma)
        sigma_temp[k] = 1
        A_k = sum(sigma_temp[i]*adj_matrix[k,i] for i=1:num_spins)

        entropy += log(2*cosh(A_k)) - A_k*tanh(A_k)
    end

    return entropy
end

# Function to calculate the entropy for a given set of queries/configurations
function entropy_configs(adj_matrix::Array,X_σ::Array,num_spins::Integer)
    #=
    Inputs:
    adjacency matrix of the graphical model (estimated) -> corresponds to your theta
    X_σ which is the set of queries
    =#

    # Define entropy
    entropy_configs = zeros(length(X_σ))

    for i in 1:length(X_σ)
        sigma_temp = 2*digits(X_σ[i], base=2, pad=num_spins) .- 1
        entropy_configs[i] = (1/length(X_σ))*entropy_query(adj_matrix, sigma_temp, num_spins)
    end

    return entropy_configs
end

# Unit test (by visualization) to check if sampling of M-regime with query was done correctly
function unit_test_sampling_M_regime_query(samples_U::Array, spin_number::Integer)
    config_number = 2^spin_number
    spin_configs = []
    for i = 0:(config_number - 1)
        spin_tmp = 2*digits(i, base=2, pad=spin_number).-1
        push!(spin_configs, string(spin_tmp))
    end
    dict_spin_configs = Dict(spin_configs[i+1] => i for i = 0:(config_number-1))

    bins_U_sigma0 = Array{Int}(undef, size(samples_U)[1],2)
    for i=1:size(samples_U)[1]
        sigma0_temp = samples_U[i,3:(n+2)]
        config0_temp = dict_spin_configs[string(sigma0_temp)]
        bins_U_sigma0[i,2] = config0_temp
    end

    bins_U_sigma0[:,1] = copy(samples_U[:,1])

    N_sigma0_U = bins_U_sigma0[:,1]
    config_sigma0_U = bins_U_sigma0[:,2]

    q_σ0 = zeros(config_number)
    for i=1:config_number
        idx_temp = findall(isequal(i-1),config_sigma0_U)
        q_σ0[i] = sum(N_sigma0_U[idx_temp])
    end
    q_σ0 = q_σ0/sum(q_σ0)

    return q_σ0
end


text:
using LogicCircuits, LinearAlgebra, Plots, Pkg
Pkg.activate(".")
using ProbabilisticCircuits

function learner(data :: PlainXData{Bool,BitArray{2}} , ess :: Float64)
    # Learning PSDD and CSDD
    psdd = learn_struct_prob_circuit(data; pseudocount = ess)[1]
    csdd = learn_struct_credal_circuit(WXData(data), ess)[1]
    spn = learn_probabilistic_circuit(data; pseudocount = ess)
    cspn =learn_credal_circuit(WXData(data), ess)
    @show length(psdd), length(csdd)
    @show length(spn), length(cspn)
    return true
end

#for db in twenty_dataset_names # ["nltcs"]
db = "nltcs"
@show db
ess = 1.0
threshold = 0.1
dataSet=dataset(twenty_datasets(db); do_shuffle=false, batch_size=-1)

learner(dataSet.train,ess)
dataSetMissing = add_missingness(dataSet.train, threshold)
csdd = learn_struct_credal_circuit(WXData(dataSet.train), ess)[1]

n = 16
@show log_marginal_upper(csdd, PlainXData(convert.(Int8, -ones(n,n) + 2*I)))

#sdd = load_logical_circuit("examples/random.sdd")
#@show typeof(sdd),length(sdd)

max_stars_repo_name: IDSIA/cjuice

max_issues_repo_name: IDSIA/cjuice

max_forks_repo_name: IDSIA/cjuice


text:
function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{T}}) where {T<:Robust32}
    return Rob32(Float64(rand(rng, Float32)))
end

function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{Complex{T}}}) where {T<:Robust32}
    re = rand(rng, Robust32)
    im = rand(rng, Robust32)
    return Complex{Robust32}(re, im)
end

function randpm(::Type{Robust32})
    r = rand(Robust32)
    r = rand(Bool) ? r : -r
    return r
end

function randpm(rng::MersenneTwister, ::Type{Robust32})
    r = rand(rng, Robust32)
    r = rand(rng, Bool) ? r : -r
    return r
end

function randpm(::Type{Complex{Robust32}})
    re = randpm(Robust32)
    im = randpm(Robust32)
    return Complex{Robust32}(re, im)
end

function randpm(rng::MersenneTwister, ::Type{Complex{Robust32}})
    re = randpm(rng, Robust32)
    im = randpm(rng, Robust32)
    return Complex{Robust32}(re, im)
end

function randpm(::Type{Robust32}, n::Int)
    rs = rand(Robust32, n)
    sgns = ones(n) .- (2 .* rand(Bool, n))
    return rs .* sgns
end

function randpm(::Type{Complex{Robust32}}, n::Int)
    res = rand(Robust32, n)
    ims = rand(Robust32, n)
    o = ones(n)
    rsgns = o .- (2 .* rand(Bool,n))
    isgns = o .- (2 .* rand(Bool,n))
    res = res .* rsgns
    ims = ims .* isgns
    return map((re,im)->Complex{Robust32}(re,im), res, ims)
end

function randpm(rng::MersenneTwister, ::Type{Robust32}, n::Int)
    rs = rand(rng, Robust32, n)
    sgns = ones(n) .- (2 .* rand(rng, Bool, n))
    return rs .* sgns
end

function randpm(rng::MersenneTwister, ::Type{Complex{Robust32}}, n::Int)
    res = rand(rng, Robust32, n)
    ims = rand(rng, Robust32, n)
    o = ones(n)
    rsgns = o .- (2 .* rand(rng, Bool,n))
    isgns = o .- (2 .* rand(rng, Bool,n))
    res = res .* rsgns
    ims = ims .* isgns
    return map((re,im)->Complex{Robust32}(re,im), res, ims)
end

# normal variates

function randn(rng::AbstractRNG, ::Type{Robust32})
    urand1, urand2 = rand(rng, Float64, 2)
    urand1 = urand1 + urand1 - 1
    urand2 = urand2 + urand2 - 1
    s = urand1*urand1 + urand2*urand2

    while s >= 1 || s === 0
        urand1, urand2 = rand(rng, Robust32, 2)
        urand1 = urand1 + urand1 - 1
        urand2 = urand2 + urand2 - 1
        s = urand1*urand1 + urand2*urand2
    end

    s = sqrt( -log(s) / s )
    s = (urand1 + urand2) * s
    return Robust32(s)
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

# TODO: i could use StaticArrays.jl for a more fixed type definition
using DocStringExtensions;

"""
    nrbmak(coefs::Union{Matrix{F}, Array{F,3}},
           knots::Union{Vector{F}, Vector{Vector{F}}}
          )::NURBS where {F<:AbstractFloat}

Construct the NURBS structure given the control points and the knots.

This function is used as a convenient means of constructing the NURBS data
structure. Many of the other functions in the toolbox rely on the NURBS
structure been correctly defined. This function not only constructs the proper
structure, but also checks for consistency. The user is still free to build his
own structure, in fact a few functions in the toolbox do this for convenience.

# Arguments:
- `coefs`: Control points, these can be either Cartesian or homogeneous
  coordinates.
  For a curve the control points are represented by a matrix of size `(dim,nu)`
  and for a surface by a multidimensional array of size `(dim,nu,nv)`. Where
  `nu` is number of points along the parametric U direction, and `nv` the
  number of points along the V direction. `dim` is the dimension valid options
  are:

    2 : (x,y)        2D Cartesian coordinates
    3 : (x,y,z)      3D Cartesian coordinates
    4 : (wx,wy,wz,w) 4D homogeneous coordinates

- `knots`: Non-decreasing knot sequence spanning the interval `[0.0;1.0]`. It
  is assumed that the curves and surfaces are clamped to the start and end
  control points by knot multiplicities equal to the spline order. For a curve
  knots form a vector and for a surface the knots are stored by two vectors for
  U and V in a vector of vectors `[uknots, vknots]`

# Output:
- `nurbs`: NURBS data structure see [`NURBS1D`](@ref) and [`NURBS2D`](@ref).

!!! info "Note:"
    The control points are always converted and stored within the NURBS
    structure as 4D homogeneous coordinates. For a surface the spline degree is
    a vector `[du;dv]` containing the order along the U and V directions
    respectively.

# Examples:
Construct a 2D line from `[0.0;0.0]` to `[1.5;3.0]`. For a straight line a
spline of order 2 is required. Note that the knot sequence has a multiplicity
of 2 at the start `[0.0;0.0]` and end `[1.0;1.0]` in order to clamp the ends.

```julia
julia> line = nrbmak([0.0 1.5; 0.0 3.0],vec([0.0 0.0 1.0 1.0]))
```

Construct a surface in the x-y plane i.e

    ^  (0.0,1.0) ------------ (1.0,1.0)
    |      |                      |
    | V    |                      |
    |      |      Surface         |
    |      |                      |
    |      |                      |
    |  (0.0,0.0) ------------ (1.0,0.0)
    |
    |------------------------------------>
                                      U

```julia
julia> coefs = cat([0.0 0.0; 0.0 1.0],[1.0 1.0; 0.0 1.0],dims=3);
julia> knots = [vec([0.0 0.0 1.0 1.]),vec([0.0 0.0 1.0 1.0])];
julia> plane = nrbmak(coefs,knots)
```

# Reference:
- related data types: [`NURBS1D`](@ref) and [`NURBS2D`](@ref)
"""
function nrbmak(coefs::Union{Matrix{F}, Array{F,3}},
                knots::Union{Vector{F}, Vector{Vector{F}}}
               )::NURBS where {F<:AbstractFloat}
# TODO: the matlab version does not check for open NURBS
# TODO: the matlab version does no further knot checking it only sorts and
#       normalizes
# TODO: have a look if [vec1,vec2] gives always vec([vec1,vec2])
# TODO: nrbmak does not check the m=p+1 start and end knot requirement which
#       could lead to a infinite loop at findspan
# TODO: nrbmak does not check for the m<=p requirement for central knots

nform = "B-NURBS";
ndim = 4;
np = size(coefs);
dim = np[1];
if typeof(knots)<:Vector{Vector{F}}
  if !(typeof(coefs)<:Array{F,3})
    throw(ArgumentError("coefs for surface has to be a Array{F,3}!"));
  end # if
  # constructing a surface
  number = [np[2];np[3]];
  if dim < 4
    ncoefs = zeros(F,(4,np[2],np[3]));
    ncoefs[4,:,:] .= 1.0;
    ncoefs[1:dim,:,:] = coefs;
  else
    ncoefs = coefs;
  end # if
  uorder = length(knots[1]) - np[2];
  vorder = length(knots[2]) - np[3];
  uknots = sort(knots[1]);
  vknots = sort(knots[2]);
  uknots .= (uknots .- uknots[1]) ./ (uknots[end] - uknots[1]);
  vknots .= (vknots .- vknots[1]) ./ (vknots[end] - vknots[1]);
  nknots = [uknots,vknots];
  norder = [uorder;vorder];
  return NURBS2D{typeof(ndim),F}(nform,ndim,number,norder,nknots,ncoefs);
else
  if !(typeof(coefs)<:Matrix{F})
    throw(ArgumentError("coefs for curve has to be a Matrix{F}!"));
  end # if
  # constructing a cruve
  number = np[2];
  if dim < 4
    ncoefs = zeros(F,(4,np[2]));
    ncoefs[4,:] .= 1.0;
    ncoefs[1:dim,:] = coefs;
  else
    ncoefs = coefs;
  end # if
  norder = length(knots) - np[2];
  nknots = sort(knots);
  nknots = (nknots .- knots[1]) ./ (knots[end] - knots[1]);
  return NURBS1D{typeof(ndim),F}(nform,ndim,number,norder,nknots,ncoefs);
end # if
end# nrbmak

"""
    NURBS{I<:Integer,F<:AbstractFloat}

Abstract supertype for NURBS data tayps see [`NURBS1D`](@ref) and
[`NURBS2D`](@ref).

Both curves and surfaces where represented by a structure that where compatible
with the Spline Toolbox from Mathworks.
"""
abstract type NURBS{I<:Integer,F<:AbstractFloat} end

"""
    NURBS1D{I<:Integer,F<:AbstractFloat} <: NURBS{I,F}

NURBS datatype representing a NURBS curve. Constructed with [`nrbmak`](@ref).

!!! info "Note:"
    The control points are always converted and stored within the NURBS
    structure as 4D homogeneous coordinates.

# Field
$(DocStringExtensions.TYPEDFIELDS)
"""
struct NURBS1D{I<:Integer,F<:AbstractFloat} <: NURBS{I,F}
"Type name 'B-NURBS'"
form::String
"Dimension of the control points"
dim::I
"Number of Control points"
number::I
"Order of the spline"
order::I
"Knot sequence"
knots::Vector{F}
"Control Points"
coefs::Matrix{F}
end # NURBS1D

"""
    NURBS2D{I<:Integer,F<:AbstractFloat} <: NURBS{I,F}

NURBS datatype representing a NURBS surface. Constructed with [`nrbmak`](@ref).

!!! info "Note:"
    The control points are always converted and stored within the NURBS
    structure as 4D homogeneous coordinates. For a surface the spline degree is
    a vector `[du;dv]` containing the order along the U and V directions
    respectively.

# Field
$(DocStringExtensions.TYPEDFIELDS)
"""
struct NURBS2D{I<:Integer,F<:AbstractFloat} <: NURBS{I,F}
"Type name 'B-NURBS'"
form::String
"Dimension of the control points"
dim::I
"Number of Control points"
number::Vector{I}
"Order of the spline"
order::Vector{I}
"Knot sequence"
knots::Vector{Vector{F}}
"Control Points"
coefs::Array{F,3}
end # NURBS2D

# pretty printing NURBS1D
Base.show(io::IO, ::MIME"text/plain", z::NURBS1D{I,F}) where{I,F} =
  print(io, "NURBS1D{$I, $F} curve nurbs structre:\n   ", z);
Base.show(io::IO, z::NURBS1D) = print(io,"order: ",z.order ,", ",
                                      "#points: ", z.number);

# pretty printing NURBS2D
Base.show(io::IO, ::MIME"text/plain", z::NURBS2D{I,F}) where{I,F} =
  print(io, "NURBS2D{$I, $F} surface nurbs structre:\n   ", z);
Base.show(io::IO, z::NURBS2D) = print(io,"order: ",z.order ,", ",
                                      "#points: ", z.number);


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
gamma = 5/3
equations = IdealGlmMhdEquations1D(gamma)

"""
    initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Extended version of the test of Shu and Osher for one dimensional ideal MHD equations.
Taken from Section 4.1 of
- Derigs et al. (2016)
  A Novel High-Order, Entropy Stable, 3D AMR MHD Solver withGuaranteed Positive Pressure
  [DOI: 10.1016/j.jcp.2016.04.048](https://doi.org/10.1016/j.jcp.2016.04.048)
"""
function initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
  # initial shock location is taken to be at x = -4
  x_0 = -4.0
  rho = x[1] <= x_0 ? 3.5 : 1.0 + 0.2 * sin(5.0 * x[1])
  v1 = x[1] <= x_0 ? 5.8846 : 0.0
  v2 = x[1] <= x_0 ? 1.1198 : 0.0
  v3 = 0.0
  p = x[1] <= x_0 ? 42.0267 : 1.0
  B1 = 1.0
  B2 = x[1] <= x_0 ? 3.6359 : 1.0
  B3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D)

Extended version of the test of Shu and Osher for one dimensional ideal MHD equations
but shock propogates from right to left.

!!! note
    This is useful to exercise some of the components of the HLL flux.
"""
function initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
  # initial shock location is taken to be at x = 4
  x_0 = 4.0
  rho = x[1] <= x_0 ? 1.0 + 0.2 * sin(5.0 * x[1]) : 3.5
  v1 = x[1] <= x_0 ? 0.0 : -5.8846
  v2 = x[1] <= x_0 ? 0.0 : -1.1198
  v3 = 0.0
  p = x[1] <= x_0 ? 1.0 : 42.0267
  B1 = 1.0
  B2 = x[1] <= x_0 ? 1.0 : 3.6359
  B3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end
initial_condition = initial_condition_shu_osher_shock_tube

boundary_conditions = BoundaryConditionDirichlet(initial_condition)

surface_flux = flux_hll
volume_flux  = flux_hindenlang_gassner
basis = LobattoLegendreBasis(4)

indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = -5
coordinates_max =  5
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=10_000,
                periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.7)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(energy_kinetic, energy_internal,
                                                               energy_magnetic, cross_helicity))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorHennemannGassner(semi,
                                          alpha_max=0.5,
                                          alpha_min=0.001,
                                          alpha_smooth=true,
                                          variable=density_pressure)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=4,
                                      max_level=7, max_threshold=0.01)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback, stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
"""
ptcBeam(n, maxit, delta=.01, lambda=20.0; precision=Float64, keepsolhist=false)

Test PTC for systems on the buckling beam problem.
Compare to Newton, which will converge to the unstable solution.
"""
function ptcBeam(
    n,
    maxit,
    delta = 0.01,
    lambda = 20.0;
    precision = Float64,
    keepsolhist = false,
    jknowsdt = false,
)
    #
    # This is a steady-state computation, so there is no dt in the problem.
    #
    bdata = beaminit(n, 0.0, lambda)
    x = bdata.x
    u0 = x .* (1.0 .- x) .* (2.0 .- x)
    u0 .*= exp.(-10.0 * u0)
    FS = copy(u0)
    FPS = precision.(copy(bdata.D2))
    if jknowsdt
        Jeval = BeamJdt!
    else
        Jeval = BeamJ!
    end
    bout = ptcsol(
        FBeam!,
        u0,
        FS,
        FPS,
        Jeval;
        #        BeamJ!;
        rtol = 1.e-10,
        pdata = bdata,
        delta0 = delta,
        maxit = maxit,
        jknowsdt = jknowsdt,
        keepsolhist = keepsolhist,
    )
    if ~jknowsdt
        qout = nsol(FBeam!, u0, FS, FPS, BeamJ!; pdata = bdata, sham = 1)
        return (bout, qout)
    else
        return bout
    end
end

"""
BeamJdt!(FP, FV, U, dt, bdata)
Evaluates the Jacobian + (1/dt) I for PTC.
"""
function BeamJdt!(FP, FV, U, dt, bdata)
    FP .= BeamJ!(FP, FV, U, bdata)
    FP .= FP + (1.0 / dt) * I
end


text:
const DEFAULT_SOURCE = 'F'
const EOL_DEFAULT = "\n"
const TZ = tz"UTC"
const DEFAULT_STORE_ALL_RECORDS = false

const SINGLE_INSTANCE_DATA_RECORDS_CHARS = Set(['A', 'H', 'I', 'J', 'C', 'D', 'G'])
const MULTIPLE_INSTANCE_DATA_RECORDS_CHARS = Set(['B', 'E', 'F', 'K', 'L'])
const DATA_RECORDS_CHARS = union(SINGLE_INSTANCE_DATA_RECORDS_CHARS, MULTIPLE_INSTANCE_DATA_RECORDS_CHARS)

module FixValidity
    @enum FixValidityEnum begin
        Fix3D = Int('A')  # 3D fix
        Fix2D = Int('V')  # 2D fix (no GNSS altitude) or for no GNSS data (pressure altitude data)
    end

    const ALLOWED = Set(Char(Int(fix)) for fix in instances(FixValidityEnum))

    function parse(val::Char)
        val ∉ ALLOWED && throw(ArgumentError("Invalid fix validity $(repr(val)). It can only be in $(repr(ALLOWED))"))
        return FixValidityEnum(Int(val))
    end

    function string(fv::FixValidityEnum)
        return "$(Char(Integer(fv)))"
    end

    function write(stream::IO, fv::FixValidityEnum)
        Base.write(stream, string(fv))
    end
end

module GpsQualifier
    @enum GpsQualifierEnum begin
        GPS = 1
        DGPS = 2
    end

    const ALLOWED = Set(Int(fix) for fix in instances(GpsQualifierEnum))

    function decode(val)
        val ∉ ALLOWED && throw(ArgumentError("Invalid GPS qualifier $(repr(val)). It can only be in $(repr(ALLOWED))"))
        return GpsQualifierEnum(val)
    end

    function parse(val::Char)
        return decode(Base.parse(Int64, val))
    end

    function string(qualifier::GpsQualifierEnum)
        return "$(Integer(qualifier))"
    end

    function write(stream::IO, qualifier::GpsQualifierEnum)
        Base.write(stream, "$(string(qualifier))")
    end
end

# ===

module ParsingMode
    @enum(ParsingModeEnum, STRICT, UNSTRICT)

    const DEFAULT = UNSTRICT
end


text:
module GridapMHD

using Random
using SparseArrays
using SparseMatricesCSR
using PartitionedArrays
using Gridap
using Gridap.Helpers
using Gridap.Algebra
using Gridap.CellData
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using GridapDistributed
using GridapPETSc
using GridapPETSc.PETSC
using FileIO
using BSON

include("Main.jl")

include("Hunt.jl")

end # module


text:
#
# Discrete Uniform Sampler
#
# parameters:
# 	(1) a (lower bound)
#	(2) b (upper bound)
#
const ONE_EIGTH_INT_RANGE = (typemax(Int)>>3) - (typemin(Int)>>3) # this formulation avoids overflow

type DiscreteUniformSampler <: DiscreteDistributionSampler
	paramranges::Vector{Tuple{Float64,Float64}}
	params::Vector{Float64}
	distribution::Distribution
	function DiscreteUniformSampler(params=Float64[])
		intbounds = (Float64(typemin(Int)), Float64(typemax(Int)))
		s = new([intbounds, intbounds])
		setparams!(s, isempty(params) ? [Float64(typemin(Int)), Float64(typemax(Int))] : params)
		s
	end
end

function setparams!(s::DiscreteUniformSampler, params)
	# note the paramranges check will trap an attempt to sample from the entire UInt range: this cannot be supported
	# since Distributions.Uniform currently only supports Int range (and then only with the workaround below using a Mixture Model)
	# therefore adjust parameters and warn
	adjparams = deepcopy(params)
	for i in 1:length(adjparams)
		if adjparams[i] < Float64(typemin(Int))
			warn("Adjusting parameters of discrete uniform sampler to typemin(Int)")
			adjparams[i] = Float64(typemin(Int))
		end
		if adjparams[i] > Float64(typemax(Int))
			warn("Adjusting parameters of discrete uniform sampler to typemax(Int)")
			adjparams[i] = Float64(typemax(Int))
		end
	end
	checkparamranges(s, adjparams)
	# swap parameters if necessary (as a silent repair during search)
	if adjparams[1] <= adjparams[2]
		a, b = adjparams[[1,2]]
	else
		a, b = adjparams[[2,1]]
	end
	# convert to Int, but with a check to avoid InexactError when at bounds since params will be Float64 and imprecision means
	# that the values can be outside range
	a = a < typemin(Int) ? typemin(Int) : round(Int, a, RoundUp) # round up so that value sampled remain in unrounded range
	b = b > typemax(Int) ? typemax(Int) : round(Int, b, RoundDown) # round down so that value sampled remain in unrounded range
	# if range is too large (more than one quarter of entire Int range), then must use Mixture
	# distribution instead to avoid InexactErrors when sampling from DiscreteUniform
	halfrange = (b>>1) - (a>>1) # >> 1 formulation avoids overflow by calculating half range
	if  halfrange < ONE_EIGTH_INT_RANGE
		s.distribution = DiscreteUniform(a, b)
	else
		# in the worst case (the entire Int range), we must split into four ranges for the
		# mixture model if we are to avoid errors, so let's always do this
		quarterrange = halfrange >> 1  # one quarter desired range
		q1 = a + quarterrange
		q2 = q1 + quarterrange
		q3 = q2 + quarterrange
		prior = [Float64(q1)-Float64(a)+1, Float64(q2)-Float64(q1), Float64(q3)-Float64(q2), Float64(b)-Float64(q3)]
		s.distribution = MixtureModel(
			[DiscreteUniform(a,q1), DiscreteUniform(q1+1,q2), DiscreteUniform(q2+1,q3),  DiscreteUniform(q3+1,b)], prior ./ sum(prior))
	end
	s.params = [Float64(a), Float64(b)]
end

getparams(s::DiscreteUniformSampler) = s.params

estimateparams!(s::DiscreteUniformSampler, traces) = estimateparams!(s, traces, 2)


text:
# Astronomy example from Hairer et al 1992 p.129 (see also Szebehely
# 1967): the orbit of a satellite in planet-moon system.  This system
# has know periodic solutions, so called “Arenstorf orbits”
# (Arenstorf, 1963); one of which is used as a test-case.
#
# This is a non-stiff problem.  Adaptive solvers will fare much better
# as there are three regions where the satellite experiences very
# strong acceleration interspersed with regions of very small
# acceleration.

export threebody

threebody = let
    tcname = :threebody
    T = Float64 # the datatype used
    Tarr = Matrix
    dof = 4 # degrees of freedom
    dae=0  # index of DAE, ==0 for ODE
    # stiffness of system, one of the three constants
    stiffness = [nonstiff, mildlystiff, stiff][1]

    # parameters
    μ = 0.012277471
    μ′ = 1-μ

    function fn!(t,y,dydt)
        # the ode function

        y1, y2, y1′, y2′ = y[:] # this allocates memory... but less than below
        D1 = ((y1+μ )^2 + y2^2)^(3/2)
        D2 = ((y1-μ′)^2 + y2^2)^(3/2)
        y1″ = y1 + 2*y2′ - μ′*(y1+μ)/D1 - μ*(y1-μ′)/D2
        y2″ = y2 - 2*y1′ - μ′*y2/D1     - μ*y2/D2
        dydt[:] = [y1′, y2′, y1″, y2″] # this allocates memory...
        return nothing
    end
    fn!( ; T_::Type=T, dof_=dof) = zeros(T_,dof_)

    jac! = nothing
    mass! = nothing

    ic = T[0.994, 0.0, 0.0, -2.00158510637908252240537862224] # vector of initial conditions
    tspan = T[0, 17.0652165601579625588917206249] # integration interval
    refsol = ic # reference solution at tspan[2]
    refsolinds = trues(dof)
    scd_absinds = ic.==0.0
    tc = TestCaseExplicit{tcname, T, Tarr}(
                             stiffness,
                             dae,
                             dof,
                             fn!,
                             jac!,
                             mass!,
                             ic,
                             tspan,
                             refsol,
                             refsolinds,
                             scd_absinds)

    # put into the right buckets
    tc_all[tcname] = tc
    tc_nonstiff[tcname] = tc
end


text:
#Neural net in julia
module mlp
#setting parameters
n_hidden = nothing
lr = nothing
max_iter = nothing
w0 = nothing
w1 = nothing
error_vector = nothing

function set(n_hidden_=3, lr_=0.4, max_iter_=500, seed=nothing)
    global n_hidden = n_hidden_
    global lr = lr_
    global max_iter = max_iter_
    if seed!=nothing
        srand(seed)
    end
end
#defining derivative, using tanh as tf
derivative(x) = 1 - (tanh.(x).^2)

#defining feedforward part
function feedforward(X)
    global w0
    global w1
    #println( size(X), size(w0), size(w1))

    hidden = tanh.(X*w0)
    #println("Hidden shape", size(hidden))
    output = tanh.(hidden*w1)
    #println("Output shape", size(output))
    return [hidden, output]
end

#defining the predict function
predict(X) = feedforward(X)[2]
#shape corrector
function y_reshaper(y)
    if size(y)[2] > size(y)[1]
        return y'
    else
        return y
    end
end

#training
function fit!(X, y)

    y = y_reshaper(y)
    global w0 = randn(size(X)[2], n_hidden)
    global w1 = randn(n_hidden, size(y)[2])
    #println(size(w0))
    #println(size(w1))
    global max_iter
    global n_hidden
    # training
    #error vector
    global error_vector = ones(max_iter)
    println("starting training")
    for i = 1:max_iter

        h, o = feedforward(X)
        o_error = y - o
        o_delta = lr * o_error .* derivative(o)

        # propagating error
        h_error = o_delta * w1'
        h_delta = lr * h_error .* derivative(h)

        err = sum(abs.(o_error))
        global error_vector[i] = err
        if (i%200)==0
            println("Accumulated error: ", err)
        end

        #updating weights
        global w0 += (X' * h_delta)
        global w1 += (h' * o_delta)
    end
    println("OVER")
end

end

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
using mlp
using PyPlot

X = [0 0; 0 1; 1 0; 1 1]
y = [0 1 1 0]

mlp.set(3, 0.4, 500, 2)
mlp.fit!(X,y)

println("\n\n%%%%%%%results%%%%%%%")
println("Total iterations: ", mlp.max_iter, "| Learning rate: ", mlp.lr, "| Nodes in the hidden: ",mlp. n_hidden)
results = mlp.predict(X)
#println(results)
println("Final accumulated error: ", sum(abs.(y'-results)))
for i=1:4
    println(X[i,:], " : ", results[i, :])
end
#plot

plot(mlp.error_vector)
title("Error x iteration")
ylabel("Accumalated error")
xlabel("Iteartion")
show()


text:
#=
Popis
  Funkcia vrati podgraf s vrcholmi ktore patria do k_shell.

Parametre
  g : datova struktura grafu :: GenericGraph
  k : stupen k_shell :: Int64
  core_numbers : slovnik obsahujuci indexy vrcholov a ich hodnotu jadra :: Dict{Int64,Int64}()

Navratova hodnota
  : k_shell podgraf vstupneho grafu :: GenericGraph

Poznamky
  Hodnota core_numbers je volitelna, ak nieje zadana tak prislusnost jednotlivych vrcholov
  do jadier je vypocitana funkciou coreness() zo vstupneho grafu.

  K_shell je podgraf, v ktorom vrcholi maju hodnotu jadra == k.
=#

function k_shell(g::GenericGraph; k=nothing, core_numbers::Dict{Int64,Int64}=Dict{Int64,Int64}())

  if core_numbers == Dict{Int64,Int64}()
    core_numbers=coreness(g)
  end

  if k == nothing
    k=maximum(values(core_numbers))
  end

  verts = collect(keys(filter((key, value) -> value == k, core_numbers)))

  return get_subgraph(verts, g)

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
module MagneticParticles

using Optim, StaticArrays, RecipesBase

include("sw_particle.jl")
include("random_sphere_points.jl")
include("usov_particle.jl")

export SwParticle, DataMH
export obtain_hysteresis_loop, calculate_particles

export get_points_spherical_random
export get_points_cartesian_random
export get_points_spherical_evenly
export random_sphere_pack
export create_for_constant_concentration, apply_field!, draw_representation, magnetize_particle

end


text:
"""
    C_iio_device_get_sample_size(device)

Get the current sample size.

# Parameters
- `device::Ptr{iio_device}` : A pointer to an iio_device structure

# Returns
- On success, the sample size in bytes
- On error, a negative errno code is returned

# NOTE
The sample size is not constant and will change when channels get enabled or disabled.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Debug.html#ga52b3e955c10d6f962b2c2e749c7c02fb)
"""
function C_iio_device_get_sample_size(device::Ptr{iio_device})
    return ccall(
        (:iio_device_get_sample_size, libIIO),
        Cssize_t, (Ptr{iio_device},),
        device
    );
end

"""
    C_iio_device_identify_filename(device, filename)

Identify the channel or debug attribute corresponding to a filename.

# Parameters
- `device::Ptr{iio_device}` : A pointer to an iio_device structure
- `filename::String` : A NULL-terminated string corresponding to the filename

# Returns
- On success, `(0, channel::Ptr{iio_channel}, attribute::String)` is returned.
- On error, `(errno, NULL, "")` is returned, where errno is a negative error code.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Debug.html#ga87ef46fa578c7be7b3e2a6f9f16fdf7e)
"""
function C_iio_device_identify_filename(device::Ptr{iio_device}, filename::String)
    channel = Ref{Ptr{iio_channel}}();
    attribute = Ref{Ptr{Cuchar}}();
    ret = ccall(
        (:iio_device_identify_filename, libIIO),
        Cint, (Ptr{iio_device}, Cstring, Ptr{Ptr{iio_channel}}, Ptr{Ptr{Cuchar}}),
        device, filename, channel, attribute
    );
    return attribute[] == C_NULL ? (ret, channel[], "") : (ret, channel[], Base.unsafe_string(attribute[]));
end


text:
"Miscellaneous utilities that are generally useful"
module MiscUtils

export @def

"""
Insert source code macro.

Produces another macro that you may use to do the insertion.

Easy code reuse; I use for type definitions with shared fields (like concrete type
inheritance).

Copied from `@def` in Optim.jl.
<https://discourse.julialang.org/t/some-clarifications-and-best-practices-on-the-oo-model-in-julia/693/8>

**Warning** Can lead to spaghetti code a la FORTRAN GOTO if used poorly.
"""
macro def(name, definition)
  esc(quote
    macro $name()
      esc($(Expr(:quote, definition)))
    end
  end)
end

end # module


text:
"""
A Variational Bayesian modeling approach. Models are defined by
bipartite graphs in which factors defining terms in the variational
objective are connected to nodes defining variables.
"""
module VinDsl

using Distributions
using PDMats
using StatsFuns
using Optim
using ForwardDiff
import Base: convert, call, show, zero, getindex, setindex!, ndims
import Distributions: params, var, entropy, cov, invcov, logpdf, logdetcov, logmvgamma
using Base.Cartesian

# general helper functions
include("utils.jl")

# new distributions and extended functionality
include("distributions.jl")

# type system
include("types.jl")

# operations and structure on types
include("structure.jl")

# inference methods in models
include("inference.jl")

# tools for working with expressions
include("expressions.jl")

# random variable calculus
include("expectations.jl")

# language macros
include("dsl.jl")

# predefined factors
include("factors.jl")

# predefined conjugate updates
include("conjugate_relations.jl")

export VBModel,
    RVType, RScalar, RVector, RMatrix, RReal, RPositive, RNegative, RBounded, RProbability, RCorrelation,
    RRealVec, RUnitVec, ROrdered, RPosOrdered, RSimplex, RCholFact, RCholCorr, RCorrMat, RCovMat, RCovLKJ,
    @advi_declarations, @advi_model, @ELBO,
    constrain, unconstrain, logdetjac, supp, parsupp, nfree,
    Factor, FactorInds, @factor, get_structure, project, project_inds, @wrapvars, value, @pmodel,
    EntropyFactor, LogNormalFactor, LogGammaCanonFactor, LogMvNormalCanonFactor,
    LogGammaFactor,
    LogMvNormalDiagCanonFactor, LogMarkovChainFactor, LogDirichletFactor,
    LogMarkovMatrixFactor, LogWishartFactor,
    _simplify, _simplify_inside, _simplify_call, _simplify_compose, get_all_syms, get_all_inds, strip_inds, @exprnode,
    @deffactor, @~,
    Node, RandomNode, ConstantNode, ExprNode, ExprDist, @simplify, nodeextract,
    register, check_conjugate, update!, unroll_pars, update_pars!, reroll_pars,
    get_par_sizes, flatten,
    E, Elog, Eloggamma, Elogmvgamma, Elogdet, V, H, C, value, naturals, @defnaturals, logpdf,
    get_node_size, get_name_mapping, HMM, MarkovChain, MarkovMatrix,
    MatrixNormal,
    nstates, naturals, naturals_to_params, uparams, forwardbackward,
    LKJ, _cpc_to_corr, meanlogdet, meanlog, entropy

end  # module


text:
module UOT

using LinearAlgebra
using JuMP
# solver (library GLPK could be considered but slower)
using Mosek, MosekTools

function euclidian(x::Array{Float64,1}, y::Array{Float64,1})
    return norm(x - y);
end

# structure to store output of compute_flat_metric
mutable struct UOT_output
    cost::Float64
    # matrix of couplings
    P::Array{Float64,2}
    # matrix of costs
    C::Array{Float64,2}
    # vector indicating creation and destruction of mass
    mass_variation::Array{Float64,1}
end

function init_UOT_output(n::Int64 = 1, m::Int64 = 1)
    return UOT_output(Inf, zeros(n, m), zeros(n, m), zeros(n + m));
end

"""
    compute_flat_metric(x, y, a, b, lambda, coupling_mass_variation = true, ground_metric = euclidian)
Compute the flat metric between two point sources distributions given by the locations `x`
with weights `a` and the locations `y` with weights `b`, for the threshold `lambda`. If
`coupling_mass_variation` is set to true then it also computes the coupling and mass variation
between the two distributions.
# Example
```jldoctest
julia> x, y, a, b = [[0.0, 0.0]], [[1.0, 1.0]], [1.0], [0.5];
julia> UOT.compute_flat_metric(x, y, a, b, 1.0)
UOT.UOT_output(1.2071067811865477, [0.5], [1.41421], [0.5, -2.49401e-12])
```
"""
function compute_flat_metric(x::Array{Array{Float64,1},1}, y::Array{Array{Float64,1},1}, a::Array{Float64,1}, b::Array{Float64,1}, lambda::Float64, coupling_mass_variation::Bool = true, ground_metric::Function = euclidian)
    n, m = length(x), length(y);
    out = init_UOT_output(n, m);
    c = cat(a, -b, dims = 1);
    p = n*m;

    # model = Model(with_optimizer(GLPK.Optimizer));
    model = Model(with_optimizer(Mosek.Optimizer, MSK_IPAR_LOG = 0)); # MSK_IPAR_NUM_THREADS number of threads https://www.gams.com/latest/docs/S_MOSEK.html

    @variable(model, f[1:n + m]);
    i = 1;
    for k in 1:n
        for l in 1:m
            e = ground_metric(x[k], y[l]);
            constr = @constraint(model, -e <= f[k] - f[l + n] <= e)
            # Lipschitz constraints of the optimization problem
            set_name(constr, "lip_constraints_$(i)")
            i += 1;
        end
    end
    # Box constraints of the optimization problem
    @constraint(model, bound_constraints, -lambda .<= f .<= lambda);

    @objective(model, Max, c'*f);
    JuMP.optimize!(model);
    out.cost = JuMP.objective_value(model);
    # println(termination_status(model))
    # println(JuMP.solve_time(model))

    # compute coupling and cost matrices from primal and dual solutions
    if coupling_mass_variation
        dual = zeros(p);
        for i in 1:p
            dual[i] = -JuMP.dual.(constraint_by_name(model, "lip_constraints_$(i)"));
        end
        out.mass_variation = -JuMP.dual.(bound_constraints);
        f_opt_value = JuMP.value.(f);

        k = 1;
        for i in 1:n
            for j in 1:m
                out.C[i, j] = f_opt_value[i] - f_opt_value[j + n];
                # to prevent numerical approximations: remove?
                if abs(dual[k]) > 1e-8
                    out.P[i, j] = dual[k];
                end
                k += 1;
            end
        end
    end
    return out;
end

end # module

max_stars_repo_name: qdenoyelle/OTMSL

max_issues_repo_name: qdenoyelle/OTMSL

max_forks_repo_name: qdenoyelle/OTMSL


text:
module CommunalHelperStationBlockTrack
using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/StationBlockTrack" StationBlockTrack(
                    x::Integer, y::Integer,
                    width::Integer = 24, height::Integer = 24,
                    horizontal::Bool = false)

const placements = Ahorn.PlacementDict(
    "Station Block Track (Vertical) (Communal Helper)" => Ahorn.EntityPlacement(
        StationBlockTrack,
        "rectangle",
        Dict{String, Any}(
            "horizontal" => false
        )
    ),
    "Station Block Track (Horizontal) (Communal Helper)" => Ahorn.EntityPlacement(
        StationBlockTrack,
        "rectangle",
        Dict{String, Any}(
            "horizontal" => true
        )
    )
)

Ahorn.minimumSize(entity::StationBlockTrack) = 24, 24

Ahorn.resizable(entity::StationBlockTrack) = Bool(get(entity.data, "horizontal", false)), !(Bool(get(entity.data, "horizontal", false)))

function Ahorn.selection(entity::StationBlockTrack)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))
    horiz = Bool(get(entity.data, "horizontal", false))

    return horiz ? Ahorn.Rectangle(x, y, width, 8) : Ahorn.Rectangle(x, y, 8, height)
end

nodeSprite = "objects/CommunalHelper/stationBlock/tracks/track/ball"
vTrack = "objects/CommunalHelper/stationBlock/tracks/track/pipeV"
hTrack = "objects/CommunalHelper/stationBlock/tracks/track/pipeH"

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::StationBlockTrack, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))

    horiz = Bool(get(entity.data, "horizontal", false))

    if horiz
        tilesWidth = div(width, 8)

        for i in 2:tilesWidth - 1
            Ahorn.drawImage(ctx, hTrack, x + (i - 1) * 8, y)
        end

        Ahorn.drawImage(ctx, nodeSprite, x, y)
        Ahorn.drawImage(ctx, nodeSprite, x + width - 8, y)
    else
        tilesHeight = div(height, 8)

        for i in 2:tilesHeight - 1
            Ahorn.drawImage(ctx, vTrack, x, y + (i - 1) * 8)
        end

        Ahorn.drawImage(ctx, nodeSprite, x, y)
        Ahorn.drawImage(ctx, nodeSprite, x, y + height - 8)
    end

end

end

max_stars_repo_name: FlynnSC/CommunalHelper

max_issues_repo_name: FlynnSC/CommunalHelper

max_forks_repo_name: FlynnSC/CommunalHelper


text:
function InterpolatingMicrosphere2D(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return InterpolatingMicrosphere2D((jint, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function copy(obj::InterpolatingMicrosphere2D)
    return jcall(obj, "copy", InterpolatingMicrosphere2D, ())
end