text:
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
using Pkg
using PyCall
using OrderedCollections
using LinearAlgebra
#using UnicodePlots
using NSGAII
using SpikeSynchrony
include("../src/SpikingNeuralNetworks.jl")
include("../src/unit.jl")
SNN = SpikingNeuralNetworks
using Evolutionary, Test, Random
using Plots
import DataStructures
using JLD
using Evolutionary, Test, Random
#Pkg.add(url="https://github.com/JuliaObjects/ConstructionBaseExtras.jl")

#Pkg.add(url="https://github.com/JuliaObjects/ConstructionBase.jl")
#using Plots

using CuArrays
using GPUifyLoops

using ConstructionBase
using ConstructionBaseExtras
using Debugger
unicodeplots()

SNN = SpikingNeuralNetworks
adparam = SNN.ADEXParameter(;a = 6.050246708405076, b = 7.308480222357973,
    cm = 803.1019662706587,
    v_rest= -63.22881649139353,
    tau_m=19.73777028610565,
    tau_w=351.0551915202058,
    v_thresh=-39.232165554444265,
    delta_T=6.37124632135508,
    v_spike = -36.58205819488362,
    v_reset = -59.18792270568965,
    spike_delta = 16.33506432689027)

E = SNN.AD(;N = 1, param=adparam)
E.I = [795.57128906]
#{'value': array(460.57128906) * pA}
# -

SNN.monitor(E, [:v,:I])
#SNN.sim!([Etest],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.sim!([E],[0],dt = 0.25ms, simulation_duration = 2000ms, delay = 500ms,stimulus_duration=2000ms)
SNN.vecplot(E, :v) |> display

py"""
from neo import AnalogSignal
"""

if isfile("ground_truth.jld")
    vmgtv = load("ground_truth.jld","vmgtv")
    ngt_spikes = load("ground_truth.jld","ngt_spikes")
    gt_spikes = load("ground_truth.jld","gt_spikes")

    ground_spikes = gt_spikes
    ngt_spikes = size(gt_spikes)[1]
    vmgtt = load("ground_truth.jld","vmgtt")

    #plot(Plot(vmgtv,vmgtt,w=1))

else

    py"""
    options = (
        325479788,
        324257146,
        476053392,
        623893177,
        623960880,
        482493761,
        471819401
    )
    specimen_id = options[1]
    target_num_spikes=7
    from neuronunit.allenapi import make_allen_tests_from_id
    sweep_numbers, data_set, sweeps = make_allen_tests_from_id.allen_id_to_sweeps(specimen_id)
    (vmm,stimulus,sn,spike_times) = make_allen_tests_from_id.get_model_parts_sweep_from_spk_cnt(
        target_num_spikes, data_set, sweep_numbers, specimen_id
    )
    """
    gt_spikes = py"spike_times"
    ground_spikes = gt_spikes

    ngt_spikes = size(gt_spikes)[1]
    vmgtv = py"vmm.magnitude"
    vmgtt = py"vmm.times"

    save("ground_truth.jld", "vmgtv", vmgtv,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes)

end

#plot(vmgtt[:],vmgtv[:])
plot(vmgtt[:],vmgtv[:]) |> display

ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms
#break_on(:error)
function get_ranges()
    ranges = DataStructures.OrderedDict{Char,Float32}()
    ranges = ("a"=>[0.002,0.3],"b"=>[0.02,0.36],"c"=>[-75,-35],"d"=>[0.005,16])#,"I"=>[100,9000])
    lower = []
    upper = []
    for (k,v) in ranges
        append!(lower,v[1])
        append!(upper,v[2])
    end
    lower,upper
end

function init_b(lower,upper)
    gene = []
    #chrome = Float32[size(lower)[1]]
    for (i,(l,u)) in enumerate(zip(lower,upper))
        p1 = rand(l:u, 1)
        append!(gene,p1)
        #chrome[i] = p1
    end
    gene
end

function initf(n)
    genesb = []
    for i in 1:n
        genes = init_b(lower,upper)
        append!(genesb,[genes])
    end
    genesb
end
lower,upper = get_ranges()

function raster_synchp(p)
    fire = p.records[:fire]
    spikes = Float32[]
    #neurons = Float32[]#, Float32[]
    for time = eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(spikes,time)
        end
    end
    spikes
end
function loss(E,ngt_spikes,ground_spikes)
    spikes = raster_synchp(E)
    spikes = [s/1000.0 for s in spikes]
    maxt = findmax(sort!(unique(vcat(spikes,ground_spikes))))[1]
    if size(spikes)[1]>1
        t, S = SPIKE_distance_profile(spikes, ground_spikes;t0=0,tf = maxt)
        spkd = SpikeSynchrony.trapezoid_integral(t, S)/(t[end]-t[1]) # == SPIKE_distance(y1, y2)
    else
        spkd = 1.0
    end
    #delta = size(spikes)[1] - ngt_spikes
    #=
    if delta == 0
        #println()
        #println(" hit ")
        comp = 10.0
    else
        #println(abs(delta))
        comp = (1.0/abs(delta))
    end
    =#
    #spkd = 1.0-spkd
    #err = spkd #+ comp
    return spkd#+spkd

end
#=
function ConstraintVioloation(x)
    if size(spikes)[1] == size(ground_spikes)[1]
        return 0
    else
        return size(spikes)[1] + size(ground_spikes)[1]
    end
end
=#
function test_current(param,current,ngs)
    #println(param)
    CH = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(;N = 1, param = CH)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    E.I = [current*nA]#[param[5]*nA]
    SNN.sim!([E], []; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+500ms)
    #v = SNN.getrecord(E, :v)
    #SNN.vecplot(E, :v) |> display

    spikes = raster_synchp(E)
    spikes = [s/1000.0 for s in spikes]
    #println(spikes)
    nspk = size(spikes)[1]
    delta = nspk - ngs
    #println(nspk)
    #if delta > 0
    #    less_current = true
    #else
    #    less_current = false
    #end
    #println(delta)
    return delta
end
function test_c(check_values,param,ngt_spikes,current_dict)
    nspk = -1.0
    for i_c in check_values
        nspk = test_current(param,i_c,ngt_spikes)
        current_dict[nspk] = i_c
        if nspk == 0
            if 0 in keys(current_dict)
                return current_dict,0
            end
        end
    end
    return current_dict,nspk
end

function current_search(param)

    current_dict = Dict()
    minc = 0.0
    maxc = 90000.0
    step_size = (maxc-minc)/10.0
    check_values = minc:step_size:maxc
    current_dict = Dict()
    cnt = 0
    while (0 in keys(current_dict))==false

        current_dict,nspk = test_c(check_values,param,ngt_spikes,current_dict)
        over_s = Dict([(k,v) for (k,v) in current_dict if k>0])
        under_s = Dict([(k,v) for (k,v) in current_dict if k<0])
        if length(over_s)>0
            # find the lowest current that caused more than one spike
            # throw away part of dictionary that has no spikes
            # find minimum of value in the remaining dictionary
            new_top = findmin(collect(values(over_s)))[1]
            #new_bottom = minc-minc*(1.0/2.0)

        else
            new_top = maxc*2

        end
        if length(under_s)>0
            #new_top = maxc*2
            # find the lowest current that caused more than one spike
            # throw away part of dictionary that has no spikes
            # find minimum of value in the remaining dictionary
            new_bottom = findmax(collect(values(under_s)))[1]

        else
            new_bottom = minc-minc*(1.0/2.0)
        end

        step_size = abs(new_top-new_bottom)/10.0
        if step_size==0
            return findmin(collect(values(current_dict)))[1]
        end
        check_values = new_bottom:step_size:new_top

        cnt+=1
        if cnt >140
            return findmin(collect(values(current_dict)))[1]
        end
    end
    return current_dict[0]
end

function zz_(param)

    #param = x
    #println(param)
    #bincoded = NSGAII.encode(seed, bc)
    #NSGAII.decode!(bincoded, bc)
    #bc = BinaryCoding(6, [:Float32,:Float32,:Float32,:Float32,:Float32], lower, upper)

    #bc = BinaryCoding(size(lower)[1], lower, upper)
    #println(param)
    CH = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(;N = 1, param = CH)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    #T = (ALLEN_DURATION+ALLEN_DELAY)*ms
    E.I = [current_search(param)*nA]
    #print(param[5])
    #print(param)
    #println("fail")
    SNN.sim!([E], []; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+443ms)
    #v = SNN.getrecord(E, :v)
    #SNN.vecplot(E, :v) |> display
    #function loss(E,ngt_spikes,ground_spikes)

    error = loss(E,ngt_spikes,ground_spikes)
    #println(error)
    error
end

function z(x::AbstractVector)
    param = x
    CH = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(;N = 1, param = CH)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    E.I = [param[5]*nA]
    SNN.sim!([E], []; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+500ms)
    error = loss(E,ngt_spikes)
    error
end

#function rastrigin(x::AbstractVector{T}) where {T <: AbstractFloat}
#    n = length(x)
#    return 10n + sum([ x[i]^2 - 10cos(convert(T,2π*x[i])) for i in 1:n ])
#end

#function rastrigin(x::AbstractVector{T}) where {T <: AbstractFloat}
#    n = length(x)
#    return 10n + sum([ x[i]^2 - 10cos(convert(T,2π*x[i])) for i in 1:n ])
#end

function checkmodel(param)
    CH = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(;N = 1, param = CH)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    #E.I = [param[5]*nA]
    E.I = [current_search(param)*nA]

    SNN.sim!([E], []; dt =0.001*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = ALLEN_DURATION+ALLEN_DELAY+500ms)
    #v = SNN.getrecord(E, :v)
    vec = SNN.vecplot(E, :v)

    py"""
    from neo import AnalogSignal
    #from neo import AnalogSignal
    import quantities as pq
    #import julia
    #julia.setup()
    vm2=[]
    vm2.append(0)
    """

    vec2=py"AnalogSignal("vec",units=pq.mV,sampling_rate=1.0*pq.Hz)"

    vec2 |> display
    vec

end

lower,upper = get_ranges()
function initd()
    population = initf(10)
    garray = zeros((length(population)[1], length(population[1])))
    for (i,p) in enumerate(population)
        garray[i,:] = p
    end
    garray[1,:]
end
cb = Evolutionary.ConstraintBounds(lower,upper,lower,upper)
selections = [:roulette=>rouletteinv, :sus=>susinv, :rank=>ranklinear(1.5)]
crossovers = [:discrete=>discrete, :intermediate0=>intermediate(0.), :intermediate0_25=>intermediate(0.5), :line=>line(0.2)]
mutations = [:domrng0_5=>domainrange(fill(0.5,4)), :uniform=>uniform(3.0), :gaussian=>gaussian(0.6)]
etas = [0.25,0.35,0.5,0.75]#,0.8]
#temp_GA
function twiny(sp::Plots.Subplot)
    sp[:top_margin] = max(sp[:top_margin], 30Plots.px)
    plot!(sp.plt, inset = (sp[:subplot_index], bbox(0,0,1,1)))
    twinsp = sp.plt.subplots[end]
    twinsp[:xaxis][:mirror] = true
    twinsp[:background_color_inside] = RGBA{Float64}(0,0,0,0)
    Plots.link_axes!(sp[:yaxis], twinsp[:yaxis])
    twinsp
end
#twiny(plt::Plots.Plot = current()) = twiny(plt[1])
#A
#CHOLECYSTITIS

function dontdo
    meta_param_dict = Dict()
    for (sn,ss) in selections, (xn,xovr) in crossovers, (mn,ms) in mutations, (ɛ) in etas
        xn == :discrete && (mn == :uniform || mn == :gaussian) && continue # bad combination
        xn == :line && mn == :gaussian && continue # bad combination
        temp_GA = GA(
            populationSize = 10,
            ɛ = ɛ,
            selection = ss,
            crossover = xovr,
            mutation = ms
            #ɛ = 0.125,
            #selection = ranklinear(1.5),
            #crossover = line(0.25),
            #mutation = gaussian(0.125)
        )
        println()
        println("GA:$(sn):$(xn):$(mn):ɛ=$ɛ) => F before crash...")
        println()

        result = Evolutionary.optimize(zz_, initd, temp_GA,
            Evolutionary.Options(iterations=12, successive_f_tol=25, show_trace=false, store_trace=false)
        )
        #print(result)
        fitness = minimum(result)
        meta_param_dict[:fitness] = Dict()
        meta_param_dict[:fitness][:"sn"] = sn
        meta_param_dict[:fitness][:"xn"] = xn
        meta_param_dict[:fitness][:"mn"] = mn
        meta_param_dict[:fitness][:"ɛ"] = ɛ

        println("GA:$(sn):$(xn):$(mn):ɛ=$ɛ) => F: $(minimum(result))")# C: $(Evolutionary.iterations(result))")

        extremum = Evolutionary.minimizer(result)
        meta_param_dict[:fitness][:"extremum"] = extremum
        save("meta_param_dict.jld", "meta_param_dict", meta_param_dict)#,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes)

        vec = checkmodel(extremum)
        #py"vvec" = vec
        #py"""
        #print(vecc)
        #print(type(vecc))
        #ass = AnalogSignal(vvec)
        #"""
        plot(vec)|> display

        plot(vmgtt[:],vmgtv[:])|> display
        #, label = "randData", ylabel = "Y axis",color = :red, legend = :topleft, grid = :off, xlabel = "Numbers Rand")
        #p = twiny()
        #plot!(p,vec, label = "log(x)", legend = :topright, box = :on, grid = :off, xlabel = "Log values") |> display

    end
end
#=
function plot_pop(pop)
    pop = filter(indiv -> indiv.rank <= 1, pop) #keeps only the non-dominated solutions
    scatter3d(map(x -> x.y[1], pop), map(x -> x.y[2], pop),  map(x -> x.y[3], pop), markersize = 1) |> display
    sleep(0.1)
end

#bc1 = BinaryCoding(2, [:Float64,:Float64], lower[1:2], upper[1:2])
bc = BinaryCoding(4, [:Float64,:Float64,:Float64,:Float64], lower, upper)
=#
#=
function CV(x)
    #sumW = dot(x, w)
    println(x[1])
    println("gets here")
    return x[1]# <= c ? 0 : sumW - c
end

function two_bits_flip!(bits)
    for i = 1:2
        n = rand(1:length(bits))
        bits[n] = 1 - bits[n]
    end
end
function one_point_crossover!(parent_a, parent_b, child_a, child_b)
    n = length(parent_a)
    cut = rand(1:n-1)

    child_a[1:cut] .= parent_a[1:cut]
    child_a[cut+1:n] .= parent_b[cut+1:n]

    child_b[1:cut] .= parent_b[1:cut]
    child_b[cut+1:n] .= parent_a[cut+1:n]
end
function two_point_crossover!(bits_a, bits_b, child1, child2)
    cut_a = cut_b = rand(1:length(bits_a)-1)
    while(cut_b == cut_a)
        cut_b = rand(1:length(bits_a))
    end
    cut_a,cut_b = minmax(cut_b,cut_a)

    copyto!(child1, 1, bits_a, 1, cut_a-1)
    copyto!(child1, cut_a, bits_b, cut_a, cut_b-cut_a+1)
    copyto!(child1, cut_b+1, bits_a, cut_b+1, length(bits_a)-cut_b)

    copyto!(child2, 1, bits_b, 1, cut_a-1)
    copyto!(child2, cut_a, bits_a, cut_a, cut_b-cut_a+1)
    copyto!(child2, cut_b+1, bits_b, cut_b+1, length(bits_a)-cut_b)
end
function crossover!(ind_a, ind_b, child_a, child_b)
    two_point_crossover!(ind_a.x, ind_b.x, child_a.x, child_b.x)
end

#default_mutation!(p::Vector{Int}) = rand_swap!(p)
#default_mutation!(b::T) where T<:AbstractVector{Bool} = rand_flip!(b)

function rand_flip!(bits)
    print(typeof(bits))
    nb = length(bits.x)
    for i = 1:nb
        if rand() < 1/nb
            @inbounds bits[i] = 1 - bits[i]
        end
    end
end
=#
#mutate!(ind::Indiv, fmut!) =
#rand_flip!(ind.x)

#nsga_max(popsize, nbgen, z, init, fCV = CV, fcross = one_point_crossover!)
#nsga_max(popsize, nbgen, z, init, fCV = CV, fplot = plot_pop, plotevery = 5)
#popsize = 10
#nbgen = 10
#seed = initf(10)
#pf = nsga_max(popsize, nbgen, z, bc, pmut = 0.85, plotevery = 5,seed=seed, fcross = crossover!)#, fCV = ConstraintVioloation)#fmut = rand_flip!)#, fplot = plot_pop, plotevery = 1)#, fCV = CV)
#model = NSGAII.decode!(pf[1].x,bc)
#println(model.p)
#checkmodel(model.p)
#init = () -> NSGAII.BinaryCodedIndiv(bitrand(bc.nbbitstotal), zeros(bc.nbvar))
#ind = init()
#re_model = NSGAII.decode!(ind.x,bc)

#indiv -> (decode!(indiv, bc) ; z(indiv.p))
#println(indiv.p)
#_z = indiv -> (decode!(indiv, bc) ; z(indiv.p))

#param = getproperty.(result[1], :x)
#NSGAII.decode!(bincoded, bc)
#println(res)
#=
options = Evolutionary.Options(iterations=10, abstol=1e-5)

#bounds = ConstraintBounds(lower,upper,[],[])

method = GA(populationSize=10, crossoverRate=0.8, mutationRate=0.05, selection=susinv)#, crossover=MILX(), mutation=MIPM(lx,ux))
Random.seed!(10);
#population = Evolutionary.initial_population(method, model.p)

# These could be the same if we could restrict g below not to be an AbstractArray
function NonDifferentiable1(f, x::AbstractArray, F::Real = real(zero(eltype(x))); inplace = true)
    xnans = x_of_nans(x)
    NonDifferentiable1{typeof(F),typeof(xnans)}(f, F, xnans, [0,])
end
function NonDifferentiable1(f, x::AbstractArray, F::AbstractArray; inplace = true)
    f = !inplace && (F isa AbstractArray) ? f!_from_f(f, F, inplace) : f
    xnans = x_of_nans(x)
    NonDifferentiable1{typeof(F),typeof(xnans)}(f, F, xnans, [0,])
end

# this is the g referred to above!
NonDifferentiable1(f, g,        x::AbstractArray, F::Union{AbstractArray, Real} = real(zero(eltype(x)))) = NonDifferentiable1(f, x, F)
NonDifferentiable1(f, g, h,     x::TX, F) where TX  = NonDifferentiable(f, x, F)
NonDifferentiable1(f, g, fg, h, x::TX, F) where TX  = NonDifferentiable(f, x, F)

# Objective function
function rastrigin(x::AbstractVector{T}) where {T <: AbstractFloat}
    n = length(x)
    return 10n + sum([ x[i]^2 - 10cos(convert(T,2π*x[i])) for i in 1:n ])
end

# Parameters
N = 3
P = 100
initState = ()->rand(N)
println(initState)
println("initstate")
println()
=#
#=
#initState = garray[1]
# Testing: (μ/μ_I,λ)-σ-Self-Adaptation-ES
# with non-isotropic mutation operator y' := y + (σ_1 N_1(0, 1), ..., σ_N N_N(0, 1))
m = ES(mu = 15, lambda = P)
m.μ == 15
m.λ == P
result = Evolutionary.optimize( rastrigin,
    initState,
    ES(
        initStrategy = AnisotropicStrategy(N),
        recombination = average, srecombination = average,
        mutation = gaussian, smutation = gaussian,
        selection=:comma,
        μ = 15, λ = P
    ),Evolutionary.Options(iterations=1000, show_trace=false)
)
println("(15,$(P))-σ-SA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
#test_result(result, N, 1e-1)
print(result)
# Testing: CMA-ES
result = Evolutionary.optimize(rastrigin, initState, CMAES(mu = 15, lambda = P))
println("(15/15,$(P))-CMA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
#test_result(result, N, 1e-1)

selections = [:roulette=>rouletteinv, :sus=>susinv, :rank=>ranklinear(1.5)]
crossovers = [:discrete=>discrete, :intermediate0=>intermediate(0.), :intermediate0_25=>intermediate(0.5), :line=>line(0.2)]
mutations = [:domrng0_5=>domainrange(fill(0.5,N)), :uniform=>uniform(3.0), :gaussian=>gaussian(0.6)]

for (sn,ss) in selections, (xn,xovr) in crossovers, (mn,ms) in mutations
    xn == :discrete && (mn == :uniform || mn == :gaussian) && continue # bad combination
    xn == :line && mn == :gaussian && continue # bad combination
    result = Evolutionary.optimize( rastrigin, initState,
        GA(
            populationSize = P,
            ɛ = 0.1,
            selection = ss,
            crossover = xovr,
            mutation = ms
        ), Evolutionary.Options(iterations=1000, successive_f_tol=25)
    )
    println(result)
    println("GA:$(sn):$(xn):$(mn)(N=$(N),P=$(P),x=.8,μ=.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
end
# Testing: GA
#m = GA(epsilon = 10.0)
#@test m.ɛ == 10

selections = [:roulette=>rouletteinv, :sus=>susinv, :rank=>ranklinear(1.5)]
crossovers = [:discrete=>discrete, :intermediate0=>intermediate(0.), :intermediate0_25=>intermediate(0.5), :line=>line(0.2)]
mutations = [:domrng0_5=>domainrange(fill(0.5,N)), :uniform=>uniform(3.0), :gaussian=>gaussian(0.6)]

for (sn,ss) in selections, (xn,xovr) in crossovers, (mn,ms) in mutations
    xn == :discrete && (mn == :uniform || mn == :gaussian) && continue # bad combination
    xn == :line && mn == :gaussian && continue # bad combination
    result = Evolutionary.optimize( rastrigin, initState,
        GA(
            populationSize = P,
            ɛ = 0.1,
            selection = ss,
            crossover = xovr,
            mutation = ms
        ), Evolutionary.Options(iterations=1000, successive_f_tol=25)
    )
    println("GA:$(sn):$(xn):$(mn)(N=$(N),P=$(P),x=.8,μ=.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
end
=#

#=

function initd()
    population = initf(1)
    population
end
=#

#println(garray)
#aa = zero(eltype(first(garray)))
 #println(aa)
#nd = NonDifferentiable1(zz,first(garray))
#println(nd)
#(method::GA, options, objfun, population)
#GA_ = Evolutionary.GA()
#initState = ()->rand(N)
#N = size(seed[1])[1]
#constraints(z) = [sum(z[1:end-1]),z[2]]
#state = Evolutionary.initial_state(GA_, options, z, garray[1])

#Evolutionary.update_state!(z, constraints, state, garray[1], GA, 1)

#selections = [:roulette=>rouletteinv, :sus=>susinv, :rank=>ranklinear(1.5)]
#crossovers = [:discrete=>discrete, :intermediate0=>intermediate(0.), :intermediate0_25=>intermediate(0.5), :line=>line(0.2)]
#mutations = [:domrng0_5=>domainrange(fill(0.5,N)), :uniform=>uniform(3.0), :gaussian=>gaussian(0.6)]

#res2 = Evolutionary.optimize(z, bounds, method, options)

#result = Evolutionary.optimize(z, lower, upper, mthd, opts)
#function optimize(f, lower, upper, method::M,
#                  options::Options = Options(;default_options(method)...)
#                 ) where {M<:AbstractOptimizer}
#    bounds = ConstraintBounds(lower,upper,[],[])
#    optimize(f, bounds, method, options)
#end

@test minimum(result) ≈ 2.0 atol=1e-1
@test Evolutionary.minimizer(result) ≈ [0.5, 1] atol=1e-1

f1(x) = 2*x[1]+x[2]

#lx = [0.0, 0]
#ux = [1.6, 1]

#cons(x) = [1.25-x[1]^2-x[2], x[1]+x[2]]
#lc = [-Inf, -Inf]
#uc = [0.0, 1.6]

tc = [Float64, Float64, Float64, Float64, Float64]
c = MixedTypePenaltyConstraints(PenaltyConstraints(1e3, cb, z), tc)
#c = MixedTypePenaltyConstraints(WorstFitnessConstraints(cb, cons), tc)
#init = ()->Real[rand(Float64), rand(0:1)]

#-----------------------------------------------------------------

f8(z) = sum((z[4:end-1].-1).^2) - log(z[end]+1) + sum((z[i]-i)^2 for i in 1:3)

@test f8(Real[0.2, 1.280624, 1.954483, 1, 0, 0, 1]) ≈ 3.55746 atol=1e-6

cons(z) = [
    sum(z[1:end-1]),
    sum(v->v^2, z[[1,2,3,6]]),
    z[1]+z[4],
    z[2]+z[5],
    z[3]+z[6],
    z[1]+z[7],
    z[2]^2+z[5]^2,
    z[3]^2+z[6]^2,
    z[2]^2+z[6]^2,
]

lx = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ux = [Inf, Inf, Inf, 1.0, 1.0, 1.0, 1.0]

lc = Float64[]
uc = [5.0, 5.5, 1.2, 1.8, 2.5, 1.2, 1.64, 4.25, 4.64]

tc = [Float64, Float64, Float64, Int, Int, Int, Int]

cb = Evolutionary.ConstraintBounds(lx,ux,lc,uc)
# c = MixedTypePenaltyConstraints(PenaltyConstraints(1E, cb, cons), tc)
c = MixedTypePenaltyConstraints(WorstFitnessConstraints(cb, cons), tc)
init = ()->Real[rand(Float64,3); rand(0:1,4)]

opts = Evolutionary.Options(iterations=2000)
mthd = GA(populationSize=300, ɛ=0.05, crossoverRate=0.8, mutationRate=0.01, selection=susinv, crossover=MILX(0.0,0.5,0.3), mutation=MIPM(lx,ux))

Random.seed!(544);
result = Evolutionary.optimize(f8, c, init, mthd, opts)
@test minimum(result) < 4

#=
opts = Evolutionary.Options(iterations=100, show_every=1, show_trace=true, store_trace=true)
#opts = Evolutionary.Options(iterations=500, abstol=1e-5)
#mthd = GA(populationSize=10, crossoverRate=0.8, mutationRate=0.05, selection=susinv, crossover=MILX(), mutation=MIPM(lx,ux))
mthd = GA(populationSize = 120, ɛ = 0.03, crossoverRate=0.8, mutationRate=0.01, selection=rouletteinv)#, crossover=LX(0.0,4.0), mutation = PM(lx,ux,1.0))

tc = [Float32]#, Float64, Float64, Int, Int, Int, Int]

cb = Evolutionary.ConstraintBounds(lx,ux,lc,uc)
# c = MixedTypePenaltyConstraints(PenaltyConstraints(1e2, cb, cons), tc)
#c = MixedTypePenaltyConstraints(PenaltyConstraints(cb, cons), tc)

Random.seed!(10);
result = Evolutionary.optimize(evaluate, population[1], initf, mthd, opts)
#res = Evolutionary.optimize(fitness, model, algo, opts)

println()
println()
println()
println()

println("error $error")
println()
println()

function init_b()
    for i in 0:9
        gene = []
        for (l,u) in zip(lower,upper)
            p1 = rand(l:u, 1)
            append!(gene,p1)
        end
        println(size(gene))
        append!(genes,gene)
        bincoded = NSGAII.encode(gene, bin_code)
        append!(genesb,bincoded)
    end
    #bincoded = NSGAII.encode(gene_out, bc)
    #return bincoded

    genesb
end
#println(size(genes))
#println(size(genesb))

#bincoded = NSGAII.encode(seed, bin_code)
#=
export transdict
function transdict(x)
    bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

    #py"""
    #def map_dict(x):
    #    genes_out = x
    #    genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
    #    return genes_out, genes_out_dic
    #"""
    decoded = x

    try
       decoded = NSGAII.decode(x,bc)
    catch
       decoded = x
    end
    _,god_dict = py"map_dict"(decoded)
    #bincoded = NSGAII.encode(genes_out, bc)
    return god_dict#, genes_out_dic
end
=#
##py"""
#lower_list = [v[0] for k,v in ranges.items()]
#upper_list = [v[1] for k,v in ranges.items()]
##"""

#export init_function
function init_function()
    #const
    bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

    py"""
    def pheno_map():

        lower_list = [v[0] for k,v in ranges.items()]
        upper_list = [v[1] for k,v in ranges.items()]
        gene_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
        #gene_out_dic = OrderedDict({k:gene_out[i] for i,k in enumerate(ranges.keys()) })
        return gene_out#, gene_out_dic
    """
    gene_out = py"pheno_map"()
    bincoded = NSGAII.encode(gene_out, bc)
    return bincoded
end

function z(x)
    float_coded = NSGAII.decode(x, bc)

    #god = transdict(x)
    #contents = py"z_py"(py"evaluate",god)
    #return contents
end

#=
function references_for_H()
    bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

    py"""
    from sklearn.model_selection import ParameterGrid
    grid = ParameterGrid(ranges)
    genes = []
    for g in grid:
        genes.append(list(g.values()))
    genes_out_dic = OrderedDict({k:genes[i] for i,k in enumerate(ranges.keys()) })

    """
    genes_out = py"genes"
    bcd = []
    for g in genes_out
         bincoded = NSGAII.encode(g, bc)
         append!(bcd,bincoded)
    end
    return bcd
end
=#

export z
function z(x)
    god = transdict(x)
    contents = py"z_py"(py"evaluate",god)
    return contents
end

#bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

module HHNSGA
    #import Pkg;

    include("../src/SpikingNeuralNetworks.jl")
    include("../src/units.jl")
    include("../src/plot.jl")

    SNN = SpikingNeuralNetworks.SNN
    using Random
    using PyCall
    py"""
    import copy
    #import matplotlib
    #import matplotlib.pyplot as plt
    import numpy as np
    #matplotlib.get_backend()
    import sys
    import os
    sys.path.append(os.getcwd())
    """
    #include("../src/SpikingNeuralNetworks.jl")
    #SNN = SpikingNeuralNetworks.SNN")
    #import pickle
    #import numpy as np
    #import random
    #from sciunit.scores.collections import ScoreArray
    #from sciunit import TestSuite
    #import sciunit
    #import random
    #from neuronunit.tests.fi import RheobaseTestP
    #from collections import OrderedDict
    """
    #from neuronunit import get_neab
    #import pdb
    #pdb.set_trace()
    #cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
    for test in cell_tests.values():
        if "Rheobase test" in test.keys():
            temp_test = {k:v for k,v in test.items()}
            break
    rt = temp_test["Rheobase test"]
    #rtp = RheobaseTestP(rt.observation)
    JUIZI = {
        'a': 0.02,
        'b': 0.2,
        'c': -65,
        'd': 8,
    }

    #ranges = OrderedDict({k:[v-0.5*np.abs(v),v+0.5*np.abs(v)] for k,v in copy.copy(JUIZI).items()})

    ranges = OrderedDict(temp)

    N = len(JUIZI)
    """
    #N = py"N"
    #ranges = py"ranges"
    #H1=[values(ranges)]
    #current_params = py"rt.params"
    #simple.attrs = py"JUIZI"
    py"""
    from izhi import IZModel

    def evaluate(test,god):
        model = IZModel(attrs=god)
        model.attrs = god
        rt = test['Rheobase test']

        rheo = rt.generate_prediction(model)

        scores_ = []
        test = {k:v for k,v in test.items() if v.observation is not None}
        for temp_test in test.values():
            if 'injected_square_current' not in temp_test.params.keys():
               temp_test.params['injected_square_current'] = {}
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "InjectedCurrentAPWidthTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "Injected CurrentAPThresholdTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "InjectedCurrentAPAmplitudeTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

            if temp_test.name == "RheobaseTest":
               temp_test.params['injected_square_current']['amplitude'] = rheo['value']
               temp_test.params['amplitude'] = rheo['value']

        to_map = [(tt,model) for tt in test.values() ]
        def map_score(content):
            temp_test = content[0]
            model = content[1]
            temp_test.score_type = sciunit.scores.ZScore
            try:
                score = temp_test.judge(model)
                score = np.abs(score.log_norm_score)
            except:
                try:
                    score = np.abs(float(score.raw_score))
                except:
                    score = 100
            if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
                score = 100
            return score

        scores_ = list(map(map_score,to_map))
        tt = TestSuite(list(test.values()))
        SA = ScoreArray(tt, scores_)
        errors = tuple(SA.values,)
        return errors
    """

    py"""
    def z_py(evaluate,god):

        cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
        tests = cell_tests['Neocortex pyramidal cell layer 5-6']

        model_type="RAW"
        #from neuronunit.optimisation import make_sim_tests
        fps = ['a','b','c','d']
        #sim_tests, OM, target = make_sim_tests.test_all_objective_test(fps,model_type=model_type)
        SA = evaluate(tests,god)
        return SA
    """

    export transdict
    function transdict(x)
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        #py"""
        #def map_dict(x):
        #    genes_out = x
        #    genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
        #    return genes_out, genes_out_dic
        #"""
        decoded = x

        try
           decoded = NSGAII.decode(x,bc)
        catch
           decoded = x
        end
        _,god_dict = py"map_dict"(decoded)
        #bincoded = NSGAII.encode(genes_out, bc)
        return god_dict#, genes_out_dic
    end
    ##py"""
    #lower_list = [v[0] for k,v in ranges.items()]
    #upper_list = [v[1] for k,v in ranges.items()]
    ##"""

    export init_function
    function init_function()
        #const
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        py"""
        def pheno_map():

            lower_list = [v[0] for k,v in ranges.items()]
            upper_list = [v[1] for k,v in ranges.items()]
            gene_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
            #gene_out_dic = OrderedDict({k:gene_out[i] for i,k in enumerate(ranges.keys()) })
            return gene_out#, gene_out_dic
        """
        gene_out = py"pheno_map"()
        bincoded = NSGAII.encode(gene_out, bc)
        return bincoded
    end

    function references_for_H()
        bc = BinaryCoding(4, [:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")

        py"""
        from sklearn.model_selection import ParameterGrid
        grid = ParameterGrid(ranges)
        genes = []
        for g in grid:
            genes.append(list(g.values()))
        genes_out_dic = OrderedDict({k:genes[i] for i,k in enumerate(ranges.keys()) })

        """
        genes_out = py"genes"
        bcd = []
        for g in genes_out
             bincoded = NSGAII.encode(g, bc)
             append!(bcd,bincoded)
        end
        return bcd
    end

    export z
    function z(x)
        god = transdict(x)
        contents = py"z_py"(py"evaluate",god)
        return contents
    end

end
#create_sysimage(:HHNSGA, sysimage_path="sys_HHNSGA.so", precompile_execution_file="precompile_HHNSGA.jl")
=#


text:
"""
    NumFactor

Output type from the numfactor function. The fields of this type provide the
number of factors for each information criterion, as well as the underlying information
criterion vectors for each criterion. A description of the information criterion
can be found in Bai, Ng (2002) section 5. \n

The estimated number of factors can be obtained using numfactor(nf::NumFactor)
to get a vector of all 6 estimates, corresponding to criterions pc1, pc2, pc3,
ic1, ic2, ic3. Alternatively, you can get just one of these criterion using
numfactor(nf::NumFactor, y::Symbol) where y can take values :pc1, :pc2, :pc3,
:ic1, :ic2, :ic3. y::Vector{Symbol} is also accepted as input. \n

The average estimate of number of factors taken across all 6 criterion, and rounded
to the nearest integer, can be obtained using avgnumfactor(nf::NumFactor) \n
"""
struct NumFactor
    pc1numfac::Int ; pc2numfac::Int ; pc3numfac::Int ; ic1numfac::Int ; ic2numfac::Int ; ic3numfac::Int
    pc1vec::Vector{Float64} ; pc2vec::Vector{Float64} ; pc3vec::Vector{Float64}
    ic1vec::Vector{Float64} ; ic2vec::Vector{Float64} ; ic3vec::Vector{Float64}
end
NumFactor() = NumFactor(-1, -1, -1, -1, -1, -1)
numfactor(nf::NumFactor)::Vector{Int} = [nf.pc1numfac, nf.pc2numfac, nf.pc3numfac, nf.ic1numfac, nf.ic2numfac, nf.ic3numfac]
function numfactor(nf::NumFactor, icsym::Symbol)::Int
    icsym == :pc1 && return nf.pc1numfac
    icsym == :pc2 && return nf.pc2numfac
    icsym == :pc3 && return nf.pc3numfac
    icsym == :ic1 && return nf.ic1numfac
    icsym == :ic2 && return nf.ic2numfac
    icsym == :ic3 && return nf.ic3numfac
    error("Invalid information criterion symbol $(icsym). Please use :pc1, :pc2, :pc3, :ic1, :ic2, or :ic3")
end
numfactor(nf::NumFactor, icsymvec::Vector{Symbol})::Vector{Int} = [ numfactor(nf, icsym) for icsym in icsymvec ]
avgnumfactor(nf::NumFactor)::Int = Int(round(mean(numfactor(nf))))
function Base.show(io::IO, x::NumFactor)
    println(io, "Estimated number of factors:")
    println(io, "    pcp1 = $(x.pc1numfac)")
    println(io, "    pcp2 = $(x.pc2numfac)")
    println(io, "    pcp3 = $(x.pc3numfac)")
    println(io, "    icp1 = $(x.ic1numfac)")
    println(io, "    icp2 = $(x.ic2numfac)")
    println(io, "    icp3 = $(x.ic3numfac)")
end

"""
    numfactor(x::Matrix{<:Number} ; kwargs...)::NumFactor \n

Estimate the number of common factors in dataset x, using the method proposed
in Bai, Ng (2002) Determining the Number of Factors in Approximate Factor Models [Econometrica]. \n

The rows in x should correspond to observations, and columns correspond to variables.
See ?NumFactor for more information on the output type. \n

This function accepts the following keyword arguments: \n
    - kmax::Int=10 <- The maximum number of common factors to test up to. Note that
    the pcp criterion tend to over-estimate the number of common factors if kmax is
    set too large, so it is not advised to set this to an arbitrarily large value. \n
    - covmatdim::Symbol=:auto <- If x is a TxJ matrix, then set to :cols to use
    a JxJ covariance matrix, :rows for a TxT covariance matrix, or :auto to use
    whichever is smaller of J or T. From a statistical perspective, the choice
    is irrelevant, so :auto is recommended. \n
    - xcent::Bool=false <- Set to true if the data in x is already centred. \n
    - standardize::Bool=true <-- Set to true to standardize x to have unit variance \n
"""
function numfactor(x::Matrix{<:Number} ; kmax::Int=10, covmatdim::Symbol=:auto,
                   xcent::Bool=false, standardize::Bool=false)::NumFactor
    #Preliminaries
    (T, J) = (size(x, 1), size(x, 2))
    if covmatdim == :auto ; min(T, J) < kmax && return NumFactor()
    elseif covmatdim == :rows ; J < kmax && return NumFactor()
    elseif covmatdim == :cols ; T < kmax && return NumFactor()
    else ; error("Invalid covmatdim: $(covmatdim)")
    end
    (T < 2 || J < 2) && return NumFactor()
    !xcent && (x = x .- mean(x, 1))
    standardize && (x = x .* (1 ./ std(x, 1)))
    !any(covmatdim .== [:auto, :rows, :cols]) && error("covmatdim must be set to :auto, :cols, or :rows. Current value is invalid: $(covmatdim)")
    if covmatdim == :auto #Choose whichever dimension is smaller for covariance matrix
        T > J ? (covmatdim = :cols) : (covmatdim = :rows)
    end
    #Get covariance matrix
    if covmatdim == :cols ; xcov = x'*x
    elseif covmatdim == :rows ; xcov = x*x'
    else ; error("covmatdim must be set to :auto, :cols, or :rows. Current value is invalid: $(covmatdim)")
    end
    #Get eigenvalues and eigenvectors for largest kmax eigenvalues
    ef = eigfact(Symmetric(xcov), size(xcov,1)-kmax+1:size(xcov,1))
    (eigvalvec, eigvecmat) = (flipdim(ef.values, 1), flipdim(ef.vectors, 2))
    #Get estimated factors and factor loadings (see Bai, Ng (2002) section 3)
    if covmatdim == :cols
        facload = sqrt(J) * eigvecmat
        estfac = (1/J) * x * facload
    elseif covmatdim == :rows
        estfac = sqrt(T) * eigvecmat
        facload = ((1/T) * x' * estfac)
    else ; error("Logic fail. It should have been impossible to reach this point. Please file an issue.")
    end
    #Get average sum of squared residuals (see Bai, Ng (2002) equation 7)
    v = [ (1/(T*J)) * sum((x' - facload[:, 1:k] * estfac[:, 1:k]').^2) for k = 1:kmax ]
    #Set some useful values
    tj1 = (T + J) / (T * J)
    tj2 = 1 / tj1
    c2 = (min(sqrt(T), sqrt(J)))^2
    #Get PC criterion vectors
    pc1vec = v + (v[end]*tj1*log(tj2))*collect(1:kmax)
    pc2vec = v + (v[end]*tj1*log(c2))*collect(1:kmax)
    pc3vec = v + (v[end]*(log(c2)/c2))*collect(1:kmax)
    #Get IC criterion vectors
    ic1vec = log.(v) + (tj1*log(tj2))*collect(1:kmax)
    ic2vec = log.(v) + (tj1*log(c2))*collect(1:kmax)
    ic3vec = log.(v) + (log(c2)/c2)*collect(1:kmax)
    #Return a NumFactor object
    (pc1nf, pc2nf, pc3nf, ic1nf, ic2nf, ic3nf) = (indmin(pc1vec), indmin(pc2vec), indmin(pc3vec), indmin(ic1vec), indmin(ic2vec), indmin(ic3vec))
    return NumFactor(pc1nf, pc2nf, pc3nf, ic1nf, ic2nf, ic3nf, pc1vec, pc2vec, pc3vec, ic1vec, ic2vec, ic3vec)
end

ext: jl


text:
module StorageServices

using ..Azure
using Azure.StorageManagementClient
using Azure.REST
using Azure.REST: AccountKey, service_response, parse_service_response, issuccess
using URIs
using MbedTLS
using Dates
using Base64
using Mmap

include("common.jl")
include("blob.jl")
include("file.jl")
include("sas.jl")

end


text:
#   This unit allows to install the PosDefManifold package
#   and to build the documentation locally.
#   v 0.1.3 - last update 10th of Mai 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home
#
#   DIRECTIONS:
#   1) If you have installed the PosDefPackage from github, uninstall it.
#   2) Change the `juliaCodeDir` path here below to the path
#           where the PosDefMaifold folder is located on your computer.
#   3) Put the cursor in this unit and hit SHIFT+CTRL+ENTER

begin
    projectName="PosDefManifold"
    juliaCodeDir=homedir()*"\\Documents\\Code\\julia\\"
    # this file will be copied once to become the logo
    logoPngFile=homedir()*"\\Pictures\\cenerentola.png"

    projectDir       =   juliaCodeDir*projectName*"\\"
    srcDir           =   projectDir*"src\\"
    testDir          =   projectDir*"test\\"
    MainModule       =   srcDir*projectName*".jl"
    docsDir          =   projectDir*"docs\\"
    docsSrcDir       =   docsDir*"src\\"
    docsSrcAssetsDir =   docsSrcDir*"assets\\"
    docsBuildDir     =   docsDir*"build\\"

    # creates directories
    ispath(srcDir)           || mkpath(srcDir)
    ispath(testDir)          || mkpath(testDir)
    ispath(docsBuildDir)     || mkpath(docsBuildDir)
    ispath(docsSrcAssetsDir) || mkpath(docsSrcAssetsDir)

    # Copy logo png file if it does not exist
    path, ext=splitext(logoPngFile)
    destination=docsSrcAssetsDir*"logo.png"
    if isfile(logoPngFile) && !isfile(destination)
        if ext==".png" cp(logoPngFile, docsSrcAssetsDir*"logo.png")
        else @warn("The logo file is not a .png file.")
        end
    end

    # just an alias for convenience
    wd=pwd()

    #--------------------------------
    # for building the documentation:
    #--------------------------------

    push!(LOAD_PATH, srcDir)

    using Documenter, DocumenterTools, Statistics,
    Base.Threads, LinearAlgebra, BenchmarkTools, Revise,
    PosDefManifold

    cd(docsDir)

    #chmod(projectDocDir, recursive=true)
    #chmod(projectDocBuildDir, recursive=true)
    if  Base.HOME_PROJECT[] !== nothing
        Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
    end

    clipboard("""makedocs(sitename="PosDefManifold", modules=[PosDefManifold])""")
    @info("\nhit CTRL+V+ENTER on the REPL for building the documentation.");
end

max_stars_repo_path: build/START HERE.jl

max_issues_repo_path: build/START HERE.jl

max_forks_repo_path: build/START HERE.jl


text: Dict{String, Any}("neg" => Dict{String, Any}("value" => "0.03", "type" => "float"), "upper" => Dict{String, Any}("value" => "300.0", "type" => "float"), "zero" => Dict{String, Any}("value" => "3.0", "type" => "float"), "pointupper" => Dict{String, Any}("value" => "310.0", "type" => "float"), "lower" => Dict{String, Any}("value" => "300.0", "type" => "float"), "minustenth" => Dict{String, Any}("value" => "-0.1", "type" => "float"), "pointlower" => Dict{String, Any}("value" => "310.0", "type" => "float"), "pos" => Dict{String, Any}("value" => "300.0", "type" => "float"))


text:
type Callback
    callback::Function
    name::String

    function Callback(callback,name)
        cb = new(callback,name)
        callbacks[name] = cb
        finalizer(cb, cb -> delete!(callbacks, cb.name))
        cb
    end
end
const callbacks = Dict{String,Callback}() # name => args -> f(args...)

# Empty callback to notify the Server that a new page is loaded and its WebSocket is ready.
Callback("connected") do
end

# Callback used to cleanup when the browser navigates away from the page.
Callback("unloaded") do session_id
    session = sessions[session_id]
    delete!(pages[session.route].sessions,session_id)
    delete!(sessions,session_id)
end

# Callback used to message a specified connected browser.
Callback("message") do args
    message(args...)
end

# Callback used to broadcast to all connected browsers.
Callback("broadcast") do args
    broadcast(args...)
end

# Callback used for blocking Julia control flow until notified by the WebSocket.
Callback("notify") do name
    if haskey(conditions,name)
        notify(conditions[name])
    else
        error("""Condition "$name" was not found.""")
    end
end


text:
export @use

const VALID_KWARGS = (:density, :gravity, :heatcapacity, :conductivity, :latent_heat, :radioactive_heat, :shearheating,
        :pwave_velocity, :swave_velocity, :meltfraction)

macro use(args...)
        checkargs(args...)
        esc(create_module(__module__, args...))
        esc(use_module(__module__, args[1]))
end

function checkargs(args...)
        if (length(args) < 2)
                error("arguments missing.")
        end # NOTE: this needs to be adapted to: 1 + nb math symbols needed
        if (args[1] != :GeoParamsAliases)
                error("the first argument must be `GeoParamsAliases`")
        end
        if !all(is_kwarg.(args[2:end]))
                error("all arguments starting from the second must be keyword arguments (see macro documentation).")
        end
end

is_kwarg(arg) = isa(arg, Expr) && (arg.head == :(=))
split_kwargs(kwargs) = Dict(x.args[1] => x.args[2] for x in kwargs)

function validate_kwargkeys(kwargs::Dict, macroname::String)
        for k in keys(kwargs)
                if !(k in VALID_KWARGS)
                        error("Invalid keyword argument in $macroname call: `$k`. Valid keyword arguments are: `$(join(valid_kwargs, "`, `"))`.")
                end
        end
end

function create_module(caller::Module, modulename::Symbol, kwargs_expr::Expr...)
        if isdefined(caller, modulename)
                error("Module $modulename already exists in $caller.")
        end
        kwargs = split_kwargs(kwargs_expr)
        validate_kwargkeys(kwargs, "@use GeoParamsAliases")
        exports = Symbol[]
        aliases = Expr[]
        for key in keys(kwargs)
                param = kwargs[key]
                param! = Symbol(param, "!")
                fname = Symbol("compute_" * String(key))
                fname! = Symbol("compute_" * String(key) * "!")

                # Exports
                push!(exports, param)
                # Aliases
                push!(aliases, quote $param(args...) = GeoParams.$fname(args...) end)
                # Only in case compute_...! method exists set alias for it as well
                if isdefined(GeoParams, fname!)
                        push!(exports, param!)
                        push!(aliases, quote $param!(args...) = GeoParams.$fname!(args...) end)
                end
        end

        module_expr = :(module $modulename # NOTE: there cannot be any newline before 'module $modulename' or it will create a begin end block and the module creation will fail.
        import GeoParams
        export $(exports...)
        $(aliases...)
        end)
        @eval(caller, $module_expr)
end

use_module(caller::Module, modulename::Symbol) = @eval(caller, using $(Symbol(caller)).$modulename)


text:
using Metis
using Random
using Test
using SparseArrays
import LightGraphs

@testset "Metis.permutation" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    perm, iperm = Metis.permutation(S)
    @test isperm(perm)
    @test isperm(iperm)
    @test S == S[perm,perm][iperm,iperm]
end

@testset "Metis.partition" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T), alg in (:RECURSIVE, :KWAY), nparts in (3, 4)
        partition = Metis.partition(G, nparts, alg = alg)
        @test extrema(partition) == (1, nparts)
        @test all(x -> findfirst(==(x), partition) !== nothing, 1:nparts)
    end
end

@testset "Metis.separator" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T)
        parts = Metis.separator(G)
        @test extrema(parts) == (1, 3)
    end
end


text:
power = 2 ^ 7
modulo = 101 % 2

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
module Conversions

export convert_to_julia, convert_to_jni, convert_to_string, convert_to_vector

using JavaCall.JNI
using JavaCall.CodeGeneration

# Conversion of primitive types.

convert_to_julia(::Type{Bool}, x::jboolean)::Bool = (x == JNI_TRUE)
convert_to_jni(::Type{jboolean}, x::Bool)::jboolean = x ? JNI_TRUE : JNI_FALSE

convert_to_julia(::Type{Int8}, x::jbyte)::Int8 = x
convert_to_jni(::Type{jbyte}, x::Int8)::jbyte = x

convert_to_julia(::Type{Char}, x::jchar)::Char = Char(x)
# May overflow as julia chars have 4 bytes and java
# chars have 2 bytes. Should not happen for UTF-16 chars.
convert_to_jni(::Type{jchar}, x::Char)::jchar = jchar(x)

convert_to_julia(::Type{Int16}, x::jshort)::Int16 = x
convert_to_jni(::Type{jshort}, x::Int16)::jshort = x

convert_to_julia(::Type{Int32}, x::jint)::Int32 = x
convert_to_jni(::Type{jint}, x::Int32)::jint = x

convert_to_julia(::Type{Int64}, x::jlong)::Int64 = x
convert_to_jni(::Type{jlong}, x::Int64)::jlong = x

convert_to_julia(::Type{Float32}, x::jfloat)::Float32 = x
convert_to_jni(::Type{jfloat}, x::Float32)::jfloat = x

convert_to_julia(::Type{Float64}, x::jdouble)::Float64 = x
convert_to_jni(::Type{jdouble}, x::Float64)::jdouble = x

convert_to_julia(::Type{Nothing}, ::jvoid)::Nothing = nothing
convert_to_jni(::Type{jvoid}, ::Nothing)::jvoid = nothing

# Conversion of vectors for primitive types
for (julia_type, jni_type, jni_arraytype, jni_function) in [
    (:Bool, :jboolean, :jbooleanArray, :get_boolean_array_elements),
    (:Int8, :jbyte, :jbyteArray, :get_byte_array_elements),
    (:Char, :jchar, :jcharArray, :get_char_array_elements),
    (:Int16, :jshort, :jshortArray, :get_short_array_elements),
    (:Int32, :jint, :jintArray, :get_int_array_elements),
    (:Int64, :jlong, :jlongArray, :get_long_array_elements),
    (:Float32, :jfloat, :jfloatArray, :get_float_array_elements),
    (:Float64, :jdouble, :jdoubleArray, :get_double_array_elements)
]
    params = [
        :(::Type{Vector{$julia_type}}),
        :(array::$jni_arraytype)
    ]
    body = quote
        array_len = JNI.get_array_length(array)
        iscopy = Ref(JNI_TRUE)
        iscopyptr = convert(Ptr{jboolean}, pointer_from_objref(iscopy))
        returnptr = JNI.$jni_function(array, iscopyptr)
        map(x -> convert_to_julia($julia_type, x), unsafe_wrap(Vector{$jni_type}, returnptr, (array_len,)))
    end
    eval(generatemethod(
        :convert_to_julia,
        params,
        body
    ))
end

for (julia_type, jni_type, jni_arraytype, jni_create_array_fn, jni_set_array_region_fn) in [
    (:Bool, :jboolean, :jbooleanArray, :new_boolean_array, :set_boolean_array_region),
    (:Int8, :jbyte, :jbyteArray, :new_byte_array, :set_byte_array_region),
    (:Char, :jchar, :jcharArray, :new_char_array, :set_char_array_region),
    (:Int16, :jshort, :jshortArray, :new_short_array, :set_short_array_region),
    (:Int32, :jint, :jintArray, :new_int_array, :set_int_array_region),
    (:Int64, :jlong, :jlongArray, :new_long_array, :set_long_array_region),
    (:Float32, :jfloat, :jfloatArray, :new_float_array, :set_float_array_region),
    (:Float64, :jdouble, :jdoubleArray, :new_double_array, :set_double_array_region)
]
    params = [
        :(::Type{$jni_arraytype}),
        :(vec::Vector{$julia_type})
    ]
    body = quote
        vec_len = length(vec)
        array = JNI.$jni_create_array_fn(vec_len)
        jni_types_vec = map(x -> convert_to_jni($jni_type, x), vec)
        JNI.$jni_set_array_region_fn(array, 0, vec_len, jni_types_vec)
        array
    end
    eval(generatemethod(
        :convert_to_jni,
        params,
        body
    ))
end

# Auxiliary conversions for simplifying types interface
# between languages

function convert_to_string(::Type{String}, str::jstring)
    iscopy = Ref(JNI_TRUE)
    iscopyptr = convert(Ptr{jboolean}, pointer_from_objref(iscopy))
    unsafe_string(JNI.get_string_utfchars(str, iscopyptr))
end

#=
This method offers an intermediate representation
for arrays of java objects that does not require
full transformation to julia types
=#
function convert_to_vector(::Type{Vector{jobject}}, array::jobjectArray)
    vector = jobject[]
    len = JNI.get_array_length(array)
    for i in 1:len
        push!(vector, JNI.get_object_array_element(array, i-1))
    end
    vector
end

end


text:
module replacement
  using state, _deme_
  export  Replacement, replace!

  """
  The base type used by all replacement schemes.
  """
  abstract Replacement

  """
  Performs replacement on each deme within the population for a given algorithm
  state, according to a provided replacement scheme.
  """
  replace!(::Replacement, ::State) =
    error("No replace! operator implemented by this replacement scheme.")

  """
  Replaces the existing members of a deme with their offspring at the end of
  each generation.
  """
  type GenerationalReplacement <: Replacement; end

  """
  Generational replacement simply replaces the existing population with their
  offspring at the each of each generation. This is the default replacement
  scheme used by simple genetic algorithms. Also known as comma replacement,
  when used in conjuction with evolution strategies.
  """
  generational() = GenerationalReplacement()
  comma = generational

  replace!(r::GenerationalReplacement, s::State) =
    for d in s.population.demes; replace!(r, d); end

  function replace!(r::GenerationalReplacement, d::Deme)
    # Truncate the number of offspring.
    if d.num_offspring > d.capacity
      deleteat!(d.offspring.fitnesses, d.capacity:d.num_offspring)
      for stg in values(d.offspring.stages)
        deleteat!(stg, d.capacity:d.num_offspring)
      end
    end

    # Swap the members collection with the offspring collection.
    d.members = d.offspring
  end
end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

"""
    StackedTimeSolver

A module that is part of StateSpaceEcon package.
Contains methods for solving the dynamic system of equations
for the model and running simulations.
"""
module StackedTimeSolver

using SparseArrays
using LinearAlgebra
using ForwardDiff
using DiffResults
using Printf

using ModelBaseEcon
using TimeSeriesEcon

using ..Plans

include("stackedtime/abstract.jl")
include("stackedtime/fctypes.jl")
include("stackedtime/misc.jl")
include("stackedtime/solverdata.jl")
include("stackedtime/simulate.jl")

end # module

using .StackedTimeSolver

export FinalCondition
export NoFinalCondition
export HasFinalCondition
export FCNone, fcnone
export FCGiven, fcgiven
export FCMatchSSLevel, fclevel
export FCMatchSSRate, fcslope, fcrate
export FCConstRate, fcnatural
export setfc

export simulate
export seriesoverlay, dictoverlay
export dict2array, array2dict, array2data


text:
module NonlinearTestProblems

# export API
export domain_dimension, range_dimension, root, starting_point, lower_bounds, upper_bounds

# export problems
export F_NWp281, Rosenbrock, PowellSingular, PowellBadlyScaled, HelicalValley, Beale

using ArgCheck: @argcheck
using DocStringExtensions: FUNCTIONNAME, SIGNATURES

####
#### generic API
####
#### A *test problem* is a `ℝⁿ → ℝᵐ` mapping with some metadata (roots, starting point,
#### dimension).
####

"""
`$(FUNCTIONNAME)(problem)`

Dimension `n` of a the domain of a ``ℝⁿ → ℝᵐ`` test problem.
"""
function domain_dimension end

"""
$(SIGNATURES)

Dimension `m` of a the range of a ``ℝⁿ → ℝᵐ`` test problem.

!!! note
    Defaults to the domain dimension, only non-square implementations should need to define
    a method.
"""
range_dimension(problem) = domain_dimension(problem)

"""
`$(FUNCTIONNAME)(problem)`

Root (solution) of a test problem. Formally,

```julia
problem(root(problem)) ≈ zeros(range_dimension(problem))
```

Problems are guaranteed to have a unique root.
"""
function root end

"""
`$(FUNCTIONNAME)(problem)`

`$(FUNCTIONNAME)(problem, α)`

Recommended starting point of a test problem for methods that need one (eg quasi-Newton).
Usually taken from the literature that defines a problem, and should be “difficult”.

When used with a (positive real) argument `α`, the method should attempt to make it more
difficult (eg move away further from the root).

!!! note
    A default method is provided for the latter, problems just need to implement the
    single-argument version.
"""
function starting_point end

function starting_point(problem, α::Real)
    @argcheck α ≥ 0
    x0 = root(problem)
    x1 = starting_point(problem)
    @. x0 + (1 + α) * (x1 - x0)
end

"""
$(SIGNATURES)

Recommended lower bounds for box-constrained methods. Lower- and upper bounds are guaranteed
to contain the root. Coordinate-wise midpoints of bounds should not be near the root.

See also [`upper_bounds`](@ref).
"""
lower_bounds(problem) = fill(-100.0, domain_dimension(problem))

"""
$(SIGNATURES)

Recommended upper bounds for box-constrained methods. See [`lower_bounds`](@ref).
"""
upper_bounds(problem) = fill(50.0, domain_dimension(problem))

####
#### specific test problems
####

###
### Nocedal and Wright p 281
###

"Problem from Nocedal and Wright, p 281."
struct F_NWp281 end

domain_dimension(::F_NWp281) = 2

root(::F_NWp281) = [0, 1]

starting_point(::F_NWp281) = [-0.5, 1.4]

function (::F_NWp281)(x)
    x1, x2 = x
    [(x1 + 3) * (x2^3 - 7) + 18, sin(x2 * exp(x1) - 1)]
end

###
### Rosenbrock
###

"Rosenbrock function. (1) from Moré et al (1981)."
struct Rosenbrock end

domain_dimension(::Rosenbrock) = 2

root(::Rosenbrock) = [1, 1]

starting_point(::Rosenbrock) = [-1.2, 1]

function (::Rosenbrock)(x)
    x1, x2 = x
    [10 * (x2 - abs2(x1)), 1 - x1]
end

###
### Powell singular function
###

"Powell singular function. (13) from Moré et al (1981)."
struct PowellSingular end

domain_dimension(::PowellSingular) = 4

root(::PowellSingular) = zeros(4)

starting_point(::PowellSingular) = [3, -1, 0, 1]

function (::PowellSingular)(x)
    x1, x2, x3, x4 = x
    [x1 + 10 * x2, √5 * (x3 - x4), abs2(x2 - 2 * x3), √10 * abs2(x1 - x4)]
end

###
### Powell badly scaled function
###

"Powell's badly scaled function. (3) from Moré et al (1981)."
struct PowellBadlyScaled end

domain_dimension(::PowellBadlyScaled) = 2

root(::PowellBadlyScaled) = [1.0981593296999222e-5, 9.106146739865656]

starting_point(::PowellBadlyScaled) = [0, 1]

function (::PowellBadlyScaled)(x)
    x1, x2 = x
    [1e4 * x1 * x2 - 1, exp(-x1) + exp(-x2) - 1.0001]
end

###
### Helical valley function
###

"Helical valley function. (7) from Moré et al (1981)."
struct HelicalValley end

domain_dimension(::HelicalValley) = 3

root(::HelicalValley) = [1, 0, 0]

starting_point(::HelicalValley) = [-1, 0, 0]

function (::HelicalValley)(x)
    x1, x2, x3 = x
    θ = 1/(2π) * atan(x2 / x1) + ifelse(x1 < 0, 0.5, 0)
    [10*(x3 - 10 * θ), 10 * (hypot(x1, x2) - 1), x3]
end

###
### Beale function
###

"Beale function. (5) from Moré et al (1981)."
struct Beale end

domain_dimension(::Beale) = 2

range_dimension(::Beale) = 3

root(::Beale) = [3, 0.5]

starting_point(::Beale) = [1, 1]

function (::Beale)(x)
    [1.5, 2.25, 2.625] .- x[1] .* (1 .- x[2] .^ (1:3))
end

end # module


text:
"""
The basic algorithm where on every step of the iteration
there is a new item applying to the registry that is being
collectively judged by agents during the majority voting
"""
function simSimple(num_of_steps, agents, benchmarks)
    registry = []
    history = []

    for round in 1:num_of_steps
        application(registry, history, agents, simpleVote)
    end

    [b(registry, agents) for b in benchmarks]
end

"""
Apart from voting on a new appliations, a random agent challenges
the item that he consider the worst in the registry
"""
function simChallenge(num_of_steps, agents, benchmarks)
    registry = []
    history = []

    for round in 1:num_of_steps
        application(registry, history, agents, simpleVote)
        challenge(registry, agents, 0, simpleVote, noneRedistribution)
    end

    [b(registry, agents) for b in benchmarks]
end

"""
This mode introduces a token as a mean to reward well-performing agents.
During the challenge phase an agent needs to stake some tokens and depending
on the results he can either earn more tokens or loose the initial deposit.
The exact logic for token redistribution is encoded as a parameter passed
to the challenge method.
"""
function simToken(num_of_steps, agents, redistribution, benchmarks)
    registry = []
    history = []

    for round in 1:num_of_steps
        application(registry, history, agents, proRataTokenVote)
        challenge(registry, agents, 10, proRataTokenVote, redistribution)
    end

    [b(registry, agents) for b in benchmarks]
end


text:
module TestConstraintRelativeEntropyToExponential

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_RelativeEntropy()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.RelativeEntropy{Float64}(mock)
    MOI.Test.test_basic_VectorOfVariables_RelativeEntropyCone(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_basic_VectorAffineFunction_RelativeEntropyCone(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_basic_VectorQuadraticFunction_RelativeEntropyCone(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    entr1 = 2 * log(2)
    entr2 = 3 * log(3 / 5)
    var_primal = [entr1 + entr2, entr1, entr2]
    exps_duals = [[-1, log(0.5) - 1, 2], [-1, log(5 / 3) - 1, 0.6]]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
            (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone) =>
                exps_duals,
        )

    MOI.Test.test_conic_RelativeEntropyCone(bridged_mock, config)
    var_names = ["u"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    greater = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.GreaterThan{Float64},
        }(),
    )
    exps = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.ExponentialCone,
        }(),
    )
    (y1, y2) = MOI.get(mock, MOI.ListOfVariableIndices())[2:3]
    MOI.set(mock, MOI.VariableName(), y1, "y1")
    MOI.set(mock, MOI.VariableName(), y2, "y2")
    @test length(greater) == 1
    MOI.set(mock, MOI.ConstraintName(), greater[1], "greater")
    @test length(exps) == 2
    MOI.set(mock, MOI.ConstraintName(), exps[1], "exps1")
    MOI.set(mock, MOI.ConstraintName(), exps[2], "exps2")

    s = """
    variables: u, y1, y2
    greater: u + -1.0y1 + -1.0y2 >= 0.0
    exps1: [-1.0y1, 2.0, 1.0] in VecMathOptInterface.ExponentialCone()
    exps2: [-1.0y2, 3.0, 5.0] in VecMathOptInterface.ExponentialCone()
    minobjective: u
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        vcat(var_names, "y1", "y2"),
        ["greater", "exps1", "exps2"],
    )
    relentr = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RelativeEntropyCone,
        }(),
    )
    @test length(relentr) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), relentr[1], "relentr")

    s = """
    variables: u
    relentr: [u, 1.0, 5.0, 2.0, 3.0] in VecMathOptInterface.RelativeEntropyCone(5)
    minobjective: u
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(bridged_mock, model, var_names, ["relentr"])
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RelativeEntropyCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value =
            (attr isa MOI.ConstraintPrimalStart) ?
            [entr1 + entr2 + 1, 1, 5, 2, 3] :
            [2, 2, 0.6, log(0.5) - 1, log(5 / 3) - 1]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), y1) == entr1
            @test MOI.get(mock, MOI.VariablePrimalStart(), y2) == entr2
            @test MOI.get(mock, attr, greater[1]) == 1
            @test MOI.get(mock, attr, exps[1]) == [-entr1, 2, 1]
            @test MOI.get(mock, attr, exps[2]) == [-entr2, 3, 5]
        else
            @test MOI.get(mock, attr, greater[1]) == 2
            for i in 1:2
                @test MOI.get(mock, attr, exps[i]) ==
                      [value[1], value[3+i], value[1+i]]
            end
        end
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        1,
        (
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0),
        ),
    )
    return
end

end  # module

TestConstraintRelativeEntropyToExponential.runtests()


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

JLLWrappers.@generate_main_file_header("NodeJS_16")
JLLWrappers.@generate_main_file("NodeJS_16", UUID("a4b94fbf-73f5-58ff-888d-48a8396e17f6"))
end  # module NodeJS_16_jll


text:
#!/usr/bin/env julia

module YAMLTests

import YAML
import Base.Filesystem
using Test

const tests = [
    "spec-02-01",
    "spec-02-02",
    "spec-02-03",
    "spec-02-04",
    "spec-02-05",
    "spec-02-06",
    "spec-02-07",
    "spec-02-08",
    "spec-02-09",
    "spec-02-10",
    "spec-02-11",
    "spec-02-12",
    "spec-02-13",
    "spec-02-14",
    "spec-02-15",
    "spec-02-16",
    "spec-02-17",
    "spec-02-18",
    "spec-02-19",
    "spec-02-20",
    "spec-02-21",
    "spec-02-22",
    "spec-02-23",
    "empty_scalar",
    "no_trailing_newline",
    "windows_newlines",
    "escape_sequences",
    "issue15",
    "issue30",
    "issue36",
    "issue39",
    "cartesian",
    "ar1",
    "ar1_cartesian",
    "merge-01",
    "version-colon"
]

# ignore some test cases in write_and_load testing
const test_write_ignored = [
    "spec-02-17",
    "escape_sequences",
    "cartesian",
    "ar1",
    "ar1_cartesian"
]

function equivalent(xs::AbstractDict, ys::AbstractDict)
    if Set(collect(keys(xs))) != Set(collect(keys(ys)))
        @info "Not equivalent" Set(collect(keys(xs))) Set(collect(keys(ys)))
        return false
    end

    for k in keys(xs)
        if !equivalent(xs[k], ys[k])
            @info "Not equivalent" xs[k] ys[k]
            return false
        end
    end

    true
end

function equivalent(xs::AbstractArray, ys::AbstractArray)
    if length(xs) != length(ys)
        @info "Not equivalent" length(xs) length(ys)
        return false
    end

    for (x, y) in zip(xs, ys)
        if !equivalent(x, y)
            @info "Not equivalent" x y
            return false
        end
    end

    true
end

function equivalent(x::Float64, y::Float64)
    isnan(x) && isnan(y) ? true : x == y
end

function equivalent(x::AbstractString, y::AbstractString)
    while endswith(x, "\n")
        x = x[1:end-1] # trailing newline characters are ambiguous
    end
    while endswith(y, "\n")
        y = y[1:end-1]
    end
    x == y
end

function equivalent(x, y)
    x == y
end

# test custom tags
function construct_type_map(t::Symbol, constructor::YAML.Constructor,
                            node::YAML.Node)
    mapping = YAML.construct_mapping(constructor, node)
    mapping[:tag] = t
    mapping
end

const more_constructors = let
    pairs = [("!Cartesian", :Cartesian),
             ("!AR1", :AR1)]
    Dict{String,Function}([(t, (c, n) -> construct_type_map(s, c, n))
                           for (t, s) in pairs])
end

# write a file, then load its contents to be tested again
function write_and_load(data::Any)
    path = Filesystem.tempname() * ".yml" # path to a temporary file
    try
        YAML.write_file(path, data)
        return YAML.load_file(path, more_constructors)
    finally
        Filesystem.rm(path, force=true)
    end
end

const testdir = dirname(@__FILE__)
@testset for test in tests
    data = YAML.load_file(
        joinpath(testdir, string(test, ".data")),
        more_constructors
    )
    expected = evalfile(joinpath(testdir, string(test, ".expected")))
    @test equivalent(data, expected)
    if !in(test, test_write_ignored)
        @test equivalent(write_and_load(data), expected)
    else
        println("WARNING: I do not test the writing of $test")
    end
end

# test that an OrderedDict is written in the correct order
using OrderedCollections, DataStructures
@test strip(YAML.yaml(OrderedDict(:c => 3, :b => 2, :a => 1))) == join(["c: 3", "b: 2", "a: 1"], "\n")

# test that arbitrary dicttypes can be parsed
const dicttypes = [
    Dict{Any,Any},
    Dict{String,Any},
    Dict{Symbol,Any},
    OrderedDict{String,Any},
    () -> DefaultDict{String,Any}(Missing),
]
@testset for dicttype in dicttypes
    data = YAML.load_file(
        joinpath(testdir, "nested-dicts.data"),
        more_constructors;
        dicttype=dicttype
    )
    if typeof(dicttype) <: Function
        dicttype = typeof(dicttype())
    end # check the return type of function dicttypes
    _key(k::String) = keytype(dicttype) == Symbol ? Symbol(k) : k # String or Symbol key
    @test typeof(data) == dicttype
    @test typeof(data[_key("outer")]) == dicttype
    @test typeof(data[_key("outer")][_key("inner")]) == dicttype
    @test data[_key("outer")][_key("inner")][_key("something_unrelated")] == "1" # for completeness

    # type-specific tests
    if dicttype <: OrderedDict
        @test [k for (k,v) in data] == [_key("outer"), _key("anything_later")] # correct order
    elseif [k for (k,v) in data] == [_key("outer"), _key("anything_later")]
        @warn "Test of OrderedDict might not be discriminative: the order is also correct in $dicttype"
    end
    if dicttype <: DefaultDict
        @test data[""] === missing
    end
end

# also check that things break correctly
@test_throws YAML.ConstructorError YAML.load_file(
    joinpath(testdir, "nested-dicts.data"),
    more_constructors;
    dicttype=Dict{Float64,Any}
)

@test_throws YAML.ConstructorError YAML.load_file(
    joinpath(testdir, "nested-dicts.data"),
    more_constructors;
    dicttype=Dict{Any,Float64}
)

@test_throws ArgumentError YAML.load_file(
    joinpath(testdir, "nested-dicts.data"),
    more_constructors;
    dicttype=(mistaken_argument) -> DefaultDict{String,Any}(mistaken_argument)
)

@test_throws ArgumentError YAML.load_file(
    joinpath(testdir, "nested-dicts.data"),
    more_constructors;
    dicttype=() -> 3.0 # wrong type
)

# issue 81
dict_content = ["key1" => [Dict("subkey1" => "subvalue1", "subkey2" => "subvalue2"), Dict()], "key2" => "value2"]
order_one = OrderedDict(dict_content...)
order_two = OrderedDict(dict_content[[2,1]]...) # reverse order
@test YAML.yaml(order_one) != YAML.yaml(order_two)
@test YAML.load(YAML.yaml(order_one)) == YAML.load(YAML.yaml(order_two))

# issue 89 - quotes in strings
@test YAML.load(YAML.yaml(Dict("a" => """a "quoted" string""")))["a"] == """a "quoted" string"""
@test YAML.load(YAML.yaml(Dict("a" => """a \\"quoted\\" string""")))["a"] == """a \\"quoted\\" string"""

end  # module


text:
# include("../musclemodel/torque_generator.jl")
function getTorque(sol, t)
    @inbounds θhe, θke, θae, θhf, θkf, θaf = sol(t)[15:end]
    @unpack he, ke, ae, hf, kf, af = sol.prob.p

    τhe = torque(he, θhe)
    τke = torque(ke, θke)
    τae = torque(ae, θae)
    τhf = torque(hf, θhf)
    τkf = torque(kf, θkf)
    τaf = torque(af, θaf)

    return τhe, τke, τae, τhf, τkf, τaf
end

getTorque(sol) = reduce(hcat, [[getTorque(sol, t)...] for t in sol.t])'

function get_torque_generator(sol, t)
    @inbounds q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7, θhe, θke, θae, θhf, θkf, θaf = sol(t)
    @unpack he, hf, ke, kf, ae, af = sol.prob.p

    # update torque generators
    θh = π + q7
    θk = π - q6
    θa = π + q5
    θm = q4
    ωh = u7
    ωk = -u6
    ωa = u5
    ωm = u4

    # get values from torque generators
    d_he = get_torque_generator(he, t, 2π - θh, -ωh, θhe)
    d_ke = get_torque_generator(ke, t, 2π - θk, -ωk, θke)
    d_ae = get_torque_generator(ae, t, 2π - θa, -ωa, θae)
    d_hf = get_torque_generator(hf, t, θh, ωh, θhf)
    d_kf = get_torque_generator(kf, t, θk, ωk, θkf)
    d_af = get_torque_generator(af, t, θa, ωa, θaf)

    # CC torque relevant values
    out = map([d_he, d_hf, d_ke, d_kf, d_ae, d_af]) do d
        [d[:cc][:α], d[:cc][:τ_θ], d[:cc][:τ_ω], d[:cc][:α] * d[:cc][:τ_θ] * d[:cc][:τ_ω]]
    end

    return out
end

get_torque_generator(sol) = [get_torque_generator(sol, t) for t in sol.t]

function get_forces(sol, t)
    @inbounds q1, q2, q3, q4, q5, q6, q7, u1, u2, u3, u4, u5, u6, u7 = sol(t)
    @unpack k1, k2, k3, k4, k5, k6, k7, k8, l2, pop1xi, pop2xi, virtual_force = sol.prob.p

    pop2y = q2 - l2 * sin(q3 - q4 - q5 - q6 - q7)

    if q2 < 0.0
        dp1x = q1 - pop1xi
        ry1 = -k3 * q2 - k4 * abs(q2) * u2
        rx1 = ry1 * (-k1 * dp1x - k2 * u1)

    else
        rx1 = 0.0
        ry1 = 0.0
    end

    if pop2y < 0.0
        pop2x = q1 - l2 * cos(q3 - q4 - q5 - q6 - q7)
        vop2x = u1 - l2 * sin(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7)
        vop2y = u2 - l2 * cos(q3 - q4 - q5 - q6 - q7) * (u3 - u4 - u5 - u6 - u7)
        dp2x = pop2x - pop2xi
        ry2 = -k7 * pop2y - k8 * abs(pop2y) * vop2y
        rx2 = ry2 * (-k5 * dp2x - k6 * vop2x)
    else
        rx2 = 0.0
        ry2 = 0.0
    end

    vrx1::Float64 = virtual_force.vrx1(t)
    vry1::Float64 = virtual_force.vry1(t)
    vrx2::Float64 = virtual_force.vrx2(t)
    vry2::Float64 = virtual_force.vry2(t)

    return rx1, ry1, rx2, ry2, vrx1, vry1, vrx2, vry2
end

get_forces(sol) = [get_forces(sol, t) for t in sol.t]
rx(sol) = rx1(sol) .+ rx2(sol)
ry(sol) = ry1(sol) .+ ry2(sol)
rx1(sol) = [x[1] for x in get_forces(sol)]
ry1(sol) = [x[2] for x in get_forces(sol)]
rx2(sol) = [x[3] for x in get_forces(sol)]
ry2(sol) = [x[4] for x in get_forces(sol)]
vry(sol) = vry1(sol) .+ vry2(sol)
vrx(sol) = vrx1(sol) .+ vrx2(sol)
vrx1(sol) = [x[5] for x in get_forces(sol)]
vry1(sol) = [x[6] for x in get_forces(sol)]
vrx2(sol) = [x[7] for x in get_forces(sol)]
vry2(sol) = [x[8] for x in get_forces(sol)]

rx1(sol, t) = get_forces(sol, t)[1]
ry1(sol, t) = get_forces(sol, t)[2]
rx2(sol, t) = get_forces(sol, t)[3]
ry2(sol, t) = get_forces(sol, t)[4]
rx(sol, t) = rx1(sol, t) + rx2(sol, t)
ry(sol, t) = ry1(sol, t) + ry2(sol, t)

vrx1(sol, t) = get_forces(sol, t)[5]
vry1(sol, t) = get_forces(sol, t)[6]
vrx2(sol, t) = get_forces(sol, t)[7]
vry2(sol, t) = get_forces(sol, t)[8]
vrx(sol, t) = vrx1(sol, t) + vrx2(sol, t)
vry(sol, t) = vry1(sol, t) + vry2(sol, t)

RX(sol) = rx(sol) .+ vrx(sol)
RY(sol) = ry(sol) .+ vry(sol)
RX(sol, t) = rx(sol, t) + vrx(sol, t)
RY(sol, t) = ry(sol, t) + vry(sol, t)

## for stance simulation
# contact_time(sol) = sol.t[end]
# function aerial_time(sol)
#     @unpack g = sol.prob.p
#     tc = contact_time(sol)
#     vcmyto = vocmy(sol, tc)
#     discrim = vcmyto^2 + 4.0g * (pocmy(sol, 0.0) - pocmy(sol, tc))
#     discrim ≥ 0.0 ? ta = (-vcmyto - sqrt(discrim)) / g : ta = 100_000.0

#     return ta
# end

# # swing time
# swing_time(sol) = contact_time(sol) + 2aerial_time(sol)

# # step averaged velocity
# function step_velocity(sol)
#     tc = contact_time(sol)
#     ta = aerial_time(sol)
#     lc = pocmx(sol, tc) - pocmx(sol, 0.0)   # contact length
#     vcmyto = vocmy(sol, tc)
#     la = vocmx(sol, tc) * ta    # aerial length

#     return (lc + la) / (ta + tc) # step averaged velocity
# end

## for stride simulation
contact_idx(sv) = findfirst(x -> x ≥ 0.0, @view sv[2:end, 2])
aerial_idx(sv) = findfirst(x -> x ≥ 0.0, @view sv[:, 4]) - contact_idx(sv)
swing_time(sol, sv) = sol.t[contact_idx(sv)] + 2 * sol.t[aerial_idx(sv)]
stride_velocity(sol) = mean(vocmx(sol))


text:
x, y = generators(MPoly{Float64}, :x, :y)

@test (x+y)^2 == x^2 + 2x*y + y^2
@test x*(2x-y+1)^2 + 1 == 4x^3 - 4x^2*y + 4x^2 + y^2*x - 2x*y + x + 1
@test x + 1 == x + 1.0
@test x[1,0] == 1.0
@test x[0,1] == 0.0
@test eltype(x + im) == Complex{Float64}

@test  (x^(-1))[-1,0] == 1.0
@test_throws ArgumentError (x*y + x)^(-1)