text:
################################################################################
##################### MODULES FOR MULTI PARTICLE COLLISIONS ####################
################################################################################
using StatsBase
#This is to get the number of box where the particle is.
function get_box(parts::Array{particle,1}, boxes::Array{box,1}, Lx::Int64)
    for box in boxes
        box.np[:] = 0
    end
    for p in parts
        p.indbox = ceil(p.pgrid[1]) + Lx * (ceil(p.pgrid[2])-1)
        if p.indbox < 0
            println(p)
        end
        boxes[p.indbox].np[p.tp] += 1
    end
end
################################################################################
#rotation function
function rotate_vec(v::Array{Float64,1}, α::Float64)
    R = [cos(α) sin(α) ; -sin(α) cos(α)]
    return R*v
end
################################################################################
#computing the momentum of the boxes
function box_vel(parts::Array{particle,1},boxes::Array{box,1})
    for (i, box) in enumerate(boxes) #this is for enumerating the boxes
        box.vel = zeros(2)
        tmass = 0 #initializating the total mass of the box
        for p in filter(x-> x.indbox == i, parts) #the loop is made in the particles inside the box i
            #println("inside the filter ", i)
            box.vel += p.mass * p.vel #sums the momentums of all particles
            tmass += p.mass #sums the mass
        end
        if tmass != 0
            box.vel /= tmass #normalize the momentum of the box with the mass of the particles.
        end
    end
end
################################################################################
#This is for the collisions of the particles of equal mass or same particles
function box_velmc(parts::Array{particle,1},boxes::Array{box,1}, m::Array{Float64,1})
    for (i, box) in enumerate(boxes) #this is for enumerating the boxes
        if countnz(box.np) <= 1; continue; end
        for j in m #cycling over the different masses.
            box.vel = zeros(2)
            tmass = 0 #initializating the total mass of the box
            for p in filter(x-> x.indbox == i && x.mass == j , parts) #the loop is made in the particles inside the box i
                #println("inside the filter ", i)
                box.vel += p.mass * p.vel #sums the momentums of all particles
                tmass += p.mass #sums the mass
            end
            if tmass != 0
            box.vel /= tmass #normalize the momentum of the box with the mass of the particles.
            end

        end
    end
end
################################################################################
#collision in the new way.
function collide_mc(parts::Array{particle,1}, angles::Array{Float64})
    α = rand(angles)rand([-1,1])
    vel = zeros(2)
    tmass = 0 #initializating the total mass of the box
    for p in parts #the loop is made in the particles inside the box
        vel += p.mass * p.vel #sums the momentums of all particles
        tmass += p.mass #sums the mass
    end
    vel /= tmass #normalize the momentum of the box with the mass of the particles.
    for p in parts #loop over all particles
        v = p.vel - vel #extraction of the velocity of the box
        #α = 90.0*rand([-1,1]) #the rotation angle
        vn = rotate_vec(v, α) #rotation of the vector
        p.vel = vn + vel #adding the vector and the velocity of the box
    end
end
function collide_sc(parts::Array{particle,1}, tp::Int64, angles::Array{Float64})
    for j =1:tp
        α = rand(angles)*rand([-1,1])
        vel = zeros(2)
        tmass = 0
        sp = filter(x-> x.tp == j, parts)
        if isempty(sp); continue; end
        for p in sp
            vel += p.mass * p.vel
            tmass += p.mass
        end
        vel /= tmass
        for p in sp #loop over all particles
            v = p.vel - vel #extraction of the velocity of the box
            #α = 90.0*rand([-1,1]) #the rotation angle
            vn = rotate_vec(v, α) #rotation of the vector
            p.vel = vn + vel #adding the vector and the velocity of the box
        end
    end
end
#function of the shifting, actually you shift the positions of the particles.
function shift_grid!(parts::Array{particle,1},a::Float64, dim::Array{Int64,1})
    δx = rand()*rand(-a/2:a/2)
    δy = rand()*rand(-a/2:a/2)
    for p in parts
        if p.pos[1] + δx < 0 || p.pos[1] + δx > dim[1]
            p.pgrid[1] = p.pos[1]
        else
            p.pgrid[1] = p.pos[1] + δx
        end
        #p.pgrid[1] = mod(p.pos[1] + δx, dim[1])
        p.pgrid[2] = mod(p.pos[2] + δy, dim[2])
    end
end
#now we need to shift back the particles.
function shiftback_grid!(parts::Array{particle,1})
    for p in parts
        p.pgrid[1] = p.pos[1]
        p.pgrid[2] = p.pos[2]
    end
end
################################################################################
#computing the new velocities of the particles, collisions.
function parts_vels!(parts::Array{particle,1}, boxes::Array{box,1}, angles::Array{Float64,1})
    for p in parts #loop over all particles
        v = p.vel - boxes[p.indbox].vel #extraction of the velocity of the box
        α = rand(angles)*rand([-1,1]) #the rotation angle
        vn = rotate_vec(v, α) #rotation of the vector
        p.vel = vn + boxes[p.indbox].vel #adding the vector and the velocity of the box
    end
end
################################################################################
#getting new positions of the particles with periodic boundary conditions (pbc)
function getpos_pbc!(parts::Array{particle,1}, τ::Float64, dim::Array{Int64,1})
    for p in parts
        p.pos[1] = mod(p.pos[1] + p.vel[1] * τ , dim[1]) #to get the periodic boundary conditions.
        p.pos[2] = mod(p.pos[2] + p.vel[2] * τ , dim[2])
    end
end
#with slip walls
function getpos_slip!(parts::Array{particle,1}, τ::Float64, dim::Array{Int64,1})
    for p in  parts
        #doing the bouncing on the x axis
        if p.pos[1] + p.vel[1] * τ > dim[1]
            dif = p.pos[1] + p.vel[1] * τ - dim[1]
            p.pos[1] = dim[1] - dif #it moves back
            p.vel[1] = -p.vel[1] #changing direction of velocity in x
        elseif p.pos[1] + p.vel[1] * τ < 0
            dif = abs(p.pos[1] + p.vel[1] * τ)
            p.pos[1] = 0 + dif
            p.vel[1] = -p.vel[1] #changin direction of velocity in x
        else
            p.pos[1] = p.pos[1] + p.vel[1] * τ
        end
        #now the bouncing on the y axis
        if p.pos[2] + p.vel[2] * τ > dim[2]
            dif = p.pos[2] + p.vel[2] * τ - dim[2]
            p.pos[2] = dim[2] - dif
            p.vel[2] = -p.vel[2]
        elseif p.pos[2] + p.vel[2] * τ < 0
            dif = abs(p.pos[2] + p.vel[2] * τ)
            p.pos[2] = 0 + dif
            p.vel[2] = -p.vel[2]
        else
            p.pos[2] = p.pos[2] + p.vel[2] * τ
        end
    end
end
################################################################################
#this is the ininitialization part, first the normalization of the total momentum.
function norm_momentum!(parts::Array{particle,1})
    vt = zeros(2) #sum of momentum
    mt = 0 #total mass
    for p in parts #summing all velocities and masses
        vt += p.mass * p.vel
        mt += p.mass
    end
    vt /= mt #normalizing the sum
    for p in parts #applying the normalization of momentum
        p.vel = p.vel - vt
    end
end
################################################################################
#now the normalization of temperature
function norm_temperature!(parts::Array{particle,1}, Tr::Float64)
    T = 0 #temperature
    for p in parts #this is the calculation of the Temperature
        T += p.mass * norm(p.vel)^2 / (2 * length(parts))
    end
    for p in parts
        p.vel = p.vel * sqrt(Tr / T) #multiplying each velocity for the factor of sqrt(Ttarget/Tactual)
    end
end
################################################################################
#this function compute probabilities of P(n(t+1)|n(t)) and returns how many particles are going to react.
function prob_box(np::Array{Int64}, kr::Float64)
    nr = min(np[1],np[2]) #estimating minimum or maximum of reactions
    if nr != 0
        p = zeros(nr+1) #array of probabilities
        p[1] = 1 - kr #no-reaction case
        for i=1:nr
            p[i+1] = kr * ( factorial(big(np[1])) / factorial(big(np[1])-i) * factorial(big(np[2])) / factorial(big(np[2])-i) * factorial(big(np[3])) / factorial(big(np[3])+i) )
        end
        s = sum(p)
        if s != 1
            p = p / s
        end
        cs = cumsum(p); a = rand()
        out = findfirst(sort([a;cs]),a)-1
    else
        out = 0
    end
    return out
end
################################################################################
function col_box(box::box,parts::Array{particle,1}, m::Array{Float64,1})
    parbox = filter(x-> x.indbox == box.ind, parts) #selecting the particles that are in the box.
    #println(i,'\t',length(parbox))
    if isempty(parbox) == false #ignore next steps if the box is empty
        collide_mc(parbox)
        if countnz(box.np) > 1
            collide_sc(parbox, m)
        end
    end
end
################################################################################
#next function is for the birth and death process of the reaction.
function reac_box(parts::Array{particle,1}, nr::Int64)
    pa = filter(x -> x.tp == 1, parts)
    pb = filter(x -> x.tp == 2, parts)
    pc = Array{particle,1}()
    for i =1:nr #nr = number of reactions
        a = rand(pa); b = rand(pb) #choosing particles a and b at random
        #println(a,b)
        cm = a.mass + b.mass #c mass
        #println(cm)
        x = a.mass/cm * a.pos[1] + b.mass/cm * b.pos[1] #position of new particle
        y = a.mass/cm * a.pos[2] + b.mass/cm * b.pos[2]
        vx = a.mass/cm * a.vel[1] + b.mass/cm * b.vel[1] #velocity of the new particle, conserving the momentum
        vy = a.mass/cm * a.vel[2] + b.mass/cm * b.vel[2]
        c = particle([1,1], cm, 3) #creating the particle c
        c.pos = [x,y]; c.vel = [vx,vy] #positions and velocities.
        push!(pc, c) #adding the particle c to the
        a.mass = 0.0; b.mass = 0.0
        filter!(x -> x.mass != 0.0, pa); filter!(x -> x.mass != 0.0, pb)  #removing the particles that reacted.
    end
    return pc
end
################################################################################
#next function is to graph the positions
function grap_pos(parts::Array{particle,1},tp::Int64)
    x = filter(x -> x.tp == tp, parts)
    xy = zeros(length(x),2)
    for (i,p) in enumerate(x)
        xy[i,:] = collect(p.pos)
    end
    return xy
end
################################################################################
#This is the nucleation function
function nucleate(parts::Array{particle,1}, ks::Int64)
    pc = filter(x -> x.tp == 3, parts) #choose the c-type particles
    dm = pc[1].mass * ks #calculate mass of d
    nd = div(length(pc),ks) #the number of nucleations
    d = [particle([1,1], dm, 4) for _ in 1:nd] #the array of d particles
    for k = 1:nd
        rc = sample(pc, ks, replace=false) #choosing randomly the c particles
        d[k].pos = 1/ks * mapreduce(x -> x.pos,+ , rc); d[k].vel = 1/ks * mapreduce(x -> x.vel,+ , rc) #calculation of the position
        for p in rc; p.mass = 0.0; end  #removing the c particles that nucleate
        filter!(x -> x.mass != 0.0, pc)
    end
    return d # returning the array of d particles.
end
################################################################################
#function for agregation in the same box.
function agg_sb!(parts::Array{particle,1})
    pc = filter(x -> x.tp == 3, parts)
    pd = filter(x -> x.tp == 4, parts)
    d = rand(pd)
    dm = pc[1].mass * length(pc) + d.mass;
    pos = 1 / dm * ( d.mass * d.pos + pc[1].mass * mapreduce(x -> x.pos,+ , pc)); vel = 1 / dm * (d.mass * d.vel + pc[1].mass * mapreduce(x -> x.vel,+ , pc))
    d.pos = pos; d.vel = vel
    for p in pc; p.mass = 0.0; end
end


text:
# File: getXLData.jl
# Author:
# Date: 29-June-2020
using XLSX
# GET XL DATA
function getXLData(fname::AbstractString, sname::AbstractString)
    # raw data
    rwdata = xlsread(fname,sname);
    # labels
    names = rwdata(1:1:end)
    names(1) = [];
    # dates
    dates = rwdata(2:1:)
    # returns
    rwdata(:,1) = [];
    returns = rwdata(2::);
    returns = cell2mat(returns);
end

ext: jl


text:
# plot/z-list.jl

include("jim.jl")
include("setup.jl")


text:
@testset "test output api" begin
    @testset "24-bus rts case" begin
        result = run_opf("../test/data/case24.m", ACPPowerModel, ipopt_solver)

        @test haskey(result, "solver") == true
        @test haskey(result, "status") == true
        @test haskey(result, "objective") == true
        @test haskey(result, "objective_lb") == true
        @test haskey(result, "solve_time") == true
        @test haskey(result, "machine") == true
        @test haskey(result, "data") == true
        @test haskey(result, "solution") == true
        @test haskey(result["solution"], "branch") == false

        @test !isnan(result["solve_time"])

        @test length(result["solution"]["bus"]) == 24
        @test length(result["solution"]["gen"]) == 33
    end
end

@testset "test branch flow output" begin
    @testset "24-bus rts case ac opf" begin
        result = run_opf("../test/data/case24.m", ACPPowerModel, ipopt_solver; setting = Dict("output" => Dict("branch_flows" => true)))

        @test haskey(result, "solver") == true
        @test haskey(result, "status") == true
        @test haskey(result, "objective") == true
        @test haskey(result, "objective_lb") == true
        @test haskey(result, "solve_time") == true
        @test haskey(result, "machine") == true
        @test haskey(result, "data") == true
        @test haskey(result, "solution") == true
        @test haskey(result["solution"], "branch") == true

        @test length(result["solution"]["bus"]) == 24
        @test length(result["solution"]["gen"]) == 33
        @test length(result["solution"]["branch"]) == 38

        branches = result["solution"]["branch"]

        @test isapprox(branches["2"]["pf"],  0.2001; atol = 1e-3)
        @test isapprox(branches["2"]["pt"], -0.1980; atol = 1e-3)
        @test isapprox(branches["2"]["qf"],  0.0055; atol = 1e-3)
        @test isapprox(branches["2"]["qt"], -0.0571; atol = 1e-3)
    end

    # A DCPPowerModel test is important because it does have variables for the reverse side of the branchs
    @testset "3-bus case dc opf" begin
        result = run_opf("../test/data/case3.m", DCPPowerModel, ipopt_solver; setting = Dict("output" => Dict("branch_flows" => true)))

        @test haskey(result, "solution") == true
        @test haskey(result["solution"], "branch") == true

        @test length(result["solution"]["bus"]) == 3
        @test length(result["solution"]["gen"]) == 3
        @test length(result["solution"]["branch"]) == 3
        @test length(result["solution"]["dcline"]) == 1

        branches = result["solution"]["branch"]

        @test isapprox(branches["3"]["pf"], -0.103497; atol = 1e-3)
        @test isapprox(branches["3"]["pt"],  0.103497; atol = 1e-3)
        @test isnan(branches["3"]["qf"])
        @test isnan(branches["3"]["qt"])
    end
end

@testset "test dual value output" begin
    settings = Dict("output" => Dict("duals" => true))
    result = run_dc_opf("../test/data/case14.m", ipopt_solver, setting = settings)

    PowerModels.make_mixed_units(result["solution"])
    @testset "14 bus - kcl duals" begin
        for (i, bus) in result["solution"]["bus"]
            @test haskey(bus, "lam_kcl_r")
            @test haskey(bus, "lam_kcl_i")
            @test isapprox(bus["lam_kcl_r"], -39.02; atol = 1e-2) # Expected result for case14
            @test isnan(bus["lam_kcl_i"])
        end
    end

    @testset "14 bus - thermal limit duals" begin
        for (i, branch) in result["solution"]["branch"]
            @test haskey(branch, "mu_sm_fr")
            @test haskey(branch, "mu_sm_to")
            @test isapprox(branch["mu_sm_fr"], 0.0; atol = 1e-2)
            @test isnan(branch["mu_sm_to"])
        end
    end

    result = run_dc_opf("../test/data/case5.m", ipopt_solver, setting = settings)

    PowerModels.make_mixed_units(result["solution"])
    @testset "5 bus - kcl duals" begin
        for (i, bus) in result["solution"]["bus"]
            @test bus["lam_kcl_r"] <= -9.00
            @test bus["lam_kcl_r"] >= -45.00
        end
    end

    @testset "5 bus - thermal limit duals" begin
        for (i, branch) in result["solution"]["branch"]
            if i != "7"
                @test isapprox(branch["mu_sm_fr"], 0.0; atol = 1e-2)
            else
                @test isapprox(branch["mu_sm_fr"], 54.70; atol = 1e-2)
            end
            @test isnan(branch["mu_sm_to"])
        end
    end
end

# recomended by @lroald
@testset "test solution feedback" begin

    function solution_feedback(case, ac_opf_obj)
        data = PowerModels.parse_file(case)
        opf_result = run_ac_opf(data, ipopt_solver)
        @test opf_result["status"] == :LocalOptimal
        @test isapprox(opf_result["objective"], ac_opf_obj; atol = 1e0)

        PowerModels.update_data(data, opf_result["solution"])

        pf_result = run_ac_pf(data, ipopt_solver)
        @test pf_result["status"] == :LocalOptimal
        @test isapprox(pf_result["objective"], 0.0; atol = 1e-3)

        for (i,bus) in data["bus"]
            @test isapprox(opf_result["solution"]["bus"][i]["va"], pf_result["solution"]["bus"][i]["va"]; atol = 1e-3)
            @test isapprox(opf_result["solution"]["bus"][i]["vm"], pf_result["solution"]["bus"][i]["vm"]; atol = 1e-3)
        end

        for (i,gen) in data["gen"]
            @test isapprox(opf_result["solution"]["gen"][i]["pg"], pf_result["solution"]["gen"][i]["pg"]; atol = 1e-3)
            # cannot check this value solution does not appeat to be unique; verify this!
            #@test isapprox(opf_result["solution"]["gen"][i]["qg"], pf_result["solution"]["gen"][i]["qg"]; atol = 1e-3)
        end

        for (i,dcline) in data["dcline"]
            @test isapprox(opf_result["solution"]["dcline"][i]["pf"], pf_result["solution"]["dcline"][i]["pf"]; atol = 1e-3)
            @test isapprox(opf_result["solution"]["dcline"][i]["pt"], pf_result["solution"]["dcline"][i]["pt"]; atol = 1e-3)
        end
    end

    @testset "3-bus case" begin
        solution_feedback("../test/data/case3.m", 5907)
    end

    @testset "5-bus asymmetric case" begin
        solution_feedback("../test/data/case5_asym.m", 17551)
    end

    @testset "5-bus with dcline costs" begin
        solution_feedback("../test/data/case5_dc.m", 17760.2)
    end

end


text:
module DemoPackageBXW
print_greeting() = print("Hello, world!")
# Write your package code here.

end


text:
function besselj(m, x)
    sn = 0
    an = 1/factorial(m)
    k = 0

    while abs(an) > eps()
        sn += an
        an = -an * x * x/4/(k+1)/(k+m+1)
        k += 1
    end

    sn *= (x/2)^m
    return sn
end

function besselj_plots()
    x = 0:0.02:20
    y1, y2, y3, y4, y5, y6 =[], [], [], [], [], []

    for i in x
        push!(y1, besselj(0, i))
        push!(y2, besselj(1, i))
        push!(y3, besselj(2, i))
        push!(y4, besselj(3, i))
        push!(y5, besselj(4, i))
        push!(y6, besselj(5, i))
    end

    p = plot(x, y1)
    plot!(p, x, y2)
    plot!(p, x, y3)
    plot!(p, x, y4)
    plot!(p, x, y5)
    plot!(p, x, y6)
    display(p)
end


text:
struct Backtest{Fl <: AbstractFloat}
    abs_errors::Matrix{Fl}
    mae::Vector{Fl}
    crps_scores::Matrix{Fl}
    mean_crps::Vector{Fl}
    function Backtest{Fl}(n::Int, steps_ahead::Int) where Fl
        abs_errors = Matrix{Fl}(undef, steps_ahead, n)
        crps_scores = Matrix{Fl}(undef, steps_ahead, n)
        mae = Vector{Fl}(undef, steps_ahead)
        mean_crps = Vector{Fl}(undef, steps_ahead)
        return new(abs_errors, mae, crps_scores, mean_crps)
    end
end

discrete_crps_indicator_function(val::Fl, z::Fl) where {Fl} = val < z
function crps(val::Fl, scenarios::Vector{Fl}) where {Fl}
    sorted_scenarios = sort(scenarios)
    m = length(scenarios)
    crps_score = zero(Fl)
    for i = 1:m
        crps_score +=
            (sorted_scenarios[i] - val) *
            (m * discrete_crps_indicator_function(val, sorted_scenarios[i]) - i + 0.5)
    end
    return (2 / m^2) * crps_score
end
evaluate_abs_error(y::Vector{Fl}, forecast::Vector{Fl}) where Fl = abs.(y - forecast)
function evaluate_crps(y::Vector{Fl}, scenarios::Matrix{Fl}) where {Fl}
    crps_scores = Vector{Fl}(undef, length(y))
    for k = 1:length(y)
        crps_scores[k] = crps(y[k], scenarios[k, :])
    end
    return crps_scores
end

"""
    backtest(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
             n_scenarios::Int = 10_000,
             filter::KalmanFilter=default_filter(model),
             optimizer::Optimizer=default_optimizer(model)) where Fl

Makes rolling window estimating and forecasting to benchmark the forecasting skill of the model
in for different time periods and different lead times. The function returns a struct with the MAE
and mean CRPS per lead time. See more on [Backtest the forecasts of a model](@ref)

# References
 * DTU course "31761 - Renewables in electricity markets" available on youtube https://www.youtube.com/watch?v=Ffo8XilZAZw&t=556s
"""
function backtest(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
                  n_scenarios::Int = 10_000,
                  filter::KalmanFilter=default_filter(model),
                  optimizer::Optimizer=default_optimizer(model))
    Fl = typeof_model_elements(model)
    num_mle = length(model.system.y) - start_idx - steps_ahead
    b = Backtest{Fl}(num_mle, steps_ahead)
    for i in 1:num_mle
        println("Backtest: step $i of $num_mle")
        y_to_fit = model.system.y[1:start_idx - 1 + i]
        y_to_verify = model.system.y[start_idx + i:start_idx - 1 + i + steps_ahead]
        model_to_fit = reinstantiate(model, y_to_fit)
        fit!(model_to_fit; filter=filter, optimizer=optimizer)
        forec = forecast(model_to_fit, steps_ahead; filter=filter)
        scenarios = simulate_scenarios(model_to_fit, steps_ahead, n_scenarios; filter=filter)
        expected_value_vector = forecast_expected_value(forec)[:]
        abs_errors = evaluate_abs_error(y_to_verify, expected_value_vector)
        crps_scores = evaluate_crps(y_to_verify, scenarios[:, 1, :])
        b.abs_errors[:, i] = abs_errors
        b.crps_scores[:, i] = crps_scores
    end
    for i in 1:steps_ahead
        b.mae[i] = mean(b.abs_errors[i, :])
        b.mean_crps[i] = mean(b.crps_scores[i, :])
    end
    return b
end


text:
# the parameters are:
#   B Boundary condition
#   S Storage data type
#   W ContWaveClass
#   N the type of value of the normalization
#   isAn boolean for whether this is an analytic transform or not
struct CWT{B,S,W<:ContWaveClass,N,isAn} <: ContWave{B,S}
    Q::S # the number of wavelets per octave, ie the scaling is s=2^(j/scalingfactor)
    β::S # the amount that Q decreases per octave
    fourierFactor::S
    coi::S
    α::Int   # the order for a Paul and the number of derivatives
    # for a DOG
    σ::Array{S} # the morlet wavelet parameters
    # (σ,κσ,cσ). NaN if not morlet.
    waveType::W        # because multiple dispatch is good TODO this is actually redundant, make a function called waveType that extracts this from above
    extraOctaves::S # adjust the last frequency used, in case you need frequencies beyond where the wavelets are well defined.
    averagingLength::S   # the number of scales to override with averaging.
    averagingType::Average # either Dirac or mother; the first uniformly
    # represents the lowest frequency information,
    # while the second constructs
    # a wavelet using Daughter that has mean frequency
    # zero, and std equal to
    # the first non-removed wavelet's mean
    frameBound::S       # if positive, set the frame bound of the
    # transform to be frameBound. Otherwise leave it
    # so that each wavelet has an L2 norm of 1
    p::N                     # the normalization that is preserved for the
    # wavelets as the scale changes. The conjugate
    # p-norm is preserved for the  signal. Should be
    # larger than 1
end

# aliased = ((:Q,:s,:scalingFactor), (:β,:decreasing), (:p, :normalization))
function processKeywordArgs(Q, β, p; kwargs...)
    keysVarg = keys(kwargs)
    if :s in keysVarg
        Q = kwargs[:s]
        if :scalingFactor in keysVarg
            error("you used two names for Q")
        end
    elseif :scalingFactor in keysVarg
        Q = kwargs[:scalingFactor]
    end

    if :decreasing in keysVarg
        β = kwargs[:decreasing]
    end

    if :normalization in keysVarg
        p = kwargs[:normalization]
    end
    return Q, β, p
end

@doc """
    CWT(wave::ContWaveClass, Q=8, boundary::WaveletBoundary=SymBoundary(),
    averagingType::Average = Father(), averagingLength::Int = 4, frameBound=1, p::N=Inf, β=4)
"""
function CWT(wave::WC, Q = 8, boundary::B = DEFAULT_BOUNDARY,
    averagingType::A = Father(),
    averagingLength::Real = 0,
    frameBound = 1, p::N = Inf,
    β = 4; extraOctaves = 0, kwargs...) where {WC<:ContWaveClass,A<:Average,B<:WaveletBoundary,N<:Real}
    Q, β, p = processKeywordArgs(Q, β, p; kwargs...) # some names are redundant
    @assert β > 0
    @assert p >= 1
    nameWavelet = name(wave)[1:3]
    tdef = calculateProperties(wave)

    if typeof(averagingType) <: NoAve
        averagingLength = 0
        averagingType = NoAve()
    end

    # S is the most permissive type of the listed variables
    S = promote_type(typeof(Q), typeof(β),
        typeof(frameBound), typeof(p),
        typeof(tdef[1]), typeof(tdef[2]))
    return CWT{B,S,WC,N,isAnalytic(wave)}(S(Q), S(β), tdef..., S(extraOctaves),
        S(averagingLength), averagingType, S(frameBound),
        S(p))
end

function calculateProperties(w::Morlet)
    σ = w.σ
    fourierFactor = (4 * π) / (σ + sqrt(2 + σ^2))
    coi = 1 / sqrt(2)
    α = -1
    σ = [w.σ, w.κσ, w.cσ]
    return fourierFactor, coi, α, σ, w
end

function calculateProperties(w::Dog)
    α = order(w)
    fourierFactor = 2 * π * sqrt(2 / (2 * α + 1))
    coi = 1 / sqrt(2)
    σ = [NaN]
    return fourierFactor, coi, α, σ, w
end

function calculateProperties(w::Paul)
    α = order(w)
    fourierFactor = 4 * π / (2 * α + 1)
    coi = sqrt(2)
    σ = [NaN]
    return fourierFactor, coi, α, σ, w
end

function calculateProperties(w::ContOrtho)
    α = -1
    fourierFactor = NaN
    coi = NaN
    σ = [NaN]
    return fourierFactor, coi, α, σ, w
end

function calculateProperties(w::Morse)
    α = -1
    fourierFactor = NaN
    coi = NaN
    σ = [NaN]
    return fourierFactor, coi, α, σ, w
end

name(s::CWT) = name(s.waveType)

function eltypes(::CWT{B,T,W,N}) where {B,T,W,N}
    T
end
function boundaryType(::CWT{B,T,W,N}) where {B,T,W,N}
    B
end
function waveletType(::CWT{B,T,W,N}) where {B,T,W,N}
    W
end

function Base.show(io::IO, cf::CWT{W,S,WT,N}) where {W,S,WT,N}
    print("CWT{$(cf.waveType), $(cf.averagingType), Q=$(cf.Q), β=$(cf.β)," *
          "aveLen=$(cf.averagingLength), frame=" * "$(cf.frameBound), norm=$(cf.p), extraOctaves=$(cf.extraOctaves)}")
end

"""
    wavelet(wave::ContWaveClass; Q=8, boundary::WaveletBoundary=DEFAULT_BOUNDARY,
    averagingType::Average = Father(), averagingLength = 4,
    frameBound=1, p=Inf, β=4, kwargs...)
A constructor for the `CWT` type, using keyword rather than positional options.
"""
function wavelet(wave::WC; Q = 8, boundary::T = DEFAULT_BOUNDARY,
    averagingType::A = Father(), averagingLength = 0,
    frameBound = 1, p::N = Inf, β = 4,
    kwargs...) where {WC<:ContWaveClass,A<:Average,T<:WaveletBoundary,N<:Real}
    return CWT(wave, Q, boundary, averagingType, averagingLength, frameBound,
        p, β; kwargs...)
end


text: f2(x,y) = x + y


text:
if VERSION >= v"0.7.0-"
    # Adding Pkg in test/REQUIRE would be an error in 0.6.  Using
    # Project.toml still has some gotchas.  So:
    Pkg = Base.require(Base.PkgId(Base.UUID(0x44cfe95a1eb252eab672e2afdf69b78f), "Pkg"))
end

# Let PyCall.jl use Python interpreter from Conda.jl
# See: https://github.com/JuliaPy/PyCall.jl
ENV["PYTHON"] = ""
Pkg.build("PyCall")

using Conda
Conda.add("scipy")


text:
function add_node(obj::Mesh, arg0::jint, arg1::jint)
    return jcall(obj, "addNode", Mesh_Node, (jint, jint), arg0, arg1)
end

function get_across_gap(obj::Mesh)
    return jcall(obj, "getAcrossGap", jdouble, ())
end

function get_along_gap(obj::Mesh)
    return jcall(obj, "getAlongGap", jdouble, ())
end

function get_closest_existing_node(obj::Mesh, arg0::Vector3D)
    return jcall(obj, "getClosestExistingNode", Mesh_Node, (Vector3D,), arg0)
end

function get_closest_existing_node(obj::Mesh, arg0::jint, arg1::jint)
    return jcall(obj, "getClosestExistingNode", Mesh_Node, (jint, jint), arg0, arg1)
end

function get_coverage(obj::Mesh)
    return jcall(obj, "getCoverage", SphericalPolygonsSet, ())
end

function get_inside_nodes(obj::Mesh)
    return jcall(obj, "getInsideNodes", List, ())
end

function get_max_across_index(obj::Mesh)
    return jcall(obj, "getMaxAcrossIndex", jint, ())
end

function get_max_along_index(obj::Mesh)
    return jcall(obj, "getMaxAlongIndex", jint, ())
end

function get_max_along_index(obj::Mesh, arg0::jint)
    return jcall(obj, "getMaxAlongIndex", jint, (jint,), arg0)
end

function get_min_across_index(obj::Mesh)
    return jcall(obj, "getMinAcrossIndex", jint, ())
end

function get_min_along_index(obj::Mesh)
    return jcall(obj, "getMinAlongIndex", jint, ())
end

function get_min_along_index(obj::Mesh, arg0::jint)
    return jcall(obj, "getMinAlongIndex", jint, (jint,), arg0)
end

function get_node(obj::Mesh, arg0::jint, arg1::jint)
    return jcall(obj, "getNode", Mesh_Node, (jint, jint), arg0, arg1)
end

function get_number_of_nodes(obj::Mesh)
    return jcall(obj, "getNumberOfNodes", jint, ())
end

function get_taxicab_boundary(obj::Mesh, arg0::jboolean)
    return jcall(obj, "getTaxicabBoundary", List, (jboolean,), arg0)
end


text:
if ccall(:jl_generating_output, Cint, ()) == 1
    let
        images2d = [rand(Float32, 100, 100),
                    rand(Float64, 100, 100),
                    rand(Gray{N0f8}, 100, 100),
                    rand(Gray{N0f16}, 100, 100),
                    rand(RGB{N0f8}, 100, 100),
        ]
        images3d = [rand(Gray{N0f16}, 100, 100, 10),
        ]

        for img in images2d
            for kern in (Kernel.gaussian((3,3)),
                        KernelFactors.gaussian((3,3)),
                        KernelFactors.IIRGaussian((3.0, 3.0)))
                for r in (CPU1(FIR()), CPUThreads(FIRTiled((5,5))), CPU1(FFT()))
                    @assert precompile(imfilter, (typeof(r), typeof(img), typeof(kern)))
                end
                @assert precompile(imfilter, (typeof(img), typeof(kern)))
            end
            if eltype(img) <: Union{Number,Gray}
                @assert precompile(mapwindow, (typeof(extrema), typeof(img), typeof((3,3))))
                @assert precompile(mapwindow, (typeof(median!), typeof(img), typeof((3,3))))
            end
        end
        for img in images3d
            for kern in (Kernel.gaussian((3,3,3)),
                         KernelFactors.gaussian((3,3,3)),
                         KernelFactors.IIRGaussian((3.0, 3.0, 3.0)))
                for r in (CPU1(FIR()), CPUThreads(FIRTiled((5,5))), CPU1(FFT()))
                    @assert precompile(imfilter, (typeof(r), typeof(img), typeof(kern)))
                end
                @assert precompile(imfilter, (typeof(img), typeof(kern)))
            end
        end
    end
end


text:
simloglik!(grad, grp::ObservationGroupEmpty, theta, sim, dograd) = nothing
grad_simloglik!(grad, grp::ObservationGroupEmpty, theta, sim) = nothing
update!(d::Union{EmptyDataSet,ObservationGroupEmpty}, args...) = nothing

# SML for unit i
function simloglik!(grad::AbstractVector, grptup::NTuple{N,ObservationGroup},
    thetas, idxs, sim::SimulationDrawsVector, dograd::Bool; kwargs...
) where {N}

    fill(_qm(sim), 0)
    logM = log(_num_sim(sim))
    zips = zip(idxs, thetas, grptup)

    for (idx,theta,grp) in zips
        @views simloglik!(grad[idx], grp, theta, sim, dograd; kwargs...)
    end

    LL = logsumexp!(_llm(sim)) - logM

    if dograd
        for (idx,theta,grp) in zips
            @views grad_simloglik!(grad[idx], grp, theta, sim)
        end
    end

    return LL
end

# SML for entire dataset
function simloglik!(grad::Vector, hess::Matrix, tmpgrad::Matrix, data::DataSetofSets,
    theta::AbstractVector, sim::SimulationDrawsMatrix, dograd::Bool; kwargs...
)
    nparm, num_i = size(tmpgrad)
    nparm == length(grad) == checksquare(hess) || throw(DimensionMismatch())

    # parameters
    thetasvw = split_thetas(data, theta)
    idxs = theta_indexes(data)

    # do updates
    rhoparm = theta_ρ(data,theta)
    update!(sim, rhoparm)
    for (dat,θ) in zip(data,thetasvw)
        update!(dat, θ)
    end

    LL = 0.0
    fill!(tmpgrad, 0)

    # iterate over observations
    for i = OneTo(num_i)
        gradi = view(tmpgrad, :, i)
        grptup = getindex.(data, i)
        simi = view(sim, i)
        fill!(_qm(sim), 0)
        LL += simloglik!(gradi, grptup, thetasvw, idxs, simi, dograd; kwargs...)
    end

    if dograd
        mul!(hess, tmpgrad, tmpgrad')
        sum!(reshape(grad, :, 1), tmpgrad)
        grad .*= -1
        hess .*= -1
    end

    return -LL
end


text:
test1 = ex1()
#nlp = convert(AbstractNLPModel, test1)
nlp = NLMPCC(test1)

@test nlp.mod.meta.ncc == 1

@test nlp.meta.nvar == 2
@test nlp.meta.lvar == [-Inf, -Inf]
@test nlp.meta.uvar == [Inf, Inf]

@test nlp.meta.ncon == 4
@test nlp.meta.lcon == [-1.0, -0.0, -0.0, -Inf]
@test nlp.meta.ucon == [Inf, Inf, Inf, 0.0]

@test obj(nlp, nlp.meta.x0) == 0.0
@test grad(nlp, nlp.meta.x0) == grad(test1,[Inf,Inf])
@test hess(nlp, nlp.meta.x0) == sparse(zeros(2,2))
@test hess(nlp, nlp.meta.x0, obj_weight = 0.0) == sparse(zeros(2,2))

t0 = ones(2)
@test cons(nlp, t0) == [-1.0, 1.0, 1.0, 1.0]
@test jac(nlp, t0) == [0.0 -1.0 ; 1.0 0.0 ; 0.0 1.0 ; 1.0 1.0]
@test jac_coord(nlp, t0) == jac(nlp, t0)[:]
@test jac_structure(nlp) == ([1, 2, 3, 4, 1, 2, 3, 4], [1, 1, 1, 1, 2, 2, 2, 2])
@test jprod(nlp, t0, zeros(2))  == zeros(4)
@test jtprod(nlp, t0, zeros(4)) == zeros(2)
@test hess(nlp, t0, ones(4)) == tril([0.0 1.0 ; 1.0 0.0])
@test length(hess_coord(nlp, t0)) == 3
@test length(hess_coord(nlp, t0, ones(4))) == 3
@test hess_structure(nlp) == ([1, 2, 2], [1, 1, 2])
@test hprod(nlp, t0, ones(4), nlp.meta.x0) == [1.0, 1.0]
@test hprod(nlp, t0, nlp.meta.x0) == [0.0, 0.0]

@test viol(nlp, t0) == [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]


text:
# ---
# title: 152. Maximum Product Subarray
# id: problem152
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/maximum-product-subarray/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, find the contiguous subarray within an array
# (containing at least one number) which has the largest product.
#
# **Example 1:**
#
#
#
#     Input: [2,3,-2,4]
#     Output: 6
#     Explanation:  [2,3] has the largest product 6.
#
#
# **Example 2:**
#
#
#
#     Input: [-2,0,-1]
#     Output: 0
#     Explanation:  The result cannot be 2, because [-2,-1] is not a subarray.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module LanguageServer
import URIParser
using JSON, REPL, CSTParser, DocumentFormat, SymbolServer, StaticLint
using CSTParser: EXPR, Tokenize.Tokens, typof, kindof, parentof, valof
using StaticLint: refof, scopeof, bindingof
using UUIDs
export LanguageServerInstance

include("exception_types.jl")
include("jsonrpcendpoint.jl")
include("uri2.jl")
include("protocol/protocol.jl")
include("document.jl")
include("languageserverinstance.jl")
include("jsonrpc.jl")
include("staticlint.jl")

include("requests/init.jl")
include("requests/misc.jl")
include("requests/textdocument.jl")
include("requests/features.jl")
include("requests/hover.jl")
include("requests/completions.jl")
include("requests/workspace.jl")
include("requests/actions.jl")
include("utilities.jl")

end


text:
# This file was generated, do not modify it. # hide
V = L * Diagonal(sign.(diag(L)))

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
using CUDA
using ITensors
using ITensorsGrad
using ITensorsGPU
using Zygote
using Random: seed!

CUDA.allowscalar(true)

seed!(1234)

device = cu
#device = cpu

i = Index(2, "i")
A = device(randomITensor(i', dag(i)))

function f(A)
  return (dag(A) * A)[]
end

@show f(A)
dA = @show f'(A)

max_stars_repo_name: mtfishman/ITensorsGrad

max_issues_repo_name: mtfishman/ITensorsGrad

max_forks_repo_name: mtfishman/ITensorsGrad


text:
module SpectralGM

  function __init__()
  end

  using Revise
  using SuiteSparse
  using DataStructures
  using SparseArrays
  using Random
  using LinearAlgebra
  using Statistics
  using Printf

  using DelimitedFiles

  using Plots

  include("cmg_precondition.jl")
  export cmg_precondition, steiner_group

  include("Calculator.jl")
  export test, greet

end


text:
import gmsh

gmsh.initialize(append!(["gmsh"], ARGS))

gmsh.model.add("t10")

lc = .15
gmsh.model.geo.addPoint(0.0, 0.0, 0, lc, 1)
gmsh.model.geo.addPoint(1, 0.0, 0, lc, 2)
gmsh.model.geo.addPoint(1, 1, 0, lc, 3)
gmsh.model.geo.addPoint(0, 1, 0, lc, 4)
gmsh.model.geo.addPoint(0.2, .5, 0, lc, 5)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 5)
gmsh.model.geo.addPlaneSurface([5], 6)

gmsh.model.geo.synchronize()

gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "NodesList", [5])
gmsh.model.mesh.field.setNumber(1, "NNodesByEdge", 100)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [2])

gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", lc / 30)
gmsh.model.mesh.field.setNumber(2, "LcMax", lc)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.15)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)

gmsh.model.mesh.field.add("MathEval", 3)
gmsh.model.mesh.field.setString(3, "F",
                                "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101")

gmsh.model.mesh.field.add("Distance", 4)
gmsh.model.mesh.field.setNumbers(4, "NodesList", [1])

gmsh.model.mesh.field.add("MathEval", 5)
gmsh.model.mesh.field.setString(5, "F", "F4^3 + " * string(lc / 100))

gmsh.model.mesh.field.add("Box", 6)
gmsh.model.mesh.field.setNumber(6, "VIn", lc / 15)
gmsh.model.mesh.field.setNumber(6, "VOut", lc)
gmsh.model.mesh.field.setNumber(6, "XMin", 0.3)
gmsh.model.mesh.field.setNumber(6, "XMax", 0.6)
gmsh.model.mesh.field.setNumber(6, "YMin", 0.3)
gmsh.model.mesh.field.setNumber(6, "YMax", 0.6)

gmsh.model.mesh.field.add("Min", 7)
gmsh.model.mesh.field.setNumbers(7, "FieldsList", [2, 3, 5, 6])

gmsh.model.mesh.field.setAsBackgroundMesh(7)

function meshSizeCallback(dim, tag, x, y, z)
    return 0.02 * x + 0.01
end
gmsh.model.mesh.setSizeCallback(meshSizeCallback)

gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)

gmsh.model.mesh.generate(2)
gmsh.write("t10.msh")

gmsh.finalize()

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using BinaryBuilder

# Collection of sources required to build Ogg
sources = [
    "https://github.com/JuliaInterop/libcxxwrap-julia/archive/v0.2.3.tar.gz" =>
    "705261b9c3e73a501284fc465d15b85bc17f2864cd7c261d5267de7350f9bce6",
]

# Bash recipe for building across all platforms
script = raw"""
# Download julia
cd /usr/local
curl -L 'https://julialang-s3.julialang.org/bin/linux/x64/0.6/julia-0.6.3-linux-x86_64.tar.gz' | tar -zx --strip-components=1

# Build libcxxwrap
cd $WORKSPACE/srcdir/libcxxwrap-julia*
mkdir build && cd build
cmake -DJulia_EXECUTABLE=$(which julia) -DCMAKE_INSTALL_PREFIX=${prefix} ..
VERBOSE=ON cmake --build . --config Debug --target install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, :glibc),
]

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libcxxwrap", :libcxxwrap),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "Cxxwrap", sources, script, platforms, products, dependencies)

max_stars_repo_name: staticfloat/CxxwrapBuilder

max_issues_repo_name: staticfloat/CxxwrapBuilder

max_forks_repo_name: staticfloat/CxxwrapBuilder


text:
include("ppo_trajectory.jl")

using Random
using Distributions: Categorical, Normal, logpdf
using StructArrays

export PPOPolicy

"""
    PPOPolicy(;kwargs)

# Keyword arguments

- `approximator`,
- `γ = 0.99f0`,
- `λ = 0.95f0`,
- `clip_range = 0.2f0`,
- `max_grad_norm = 0.5f0`,
- `n_microbatches = 4`,
- `n_epochs = 4`,
- `actor_loss_weight = 1.0f0`,
- `critic_loss_weight = 0.5f0`,
- `entropy_loss_weight = 0.01f0`,
- `dist = Categorical`,
- `rng = Random.GLOBAL_RNG`,

By default, `dist` is set to `Categorical`, which means it will only works
on environments of discrete actions. To work with environments of
"""
mutable struct PPOPolicy{A<:ActorCritic,D,R} <: AbstractPolicy
    approximator::A
    γ::Float32
    λ::Float32
    clip_range::Float32
    max_grad_norm::Float32
    n_microbatches::Int
    n_epochs::Int
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    rng::R
    # for logging
    norm::Matrix{Float32}
    actor_loss::Matrix{Float32}
    critic_loss::Matrix{Float32}
    entropy_loss::Matrix{Float32}
    loss::Matrix{Float32}
end

function PPOPolicy(;
    approximator,
    γ = 0.99f0,
    λ = 0.95f0,
    clip_range = 0.2f0,
    max_grad_norm = 0.5f0,
    n_microbatches = 4,
    n_epochs = 4,
    actor_loss_weight = 1.0f0,
    critic_loss_weight = 0.5f0,
    entropy_loss_weight = 0.01f0,
    dist = Categorical,
    rng = Random.GLOBAL_RNG,
)
    PPOPolicy{typeof(approximator),dist,typeof(rng)}(
        approximator,
        γ,
        λ,
        clip_range,
        max_grad_norm,
        n_microbatches,
        n_epochs,
        actor_loss_weight,
        critic_loss_weight,
        entropy_loss_weight,
        rng,
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
        zeros(Float32, n_microbatches, n_epochs),
    )
end

function RLBase.get_prob(
    p::PPOPolicy{<:ActorCritic{<:NeuralNetworkApproximator{<:GaussianNetwork}},Normal},
    state::AbstractArray,
)
    p.approximator.actor(send_to_device(device(p.approximator), state)) |>
    send_to_host |>
    StructArray{Normal}
end

function RLBase.get_prob(p::PPOPolicy{<:ActorCritic,Categorical}, state::AbstractArray)
    logits =
        p.approximator.actor(send_to_device(device(p.approximator), state)) |>
        softmax |>
        send_to_host
    [Categorical(x; check_args = false) for x in eachcol(logits)]
end

RLBase.get_prob(p::PPOPolicy, env::MultiThreadEnv) = get_prob(p, get_state(env))

function RLBase.get_prob(p::PPOPolicy, env::AbstractEnv)
    s = get_state(env)
    s = Flux.unsqueeze(s, ndims(s) + 1)
    get_prob(p, s)[1]
end

(p::PPOPolicy)(env::MultiThreadEnv) = rand.(p.rng, get_prob(p, env))
(p::PPOPolicy)(env::AbstractEnv) = rand(p.rng, get_prob(p, env))

function RLBase.update!(p::PPOPolicy, t::PPOTrajectory)
    isfull(t) || return

    states = t[:state]
    actions = t[:action]
    action_log_probs = t[:action_log_prob]
    rewards = t[:reward]
    terminals = t[:terminal]
    states_plus = t[:full_state]

    rng = p.rng
    AC = p.approximator
    γ = p.γ
    λ = p.λ
    n_epochs = p.n_epochs
    n_microbatches = p.n_microbatches
    clip_range = p.clip_range
    w₁ = p.actor_loss_weight
    w₂ = p.critic_loss_weight
    w₃ = p.entropy_loss_weight
    D = device(AC)

    n_envs, n_rollout = size(terminals)
    @assert n_envs * n_rollout % n_microbatches == 0 "size mismatch"
    microbatch_size = n_envs * n_rollout ÷ n_microbatches

    states_flatten = flatten_batch(states)
    states_plus_flatten = flatten_batch(states_plus)
    states_plus_values =
        reshape(send_to_host(AC.critic(send_to_device(D, states_plus_flatten))), n_envs, :)
    advantages =
        generalized_advantage_estimation(rewards, states_plus_values, γ, λ; dims = 2)
    returns = advantages .+ select_last_dim(states_plus_values, 1:n_rollout)

    # TODO: normalize advantage
    for epoch in 1:n_epochs
        rand_inds = shuffle!(rng, Vector(1:n_envs*n_rollout))
        for i in 1:n_microbatches
            inds = rand_inds[(i-1)*microbatch_size+1:i*microbatch_size]
            s = send_to_device(D, select_last_dim(states_flatten, inds))
            if haskey(t, :legal_actions_mask)
                lam = send_to_device(
                    D,
                    select_last_dim(flatten_batch(t[:legal_actions_mask]), inds),
                )
            end
            a = vec(actions)[inds]
            r = send_to_device(D, vec(returns)[inds])
            log_p = send_to_device(D, vec(action_log_probs)[inds])
            adv = send_to_device(D, vec(advantages)[inds])

            ps = Flux.params(AC)
            gs = gradient(ps) do
                v′ = AC.critic(s) |> vec
                if AC.actor isa NeuralNetworkApproximator{<:GaussianNetwork}
                    μ, σ = AC.actor(s)
                    log_p′ₐ = normlogpdf(μ, σ, a)
                    entropy_loss = mean((log(2.0f0π) + 1) / 2 .+ log.(σ))
                else
                    # actor is assumed to return discrete logits
                    logit′ = AC.actor(s)
                    p′ = softmax(logit′)
                    log_p′ = logsoftmax(logit′)
                    log_p′ₐ = log_p′[CartesianIndex.(a, 1:length(a))]
                    entropy_loss = -sum(p′ .* log_p′) * 1 // size(p′, 2)
                end

                ratio = exp.(log_p′ₐ .- log_p)
                surr1 = ratio .* adv
                surr2 = clamp.(ratio, 1.0f0 - clip_range, 1.0f0 + clip_range) .* adv

                actor_loss = -mean(min.(surr1, surr2))
                critic_loss = mean((r .- v′) .^ 2)
                loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss

                ignore() do
                    p.actor_loss[i, epoch] = actor_loss
                    p.critic_loss[i, epoch] = critic_loss
                    p.entropy_loss[i, epoch] = entropy_loss
                    p.loss[i, epoch] = loss
                end

                loss
            end

            p.norm[i, epoch] = clip_by_global_norm!(gs, ps, p.max_grad_norm)
            update!(AC, gs)
        end
    end
end

function (agent::Agent{<:Union{PPOPolicy,RandomStartPolicy{<:PPOPolicy}}})(
    ::Training{PreActStage},
    env::MultiThreadEnv,
)
    state = get_state(env)
    dist = get_prob(agent.policy, env)

    # currently RandomPolicy returns a Matrix instead of a (vector of) distribution.
    if dist isa Matrix{<:Number}
        dist = [Categorical(x; check_args = false) for x in eachcol(dist)]
    elseif dist isa Vector{<:Vector{<:Number}}
        dist = [Categorical(x; check_args = false) for x in dist]
    end

    # !!! a little ugly
    rng = if agent.policy isa PPOPolicy
        agent.policy.rng
    elseif agent.policy isa RandomStartPolicy
        agent.policy.policy.rng
    end

    action = [rand(rng, d) for d in dist]
    action_log_prob = [logpdf(d, a) for (d, a) in zip(dist, action)]
    push!(
        agent.trajectory;
        state = state,
        action = action,
        action_log_prob = action_log_prob,
    )
    update!(agent.policy, agent.trajectory)

    # the main difference is we'd like to flush the buffer after each update!
    if isfull(agent.trajectory)
        empty!(agent.trajectory)
        push!(
            agent.trajectory;
            state = state,
            action = action,
            action_log_prob = action_log_prob,
        )
    end

    action
end


text:
module TweetsChannel

using Channels, App
@dependencies

function subscribe()
  haskey(@params(:payload), "search_term") &&
    ! isempty(@params(:payload)["search_term"] |> strip) &&
      Channels.subscribe(wsclient(@params), "tweets-" * @params(:payload)["search_term"])
  "OK"
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
function serve(dispatcher, io::IO, serializer)
    @unpack deserialize, serialize = Serializer(serializer)
    @debug "Starting RPC connection"
    while isopen(io)
        request = deserialize(io)
        @debug "Received:" request
        request === nothing && break
        jsonrpc(dispatcher, serialize, io, request)
    end
    return io
end

function serve(dispatcher, server::IOServer, serializer)
    tasks = IdDict()
    try
        while isopen(server)
            @debug "Waiting for a client"
            sock = try
                accept(server)
            catch err
                @debug "accept(server)" exception=(err, catch_backtrace())
                if err isa IOError
                    break
                end
                rethrow()
            end
            @debug "Accepted a client"
            @async try
                tasks[current_task()] = sock
                serve(dispatcher, sock, serializer)
                pop!(tasks, current_task())
            catch err
                close(server)
                @debug "`serve` failed" exception=(err, catch_backtrace())
                rethrow()
            end
        end
    finally
        let tasks = copy(tasks)
            foreach(close, values(tasks))
            @sync for t in keys(tasks)
                @async wait(t)
            end
        end
    end
end

serve(dispatcher, listenable, serializer) =
    serve(dispatcher, listen(listenable), serializer)


text:
export
    q2m,
    qxq,
    qdq2av

"""
    q2m(q)

Find the rotation matrix corresponding to a specified unit quaternion.

### Arguments ###

- `q`: A unit quaternion

### Output ###

A rotation matrix corresponding to `q`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/q2m_c.html)
"""
function q2m(q)
    @checkdims 4 q
    r = Matrix{SpiceDouble}(undef, 3, 3)
    ccall((:q2m_c, libcspice), Cvoid, (Ref{SpiceDouble}, Ref{SpiceDouble}), q, r)
    permutedims(r)
end

"""
    qxq(q1, q2)

Multiply two quaternions.

### Arguments ###

- `q1`: First SPICE quaternion factor (as any kind of iterable with four elements)
- `q2`: Second SPICE quaternion factor (as any kind of iterable with four elements)

### Output ###

A quaternion corresponding to the product of `q1' and `q2'

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/qxq_c.html)
"""
function qxq(q1, q2)
    @checkdims 4 q1 q2
    q = Array{SpiceDouble}(undef, 4)
    ccall((:qxq_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          q1, q2, q)
    q
end

"""
    qdq2av(q, dq)

Derive angular velocity from a unit quaternion and its derivative
with respect to time.

### Arguments ###

- `q`: Unit SPICE quaternion (as any kind of iterable with four elements)
- `dq`: Derivative of `q' with respect to time

### Output ###

Angular velocity vector defined by `q' and `dq'

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/qdq2av_c.html)
"""
function qdq2av(q,dq)
    @checkdims 4 q dq
    av = Array{SpiceDouble}(undef, 3)
    ccall((:qdq2av_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          q, dq, av)
    av
end


text:
#=
NacaThicknessFunction.jl
Jacob Child
February 3, 2022
Description:
NACA 4-series Thickness Formula The NACA 4-series formula for airfoil thickness
along the chord is given by the following polynomial function, where m is the value of maximum
thickness (as a percentage of the chord), and x is the x-position along the chord. (This formulation
is for airfoils with sharp trailing edges.)
t = 10m(.2969sqrt(x) - .1620x - .3537x^2 + .2843x^3 - .1015x^4)

Pseudo code- Have a function to ask for and generate x values in an array
Then have a function to calculate various thicknesses
Format as table? and output to screen
***use docstrings***
=#

"""
xinput(x, step)
Takes in the desired length, l, and step size and outputs an x array

"""
function x_input(L, step)
    x = [0:step:L;] #makes array stopping at L with step size step
    return x
end

"""
thickness_calculations(x)
Takes in the x coordinates, calculates thickness, and ouputs the array

"""
function thickness_calculations(x)
    m = 0.10
    t = 10*m.*(.2969*sqrt.(x) - .1260*x - .3537*x.^2 + .2843*x.^3 - .1015*x.^4)
    return t
end

#using Pkg  #Just for the first run through
#Pkg.add("DataFrames")  #just for the first run through

println("Hello! We are going to calculate the thickness along an airfoil \n")
println("\n Enter the Chord length")
L = readline()
L = parse(Int64, L)

println("\n Enter the desired step size to calculate the thickness at")
step = readline()
step = parse(Float64, step)

#println("Type of Inputs", typeof(L), typeof(step))

x = x_input(L, step)
t = thickness_calculations(x)

println(t) #If Running in VS Code Uncomment this as I couldn't get DataFrame to work

#using DataFrames #Only seems to work in Julia Terminal, not VS Code Terminal? ***research

#DataFrame()

#DataFrame("X" => x, "Thickness" => t)


text:
function _export_lp_if_allowed(rd::RoutingData, m::Model, fn::AbstractString)
    if rd.export_lps
        write_LP("$(rd.output_folder)/$(fn).lp", m)
    end
end

function _export_lp_if_failed(rd::RoutingData, status::MOI.TerminationStatusCode, m::Model, fn::AbstractString, errmsg::AbstractString)
    if status != MOI.OPTIMAL
        if rd.export_lps_on_error
            try
                write_LP("$(rd.output_folder)/$(fn).lp", m)
            catch
                write_MOF("$(rd.output_folder)/$(fn).json", m)
            end
        end

        rd.logmessage("$(errmsg) Error code: $(status)")
        error("$(errmsg) Error code: $(status)")
    end
end


text:
"""
A collection of time marching algorithms for vortex modeling

## Exports

- Forward Euler (`forward_euler!`)
- 4th-order Runge-Kutta (`rk4!`)

All schemes have the following function signature

    scheme_name!(x₊, x, t, Δt, compute_ẋ!, update_x!, k)

where

- `x₊`: where the next state is stored
- `x`: current state
- `t`: current time
- `Δt`: time step
- `compute_ẋ!(k, x, t)`: a function with computes the velocity of the state `x` at time `t` and stores the result in `k`
- `update_x!(x₊, x, Δt, k...)`: a function that updates the current state `x` with velocity `k` over the time interval `Δt` and stores the result in `x₊`
- `k...`: one or more temporary variables to store the velocity (see the docstrings for specific schemes)
"""
module TimeMarching

export forward_euler!, rk4!

"""
    forward_euler!(x₊, x, t, Δt, compute_ẋ!, update_x!, k)

Update the state `x` at time `t` using forward Euler

Expects only a single `k`, representing the velocity at time `t`.
"""
function forward_euler!(x₊, x, t, Δt, compute_ẋ!, update_x!, k)
    compute_ẋ!(k, x, t)
    update_x!(x₊, x, k, Δt)
    return x₊, Δt
end

"""
    rk4!(x₊, x, t, Δt, compute_ẋ!, update_x!, k)

Update the state `x` at time `t` using 4th-order Runge-Kutta

Expects `k` to be indexable from 1 to 4, representing velocities at the Runge-Kutta half-steps.
"""
function rk4!(x₊, x, t, Δt, compute_ẋ!, update_x!, k)
    compute_ẋ!(k[1], x, t)

    update_x!(x₊, x, k[1], 0.5Δt)
    compute_ẋ!(k[2], x₊, t + 0.5Δt)

    update_x!(x₊, x, k[2], 0.5Δt)
    compute_ẋ!(k[3], x₊, t + 0.5Δt)

    update_x!(x₊, x, k[3], Δt)
    compute_ẋ!(k[4], x₊, t + Δt)

    update_x!(x₊, x,  k[1], Δt/6)
    update_x!(x₊, x₊, k[2], Δt/3)
    update_x!(x₊, x₊, k[3], Δt/3)
    update_x!(x₊, x₊, k[4], Δt/6)

    return x₊, Δt
end

end


text:
#First add depencies for the example
using Pkg; Pkg.add.(["Plots", "DSP"])
using Weave

list_out_formats()

weave("FIR_design.md", doctype=:"github", out_path=:pwd)

max_stars_repo_path: Deep Into to Julia/00_Config/build.jl

max_stars_repo_name: bmoretz/Statistical-Computing

max_issues_repo_path: Deep Into to Julia/00_Config/build.jl

max_issues_repo_name: bmoretz/Statistical-Computing

max_forks_repo_path: Deep Into to Julia/00_Config/build.jl

max_forks_repo_name: bmoretz/Statistical-Computing


text:
"""
    ConfidenceInterval

An estimated interval around an unknown quantity at a given confidence level
"""
struct ConfidenceInterval{T <: AbstractFloat}
    lower::T
    upper::T
    confidence::T

    function ConfidenceInterval(lower::AbstractFloat, upper::AbstractFloat, confidence::AbstractFloat)
        T = promote_type(typeof(lower), typeof(upper), typeof(confidence), Float32)
        lower <= upper || error("Inconsistent confidence interval: $lower ≰ $upper")
        return new{T}(lower, upper, confidence)
    end
end
lower(interval::ConfidenceInterval) = interval.lower
upper(interval::ConfidenceInterval) = interval.upper
confidence(interval::ConfidenceInterval) = interval.confidence
length(interval::ConfidenceInterval) = upper(interval) - lower(interval)
function in(x::T, interval::ConfidenceInterval) where T <: Real
    isinf(x) && return false
    return interval.lower <= x <= interval.upper
end
function issubset(a::ConfidenceInterval, b::ConfidenceInterval)
    return a.lower >= b.lower && a.upper <= b.upper
end
function show(io::IO, interval::ConfidenceInterval)
    intervalstr = @sprintf "Confidence interval (p = %.0f%%): [%.2f − %.2f]" 100*interval.confidence interval.lower interval.upper
    print(io, intervalstr)
    return io
end


text:
# CUFFT API function return values
typealias cufftResult Cint
const CUFFT_SUCCESS        = 0x0
const CUFFT_INVALID_PLAN   = 0x1
const CUFFT_ALLOC_FAILED   = 0x2
const CUFFT_INVALID_TYPE   = 0x3
const CUFFT_INVALID_VALUE  = 0x4
const CUFFT_INTERNAL_ERROR = 0x5
const CUFFT_EXEC_FAILED    = 0x6
const CUFFT_SETUP_FAILED   = 0x7
const CUFFT_INVALID_SIZE   = 0x8
const CUFFT_UNALIGNED_DATA = 0x9

typealias cufftReal Float32
typealias cufftDoubleReal Float64

typealias cufftComplex Complex64
typealias cufftDoubleComplex Complex128

# CUFFT transform directions
const CUFFT_FORWARD = -1 # Forward FFT
const CUFFT_INVERSE =  1 # Inverse FFT

# CUFFT supports the following transform types
typealias cufftType Cint
const CUFFT_R2C = 0x2a     # Real to Complex
const CUFFT_C2R = 0x2c     # Complex to Real
const CUFFT_C2C = 0x29     # Complex to Complex
const CUFFT_D2Z = 0x6a     # Double to Double-Complex
const CUFFT_Z2D = 0x6c     # Double-Complex to Double
const CUFFT_Z2Z = 0x69     # Double-Complex to Double-Complex

typealias cufftCompatibility Cint
const   CUFFT_COMPATIBILITY_NATIVE          = 0x00
const   CUFFT_COMPATIBILITY_FFTW_PADDING    = 0x01
const   CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC = 0x02
const   CUFFT_COMPATIBILITY_FFTW_ALL        = 0x03

typealias cufftHandle Cint


text:
module DoubleExponentialFormulas

using LinearAlgebra: norm

export
    QuadTS,
    QuadES,
    QuadSS,
    QuadDE,
    quadde,
    quaddeo

"""
    startindex(f::Function, weights, istart::Integer)

Returns the first index i which `f(weights[i][1])` is not NaN,
where `i >= istart`.

This function is employed to avoid sampling `f(x)` with too large `abs(x)` in
trapezoidal rule. For example, the limit of `x*exp(-x^2)` may result in NaN
when `x` approaches infinity since the expression reduced to be `Inf*0` in
float point number computations.
"""
function startindex(f::Function, weights, istart::Integer)
    iend = length(weights)
    for i in istart:iend
        x, _ = @inbounds weights[i]
        y = try
            f(x)
        catch
            continue
        end
        all(.!isnan.(y)) && return i
    end
    return iend
end

include("mapsum.jl")

"""
    estimate_error(T::Type{<:AbstractFloat}, prevI, I)

Estimate an error from the true integral(, though it may not be very accurate).

With double exponential formula, the error from the true value `I` decays
exponentially.
    ΔI(h) = I - I(h) ≈ exp(-C/h)
`I(h)` is the numerical integral with a step size `h`, `C` is a
constant depending on the integrand. If the step size is halved,
    ΔI(h/2) = I - I(h/2) ≈ exp(-2C/h) ≈ ( ΔI(h) )²,
the significant digits gets almost twice. Therefore, `I(h/2)` should be
considerably closer to the true value than `I(h)`.
    ΔI(h/2) ≈ ( ΔI(h) )² ≈ ( I - I(h) )² ≈ ( I(h/2) - I(h) )²
Introduce a magnification coefficient `M` for safety.
    ΔI(h/2) ≈ M*( I(h/2) - I(h) )²
Furthermore, the floating point numbers generally have a finite
significant digits; `ΔI(h/2)` would asymptotically get closer to `I*ε`
at minimum, rather than 0.
    ΔI(h/2) ≈ M*( I(h/2) - I(h) )² + I*ε
            ≈ M*( I(h/2) - I(h) )² + I(h/2)*ε

FIXME: The logical ground of the safety factor is very weak. Any better way?
"""
function estimate_error(T::Type{<:AbstractFloat}, prevI, I)
    ε = eps(T)
    M = 20
    return M*norm(I - prevI)^2 + norm(I)*ε
end

include("quadts.jl")
include("quades.jl")
include("quadss.jl")
include("quadde.jl")

include("quaddeo.jl")

end # module


text:
function rosenbrock2d(x)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function rosenbrock(x)
  return( sum( 100*( x[2:end] - x[1:end-1].^2 ).^2 + ( x[1:end-1] - 1 ).^2 ) )
end

facts("bboptimize") do
  context("example 1 from README") do
    best, fitness = bboptimize(rosenbrock2d; search_range = (-5.0, 5.0), dimensions = 2, parameters = {:ShowTrace => false})
    @fact fitness < 0.001 => true
  end

  context("example 2 from README") do
    best, fitness = bboptimize(rosenbrock2d; search_range = [(-5.0, 5.0), (-2.0, 2.0)], parameters = {:ShowTrace => false})
    @fact fitness < 0.001 => true
  end

  context("example 3 from README") do
    best, fitness = bboptimize(rosenbrock2d; search_range = (-5.0, 5.0), dimensions = 2, method = :de_rand_1_bin, parameters = {:ShowTrace => false})
    @fact fitness < 0.001 => true
  end

  context("example 4 from README") do
    best, fitness = bboptimize(rosenbrock2d; search_range = (-5.0, 5.0), dimensions = 2,
      method = :random_search, max_time = 10.0, parameters = {:ShowTrace => false})
    @fact fitness < 0.2 => true
  end

  context("example 5 from README") do
    BlackBoxOptim.compare_optimizers(rosenbrock; search_range = (-5.0, 5.0), dimensions = 30, max_time = 5.0, parameters = {:ShowTrace => false})
  end

  context("run one longer example in case there is problem with the reporting in long runs") do
    best, fitness = bboptimize(rosenbrock2d; search_range = (-5.0, 5.0), dimensions = 2,
      method = :de_rand_1_bin, parameters = {:ShowTrace => true, :MaxSteps => 25001})
    @fact fitness < 0.001 => true
  end
end


text:
"""
    bootstrap_genie() :: Void

Bootstraps the Genie framework setting up paths and workers. Invoked automatically.
"""
function bootstrap_genie() :: Void
  dirname(@__FILE__) |> cd
  include(joinpath(Pkg.dir("Genie"), "src", "branding.jl"))

  const DEFAULT_NWORKERS_REPL = 1
  const DEFAULT_NWORKERS_SERVER = 1

  isfile("env.jl") && include("env.jl")
  ! haskey(ENV, "GENIE_ENV") && (ENV["GENIE_ENV"] = "dev")
  in("s", ARGS) && ! haskey(ENV, "NWORKERS") ? ENV["NWORKERS"] = DEFAULT_NWORKERS_SERVER : ( haskey(ENV, "NWORKERS") ? ENV["NWORKERS"] : ENV["NWORKERS"] = DEFAULT_NWORKERS_REPL )
  print_with_color(:green, "\nStarting Genie in >> $(ENV["GENIE_ENV"] |> uppercase) << mode using $(ENV["NWORKERS"]) worker(s) \n\n")

  nworkers() < parse(Int, ENV["NWORKERS"]) && addprocs(parse(Int, ENV["NWORKERS"]) - nworkers())

  @everywhere push!(LOAD_PATH,  joinpath("lib"),
                                joinpath(Pkg.dir("Genie"), "src"),
                                joinpath(Pkg.dir("SearchLight"), "src"),
                                joinpath(Pkg.dir("Flax"), "src"),
                                abspath(pwd()))
end

@everywhere bootstrap_genie()
@everywhere import Genie, App
using App
App.@dependencies

@everywhere Genie.run()

try
  using OhMyREPL

  OhMyREPL.input_prompt!("genie>", :magenta)
  OhMyREPL.output_prompt!("genie>", :white)
end


text:
mutable struct Node{N, T<:Real, D<:AbstractData{N,T}}
    center::SVector{N,T}
    length::SVector{N,T}
	p::Union{D, Nothing}           #initilize when an empty leaf becomes nonempty (0->1)
	n::Union{D, Nothing}           #initilize when a nonempty leaf becomes a node (1->2)
	child::Union{Array{Node{N,T,D},1}, Nothing}  #initilize when a nonempty leaf becomes a node (1->2)
	parent::Union{Node{N,T,D}, Nothing}  #only needed for bottom-up treewalk (==nothing for root)
	ID::Int #from 1 to 2^N
end

Node{N,T,D}(center::SVector{N,T}, length::SVector{N,T}) where {N,T,D<:AbstractData{N,T}} =
Node{N,T,D}(center::SVector{N,T}, length::SVector{N,T}, nothing, nothing, nothing, nothing, zero(Int))

function isLeaf(node::Node{N,T,D}) where{N,T,D}
    return node.child == nothing ? true : false
end

function getChildIndex(pos::SVector{N,T}, node::Node{N,T,D}) where {N,T,D<:AbstractData{N,T}}
    idx = 0
    @inbounds for i in 1:N
        if pos[i] > node.center[i]
            idx |= 1<<(i-1)
            #idx |= 2^(i-1)
        end
        #println(i, 2^(i-1), idx)
    end
    return idx+1
end

function getoffset_precomp(N::T) where {T}
	offset_child = zeros(T,N,2^N)
	for i in 1:2^N
	    a = bitstring(i-1)
	    #offset = SVector(2*parse(Int, a[end])-1, 2*parse(Int, a[end-1])-1, 2*parse(Int, a[end-2])-1)
	    offset = T[]
	    @inbounds for j in 1:N
	        push!(offset,2*parse(T, a[end-(j-1)])-1)
	    end
	    #println(offset)
		offset_child[:,i] = offset
	end
	SMatrix{N,2^N,T}(offset_child)
end

const offset = getoffset_precomp(3) #precomputed offset array for 3D

#return the bit of num at jth digit
getbit(num, j) = ( num >> (j-1) ) & 1

#getoffset(i::T, N::T) where {T<:Int} = ([2*parse(T, bitstring(i-1)[end-(j-1)])-1 for j in 1:N])
getoffset(i::T, N::T) where {T<:Int} = getbit.(i-1, (1:N)) .*2 .-1

function insertpart!(p::D, node::Node{N,T,D}) where {N,T,D<:AbstractData{N,T}}
    if isLeaf(node)
        if node.p == nothing
            #println("at an empty leaf, insert particle and we're done")
			@assert node.n == nothing
			node.p = p
        else
            #println("at a leaf with a preexisting particle, so we split it")
			#initializing the child nodes and calculate their centers
			split_node!(node)
			#println("insert back the preexsisting particle...")
			insertpart!(node.p, node.child[getChildIndex(node.p.pos, node)])
			#push!(node.n.idxs, node.p.idx, p.idx)
			assign_node_data!(node.n, node.p, p)
            #println("insert the new particle...")
            insertpart!(p, node.child[getChildIndex(p.pos, node)])
        end
    else
        #println("open node")
		#push!(node.n.idxs, p.idx)
		assign_node_data!(node.n, node.n, p)
        insertpart!(p, node.child[getChildIndex(p.pos,node)])
    end
end

function split_node!(node::Node{N,T,D}) where{N,T,D}
    node.child = Array{Node{N,T,D},1}(undef,2^N)
    @inbounds for i in 1:2^N
		#childcenter = node.center + offset[:,i] * 0.25 .* node.length #precomputed offset. This leads to faster tree build but no flexible dimension
        childcenter = node.center + getoffset(i,N) .* 0.25 .* node.length
        #println("childcenter=", childcenter, "  getoffset(i,N)=", getoffset(i,N))
        node.child[i] = Node{N,T,D}(childcenter, 0.5 .* node.length)
        #println("center = ", childcenter, "  edge_min = ", childcenter - 0.25*node.length,
        #        "  edge_max = ", childcenter + 0.25*node.length)
		node.child[i].parent = node
		node.child[i].ID = i
    end
	#since it's not a leaf node anymore, initialize node data
    node.n = D()
end

#build a tree for given particle locations X
#we distinguish between topnode_length (for the tree) and boxsizes (for PBC wrap)
function buildtree(data::Vector{D}, center::SVector{N,T}, topnode_length::SVector{N,T}) where {N,T,D<:AbstractData{N,T}}

    #the top tree
	tree = Node{N,T,D}(center, topnode_length)

	#insert particles one by one
	length_min = Inf
    @inbounds for i in eachindex(data)
		insertpart!(data[i], tree)
    end
    return tree
end

#PBC wrap: typically used as nearest.(dx) where dx is a separation vector connecting two points
function nearest(x::T, boxsize::T) where {T}
	#can be negative as this is just a vector component
    return (x > 0.5*boxsize) ? (x - boxsize) : ( (x < -0.5*boxsize) ? (x + boxsize) : x )
end

function nearest_plain(x::T) where {T}
    if x < -boxHalf_X
        x = x + boxHalf_X*2
    else
        x = x
    end
    if x > boxHalf_X
        x = x - boxHalf_X*2
    else
        x = x
    end
    return x
end

function get_distance2(a::SVector{N,T}, b::SVector{N,T}, boxsizes::SVector{N,T}, periodic::Bool) where {N,T}
    c = a - b
    if periodic
        c = nearest.(c, boxsizes)
    end
    return sum(c .^ 2)
end


text:
module Khepri
using Base.Iterators
using Sockets
#using PyCall

include("Parameters.jl")
include("Utils.jl")
include("Coords.jl")
include("Paths.jl")
include("Shapes.jl")
include("Geometry.jl")
include("Encoders.jl")
include("Primitives.jl")
include("Tikz.jl")
include("AutoCAD.jl")
include("Rhinoceros.jl")
include("Revit.jl")
include("Grasshopper.jl")
#include("ArchiCAD.jl")
include("Radiance.jl")
include("Unity.jl")
include("Camera.jl")
#const com = PyNULL()

#include("Robot.jl")
end


text:
import GAMS
using JuMP
using Test

@testset "complementarity" begin

   # An MPEC from F. Facchinei, H. Jiang and L. Qi, A smoothing method for
   # mathematical programs with equilibrium constraints, Universita di Roma
   # Technical report, 03.96. Problem number 8
   @testset "gnash1m" begin
      gnash1m = Model(GAMS.Optimizer)
      set_optimizer_attribute(gnash1m, MOI.Silent(), true)

      c = [10, 8, 6, 4, 2]
      K = [5, 5, 5, 5, 5]
      b = [1.2, 1.1, 1.0, 0.9, 0.8]
      L = 20
      g = 1.7

      gg = 5000^(1/g)

      @variable(gnash1m, 0 <= x <= L)
      @variable(gnash1m, 0 <= y[1:4] <= L)
      @variable(gnash1m, l[1:4])
      @variable(gnash1m, Q >= 0)

      @NLobjective(gnash1m, Min, c[1] * x + b[1] / (b[1] + 1) * K[1]^(-1 / b[1]) *
         x^((1 + b[1]) / b[1]) - x * (gg * Q^(-1 / g)))

      @constraint(gnash1m, Q == x + y[1] + y[2] + y[3] + y[4])
      @NLconstraint(gnash1m, [i=1:4], 0 == (c[i+1] + K[i+1]^(-1 / b[i+1]) * y[i]) -
         (gg*Q^(-1 / g)) - y[i] * (-1 / g * gg * Q^(-1 - 1 / g) ) - l[i])
      for i in 1:4
         @constraint(gnash1m, l[i] ⟂ y[i])
      end

      # Initial solutions to help reaching the optimality
      set_start_value.(y, [11, 19, 20, 20])
      set_start_value.(l, [0, 0, -3, -5])
      set_start_value(x, 5)
      set_start_value(Q, 70)

      JuMP.optimize!(gnash1m)

      @test isapprox(JuMP.objective_value(gnash1m), -6.11671, atol=1e-4)
      @test isapprox( JuMP.value.(l)' * (L .- JuMP.value.(y)), 0, atol=1e-5)
   end

   # An MPEC from J.F. Bard, Convex two-level optimization,
   # Mathematical Programming 40(1), 15-27, 1988.
   @testset "bard1" begin
      bard1 = Model(GAMS.Optimizer)
      set_optimizer_attribute(bard1, GAMS.ModelType(), "MPEC")
      set_optimizer_attribute(bard1, MOI.Silent(), true)

      @variable(bard1, x >= 0)
      @variable(bard1, y >= 0)
      @variable(bard1, l[1:3] >= 0)

      @objective(bard1, Min, (x - 5)^2 + (2*y + 1)^2)

      @constraint(bard1, 2*(y-1) - 1.5*x + l[1] - l[2]*0.5 + l[3] == 0)
      @constraint(bard1, 0 <= 3*x - y - 3)
      @constraint(bard1, 0 <= - x + 0.5*y + 4)
      @constraint(bard1, 0 <= - x - y + 7)
      @constraint(bard1, 3*x - y - 3 ⟂ l[1])
      @constraint(bard1, - x + 0.5*y + 4 ⟂ l[2])
      @constraint(bard1, - x - y + 7 ⟂ l[3])

      JuMP.optimize!(bard1)

      @test isapprox(JuMP.objective_value(bard1), 17.0000, atol=1e-4)
   end

   # An MPEC from S. Dempe, "A necessary and sufficient optimality
   # condition for bilevel programming problems", Optimization 25,
   # pp. 341-354, 1992.
   @testset "dempe" begin
      dempe = Model(GAMS.Optimizer)
      set_optimizer_attribute(dempe, MOI.Silent(), true)

      @variable(dempe, x)
      @variable(dempe, z)
      @variable(dempe, w >= 0)

      @objective(dempe, Min, (x - 3.5)^2 + (z + 4)^2)

      @constraint(dempe, z - 3 + 2 * z * w == 0)
      @constraint(dempe, 0 >= z^2 - x)
      @constraint(dempe,  z^2 - x ⟂ w)

      # Initial solutions to help reaching the optimality
      set_start_value(x, 0)
      set_start_value(z, 0)
      set_start_value(w, 1e7)

      JuMP.optimize!(dempe)

      @test isapprox(JuMP.objective_value(dempe), 28.25, atol=1e-4)

      xx = JuMP.value(x)
      zz = JuMP.value(z)
      ww = JuMP.value(w)
      @test isapprox(zz - 3 + 2 * zz * ww, 0, atol=1e-4)
   end

   @testset "MCP" begin
      mcp = Model(GAMS.Optimizer)
      set_optimizer_attribute(mcp, GAMS.ModelType(), "MCP")
      set_optimizer_attribute(mcp, MOI.Silent(), true)

      @variable(mcp, x[1:4] >= 0)

      @constraint(mcp, x[1] + 2x[2] - 2x[3] + 4x[4] - 6 ⟂ x[4])
      @constraint(mcp, -x[3] - x[4] + 2 ⟂ x[1])
      @constraint(mcp, x[1] - x[2] + 2x[3] - 2x[4] - 2 ⟂ x[3])
      @constraint(mcp, x[3] - 2x[4] + 2 ⟂ x[2])

      JuMP.optimize!(mcp)
      @test isapprox(value.(x), [2.8, 0.0, 0.8, 1.2])
   end
end

return


text:
# ==============================================================================
# Docstring for functions
# ==============================================================================

"""
    MP

Immutable Julia structure for Max-Plus scalar. Promote a number of type Float64
or Int64 to a number in the tropical semi-ring (max, +) (ℝ ∪ {-∞}, ⨁, ⨂) where ℝ
is the domain of reals, ⨁ is the usual multiplication and ⨂ is the usual
maximum.

# Notes
`MP(3)` is the equivalent to ScicosLab code: `#(3)`

# Examples
```julia-repl
julia> a = MP(3.5)
Max-Plus 3.5

julia> typeof(a)
MP

julia> a = MP(3)
Max-Plus 3

julia> typeof(a)
MP
```
"""
MP

"""
    MP(x::Bool)

Force conversion from Bool to Max-Plus needed for example by the identity matrix
operator `LinearAlgebra.I` else identity matrices are not well created.

# Examples
```julia-repl
julia> MP(true)
Max-Plus 0.0

julia> MP(false)
Max-Plus -Inf

julia> using LinearAlgebra

julia> Matrix{MP}(I, 2, 2)
2×2 Max-Plus dense matrix:
   0.0   -Inf
  -Inf    0.0
```
"""
MP(x::Bool)

"""
    MP(A::Array)

Convert a dense array from classic algebra to a Max-Plus dense array.
Note: if the array alaready contains Max-Plus number then the `MP()`
is not necessary since Max-Plus number contaminate other numbers.

# Examples (constructor)
```julia-repl
julia> A = MP([1.0 -Inf; 0.0 4])
2×2 Max-Plus dense matrix:
  1.0   -Inf
  0.0    4.0

julia> typeof(A)
Matrix{MP} (alias for Array{MP, 2})
```

# Examples (constructor not needed)
```julia-repl
julia> A = [MP(1.0) -Inf; 0.0 4]
2×2 Max-Plus dense matrix:
  1.0   -Inf
  0.0    4.0

julia> typeof(A)
Matrix{MP} (alias for Array{MP, 2})
```
"""
MP(A::Array)

"""
    MP(A::SparseMatrixCSC)

Convert a sparse array from classic algebra to a Max-Plus sparse array. By
default, explicit Max-Plus zeros (`ε`, `mp0`, `MP(-Inf)`) are removed except if
the parameter `keepzeros` is set to `true`.

# Examples
```julia-repl
julia> using SparseArrays

julia> A = MP(sparse([1, 2, 3], [1, 2, 3], [-Inf, 2, 0]))
3×3 Max-Plus sparse matrix with 2 stored entries:
  .   .   .
  .   2   .
  .   .   0

julia> A.nzval
2-element Max-Plus vector:
  2.0
  0.0
```
"""
MP(S::SparseMatrixCSC)

"""
    MP(I::AbstractVector, J::AbstractVector, V::AbstractVector)

Construct a sparse Max-Plus matrix such as S[I[k], J[k]] = V[k]

# Examples
```julia-repl
julia> MP([1; 2; 3], [1; 2; 3], [42; 43; 44])
3×3 Max-Plus sparse matrix with 3 stored entries:
  42    .    .
   .   43    .
   .    .   44
```
"""
MP(I::AbstractVector, J::AbstractVector, V::AbstractVector)

"""
    MP(V::SparseVector)

Convert a sparse vector from classic algebra to a Max-Plus sparse vector. By
default, explicit Max-Plus zeros (`ε`, `mp0`, `MP(-Inf)`) are removed except if
the parameter `keepzeros` is set to `true`.

# Examples
```julia-repl
julia> MP(sparse([1.0, 0.0, 1.0]))
3-element SparseVector{MP{Float64},Int64} with 2 stored entries:
  [1]  =  1.0
  [3]  =  1.0

julia> A = MP(sparse([1; -Inf; 3]))
3-element SparseVector{MP, Int64} with 2 stored entries:
  [1]  =  1.0
  [3]  =  3.0
```
"""
MP(V::SparseVector)

"""
    MP(x::UnitRange)

Create a Max-Plus dense column vector from a given range.

# Examples
```julia-repl
julia> MP(1:3)
3-element Max-Plus Vector:
 1.0
 2.0
 3.0
```
"""
MP(x::UnitRange)

"""
    MP(x::StepRangeLen)

Create a Max-Plus dense column vector from a given range.

# Examples
```julia-repl
julia> MP(1.0:0.5:3.0)
5-element Max-Plus Vector:
  1.0
  1.5
  2.0
  2.5
  3.0
```
"""
MP(x::StepRangeLen)

"""
    zero(::Type{MP})

Create the constant Max-Plus zero (equals to `-∞`, minus infinity) in classic
algebra which is the neutral for the ⨁ operator. See also `mp0` and `ε`.
"""
Base.zero(::Type{MP})

"""
    zero(::MP)

Create the constant Max-Plus zero (equals to `-∞`, minus infinity) in classic
algebra which is the neutral for the ⨁ operator. See also `mp0` and `ε`.
"""
Base.zero(x::MP)

"""
    mpzero()

Create the constant Max-Plus zero (equals to `-∞`, minus infinity) in classic
algebra which is the neutral for the ⨁ operator. See also `mp0` and `ε`.

# Examples
```julia-repl
julia> mpzero()
Max-Plus -Inf
```
"""
mpzero()

"""
    one(::MP)

Create the constant Max-Plus one (equals to `0`, zero) in classic algebra which is
    the neutral for operators ⨁ and ⨂. See also `mp1` and `mpe`.
"""
Base.one(::Type{MP})

"""
    one(::MP)

Create the constant Max-Plus one (equals to `0`, zero) in classic algebra which is
    the neutral for operators ⨁ and ⨂. See also `mp1` and `mpe`.
"""
Base.one(x::MP)

"""
    mpone()

Create the constant Max-Plus one (equals to `0`, zero) in classic algebra which is
the neutral for operators ⨁ and ⨂. See also `mp1` and `mpe`.

# Examples
```julia-repl
julia> mpone()
Max-Plus 0.0
```
"""
mpone()

"""
    mpI

Is the equivalent of `LinearAlgebra.I` but for Max-Plus type. Allow to create
identity Max-Plus matrices.

# Examples
```julia-repl
julia> typeof(mpI)
LinearAlgebra.UniformScaling{MP}

julia> Matrix(mpI, 2, 2)
2×2 Max-Plus dense matrix:
   0.0   -Inf
  -Inf    0.0

julia> using LinearAlgebra

julia> Matrix(mpI, 2, 2) == Matrix{MP}(I, 2, 2)
true
```
"""
mpI

"""
    +(x::MP, y::MP)

Max-Plus operator ⨁. Return the maximum of `x` and `y` as Max-Plus type. At
least one parameter shall be a Max-Plus number (its conversion to Max-Plus is
automatic).

# Examples
```julia-repl
julia> MP(1.0) + MP(3.0)
Max-Plus 3.0

julia> MP(1.0) + 3
Max-Plus 3.0

julia> 1 + MP(3.0)
Max-Plus 3.0

julia> MP(3.0) + -Inf
Max-Plus 3.0
```
"""
Base.:(+)(x::MP, y::MP)

"""
    *(x::MP, y::MP)

Max-Plus operator ⨂. Return the sum of `x` and `y` as Max-Plus type. At least
one parameter shall be a Max-Plus number (its conversion to Max-Plus is
automatic).

# Examples
```julia-repl
julia> MP(1.0) * MP(3.0)
Max-Plus 4.0

julia> MP(1.0) * 3
Max-Plus 4.0

julia> 1 * MP(3.0)
Max-Plus 4.0

julia> MP(1.0) * -Inf
Max-Plus -Inf
```
"""
Base.:(*)(x::MP, y::MP)

"""
    /(x::MP, y::MP)

Divisor operator. Return the difference between `x` and `y` in classic algebra.

x ⨸ y = x ⨂ y^-1

# Examples
```julia-repl
julia> MP(1.0) / MP(2.0)
Max-Plus -1.0
```
"""
Base.:(/)(x::MP, y::MP)

"""
    -(x::MP, y::MP)

Minus operator (not used for Max-Plus but for Min-Plus).
Return the difference between `x` and `y`.

# Examples
```julia-repl
julia> MP(1.0) - MP(2.0)
Max-Plus -1.0

julia> MP(1.0) / MP(0.0)
Max-Plus 1.0

julia> MP(1.0) / mp0
Max-Plus Inf
```
"""
Base.:(-)(x::MP, y::MP)

"""
    inv(A::ArrMP)

Return the inverse of a scalar or a Max-Plus matrix which is `transpose(-A)` if
the matrix can be inversed. Throw an error if the matrix cannot be inversed.

# Example 1: Scalar
```julia-repl
julia> inv(MP(5))
Max-Plus -5.0

julia> MP(5)^-1
Max-Plus -5.0
```

# Example 2: Square matrix inversible
```julia-repl
julia> A = [mp0 1 mp0; 2 mp0 mp0; mp0 mp0 3]
3×3 Max-Plus dense matrix:
  .   1   .
  2   .   .
  .   .   3

julia> A^-1
3×3 Max-Plus transpose(dense matrix):
   .   -2    .
  -1    .    .
   .    .   -3

julia> inv(A)
3×3 Max-Plus transpose(dense matrix):
   .   -2    .
  -1    .    .
   .    .   -3

julia> A * inv(A)
3×3 Max-Plus dense matrix:
  0   .   .
  .   0   .
  .   .   0

julia> A * inv(A) == inv(A) * A == mpeye(A)
true
```

# Example 3: Square matrix inversible
```julia-repl
julia> A = [mp0 1 mp0; 2 mp0 mp0]
2×3 Max-Plus dense matrix:
  .   1   .
  2   .   .

julia> A^-1
3×2 Max-Plus transposed dense matrix:
   .   -2
  -1    .
   .    .

julia> inv(A)
3×2 Max-Plus transposed dense matrix:
   .   -2
  -1    .
   .    .

julia> A * inv(A)
2×2 Max-Plus dense matrix:
  0   .
  .   0

julia> inv(A) * A
3×3 Max-Plus dense matrix:
  0   .   .
  .   0   .
  .   .   .

julia> (A * inv(A)) != (inv(A) * A)
true
```

# Example 4: Non inversible
```julia-repl
julia> A = MP([1 2; 3 4])
2×2 Max-Plus dense matrix:
  1.0   2.0
  3.0   4.0

julia> inv(A)
ERROR: The matrix cannot be inversed
```
"""
Base.inv(A::ArrMP)

"""
    \\(A::ArrMP, b::ArrMP)

`x = A \\ b` is a solution to `A ⨂ x = b` and its simply computed as
`inv(A) * b`.

# Examples
```julia-repl
julia> A = [mp0 1 mp0; 2 mp0 mp0; mp0 mp0 3]
3×3 Max-Plus dense matrix:
  .   1   .
  2   .   .
  .   .   3

julia> B = [3 mp0 mp0; mp0 mp0 4; mp0 5 mp0]
3×3 Max-Plus dense matrix:
  3   .   .
  .   .   4
  .   5   .

julia> x = A \\ B
3×3 Max-Plus dense matrix:
  .   .   2
  2   .   .
  .   2   .

julia> A * x == B
true

julia> A \\ A
3×3 Max-Plus dense matrix:
  0   .   .
  .   0   .
  .   .   0
```
"""
Base.:(\)(A::ArrMP, b::ArrMP)

"""
    mp0

Create the constant Max-Plus zero (equals to `-∞`, minus infinity) in classic
algebra which is the neutral for the ⨁ operator. See also `mpzero`.

Equivalent to ScicosLab code: `%0` sugar notation for `#(-%inf)`

# Examples
```julia-repl
julia> mp0
Max-Plus -Inf

julia> mp0 * 5
Max-Plus -Inf

julia> mp0 + 5
Max-Plus 5.0
```
"""
mp0

"""
    ε (\\varepsilon)

Create the constant Max-Plus zero (equals to `-∞`, minus infinity) in classic
algebra which is the neutral for the ⨁ operator. See also `mpzero`.

Equivalent to ScicosLab code: `%0` sugar notation for `#(-%inf)`

# Examples
```julia-repl
julia> ε
Max-Plus -Inf

julia> ε * 5
Max-Plus -Inf

julia> ε + 5
Max-Plus 5.0
```
"""
ε

"""
    mp1

Create the constant Max-Plus one (equals to `0`, zero) in classic algebra which is
the neutral for operators ⨁ and ⨂. See also `mpone`.

Equivalent to ScicosLab code: `%1` sugar notation for `#(1)`

# Examples
```julia-repl
julia> mp1
Max-Plus 0.0

julia> mp1 * 5
Max-Plus 5.0

julia> mp1 + 5
Max-Plus 5.0
```
"""
mp1

"""
    mpe

Create the constant Max-Plus one (equals to `0`, zero) in classic algebra which is
the neutral for operators ⨁ and ⨂. See also `mpone`.

Equivalent to ScicosLab code: `%1` sugar notation for `#(1)`

# Examples
```julia-repl
julia> mpe
Max-Plus 0.0

julia> mpe * 5
Max-Plus 5.0

julia> mpe + 5
Max-Plus 5.0
```
"""
mpe

"""
    mptop

Create the constant Min-Plus one (equals to `+∞`, infinity) in classic algebra which is
the neutral for operators ⨁ and ⨂.

Equivalent to ScicosLab code: `%top = #(%inf)`

# Examples
```julia-repl
julia> mptop
Max-Plus Inf

julia> mptop * 5
Max-Plus Inf

julia> mptop + 5
Max-Plus Inf
```
"""
mptop

"""
    mpsparse_map(f, M::SparseMatrixCSC{MP,U})

Map a function ot each element of the Max-Plus sparse matrix.

# Examples
```julia-repl
julia> mpsparse_map(x -> x.λ, sparse(MP([1.0 0; 0 1.0])))
2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0  0.0
 0.0  1.0
```
"""
mpsparse_map(f, M::SpaMP)

"""
    plustimes(x::MP)

Convert a Max-Plus number to a number in standard algebra. An alternative way
could be `x.λ`.

# Examples
```julia-repl
julia> plustimes(MP(1.0))
1.0

julia> typeof(plustimes(MP(1)))
Float64
```
"""
plustimes(n::MP)

"""
    plustimes(A::ArrMP)

Convert a Max-Plus dense matrix to a dense matrix in standard algebra.

# Examples
```julia-repl
julia> A=[MP(1.0) 2.0; ε mpe]
2×2 Max-Plus dense matrix:
   1.0   2.0
  -Inf   0.0

julia> plustimes(A)
2×2 Matrix{Float64}:
   1.0  2.0
  -Inf  0.0
```
"""
plustimes(A::ArrMP)

"""
    plustimes(A::SpaMP)

Convert a Max-Plus sparse matrix to an sparse matrix in standard algebra.

# Examples
```julia-repl
julia> S = sparse(mpeye(2,2))
2×2 Max-Plus sparse matrix with 2 stored entries:
 0.0   ⋅
  ⋅   0.0

julia> findnz(S)
([1, 2], [1, 2], MP[0.0, 0.0])

julia> plustimes(S)
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
 0.0   ⋅
  ⋅   0.0

julia> findnz(plustimes(S))
([1, 2], [1, 2], [0.0, 0.0])
```
"""
plustimes(S::SpaMP)

"""
    full(::SpaMP})

Convert a sparse Max-Plus array to a dense Max-Plus array.
Alternative function name: dense.

# Examples
```julia-repl
julia> full(mpzeros(2,5))
2×5 Max-Plus dense array:
 -Inf  -Inf  -Inf  -Inf  -Inf
 -Inf  -Inf  -Inf  -Inf  -Inf
```
"""
full(S::SpaMP)

"""
    dense(::SpaMP})

Convert a sparse Max-Plus array to a dense Max-Plus array.
Alternative function name: full.

# Examples
```julia-repl
julia> dense(mpzeros(2,5))
2×5 Max-Plus dense array:
 -Inf  -Inf  -Inf  -Inf  -Inf
 -Inf  -Inf  -Inf  -Inf  -Inf
```
"""
dense(S::SpaMP)

"""
    array(::SpaMP})

Convert a sparse Max-Plus array to a dense non Max-Plus array.

# Examples
```julia-repl
julia> array(mpzeros(2,2))
2×2 Matrix{Float64}:
 -Inf  -Inf
 -Inf  -Inf
```
"""
array(S::SpaMP)

"""
    min(x::MP, y::MP)

Return the minimun of Max-Plus numbers `x` and `y`.

min(x, y) = (x ⨂ y) ⨸ (x ⨁ y)

# Examples
```julia-repl
julia> min(MP(1), 3)
Max-Plus 1

julia> min(MP([10 1; 10 1]), MP([4 5; 6 5]))
2×2 Max-Plus dense matrix:
  4   1
  6   1
```
"""
Base.:min(x::MP, y::MP)

"""
    mpeye(n::Int64)

Construct a Max-Plus identity dense n-by-n matrix.

# Examples
```julia-repl
julia> mpeye(2)
2×2 Max-Plus dense matrix:
   0.0   -Inf
  -Inf    0.0
```
"""
mpeye(n::Int64)

"""
    mpeye(m::Int64, n::Int64)

Construct a Max-Plus identity dense m-by-n matrix.

# Examples
```julia-repl
julia> mpeye(2, 3)
2×3 Max-Plus dense matrix:
   0.0   -Inf   -Inf
  -Inf    0.0   -Inf
```
"""
mpeye(m::Int64, n::Int64)

"""
    mpeye(A::Array{T})

Construct a Max-Plus identity dense matrix of same dimension
and of the same type that the matrix given as parameter.

# Examples
```julia-repl
julia> A=[1.0 2; 3 4]
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0

julia> mpeye(A)
2×2 Max-Plus dense matrix:
  0.0  -Inf
 -Inf   0.0

julia> mpeye(MP(A))
2×2 Max-Plus dense matrix:
  0.0  -Inf
 -Inf   0.0
```
"""
mpeye(A::Array)

"""
    mpzeros(n::Int64)

Construct a Max-Plus zero n-by-m sparse matrix.

# Examples
```julia-repl
julia> mpzeros(2)
2-element SparseVector{MP,Int64} with 0 stored entries
```
"""
mpzeros(n::Int64)

"""
    mpzeros(n::Int64)

Construct a sparse Max-Plus zero m-by-n matrix.

# Examples
```julia-repl
julia> mpzeros(2,5)
2×5 SparseMatrixCSC{MP,Int64} with 0 stored entries
  .   .   .   .   .
  .   .   .   .   .

julia> full(mpzeros(2,5))
2×5 Max-Plus dense matrix:
  -Inf   -Inf   -Inf   -Inf   -Inf
  -Inf   -Inf   -Inf   -Inf   -Inf
```
"""
mpzeros(m::Int64, n::Int64)

"""
    mpzeros(A::Array{T})

Construct a sparse  Max-Plus zero matrix of same dimension
and of the same type that the matrix given as parameter.

# Examples
```julia-repl
julia> A=[1.0 2; 3 4]
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0

julia> mpzeros(A)
2×2 Max-Plus sparse matrix with 0 stored entries:
  .   .
  .   .

julia> mpzeros(MP(A))
2×2 Max-Plus sparse matrix with 0 stored entries:
  .   .
  .   .
```
"""
mpzeros(A::Array)

"""
    mpones(n::Int64)

Construct a Max-Plus one n-by-1 matrix.

# Examples
```julia-repl
julia> mpones(2)
2-element Max-Plus vector:
  0.0
  0.0
```
"""
mpones(n::Int64)

"""
    mpones(m::Int64, n::Int64)

Construct a Max-Plus one m-by-n matrix.

# Examples
```julia-repl
julia> mpones(3,2)
3×2 Max-Plus dense matrix:
  0.0   0.0
  0.0   0.0
  0.0   0.0
```
"""
mpones(m::Int64, n::Int64)

"""
    mpones(A::Array{T})

Construct a sparse Max-Plus ones matrix of same dimension
and of the same type that the matrix given as parameter.

# Examples
```julia-repl
julia> A=[1.0 2; 3 4]
2×2 Matrix{Float64}:
 1.0  2.0
 3.0  4.0

julia> mpones(A)
2×2 Max-Plus dense matrix:
  0.0   0.0
  0.0   0.0
```
"""
mpones(A::Array)

"""
    mptrace(A::Array)

Compute the trace of the matrix (summation of diagonal elements).

# Examples
```julia-repl
julia> mptrace([MP(1) 2; 3 4])
Max-Plus 4.0
```
"""
mptrace(A::ArrMP)

"""
    mpnorm(A)

Compute the norm of the full or sparce matrix A.
Return the largest entry minus smallest entry of A.

# Examples
```
julia-repl
julia> A = MP([1 20 2;30 400 4;4 50 10])
3×3 Max-Plus dense matrix:
   1.0    20.0    2.0
  30.0   400.0    4.0
   4.0    50.0   10.0

julia> S = MP(sparse(A))
3×3 Max-Plus sparse matrix with 9 stored entries:
   1    20    2
  30   400    4
   4    50   10

julia> mpnorm(A)
Max-Plus 399

julia> mpnorm(S)
Max-Plus 399
```
"""
mpnorm(A::Array)

"""
    B = mpstar(A::ArrMP)

Solve `x = Ax + I` in the Max-Plus algebra. When there is no circuits with
positive weight in G(A) (the incidence graph of A) `B = I + A + ... + A^(n-1)`
where n denotes the order of the square matrix A.

See also mpplus.

# Arguments
- A : Max-Plus full square matrix.
- B : Id + A + A^2 + ...

# Examples
```julia-repl
julia> mpstar(MP(1.0))
Max-Plus Inf

julia> mpstar(MP(-1.0))
Max-Plus 0

julia> mpstar(MP([1.0 2; 3 4]))
2×2 Max-Plus dense matrix:
  Inf   Inf
  Inf   Inf

julia> A = mpstar(MP([-3.0 -2;-1 0]))
2×2 Max-Plus dense matrix:
   0   -2
  -1    0

julia> B = mpstar(A)
2×2 Max-Plus dense matrix:
   0   -2
  -1    0

julia> B == (mpeye(2,2) + A)
true

julia> B == (B + A * A)
true
```
"""
mpstar(A::ArrMP)

"""
    mpstar(x::MP)

Make x a 1x1 matrix then call mpstar(A::ArrMP).
See mpstar(A::ArrMP) for more information.
"""
mpstar(x::MP)

"""
    mpastarb(A::ArrMP, b::ArrMP)

Max-Plus linear system solution.

Solve `x = Ax + b` in the Max-Plus algebra when there is no circuits with
positive weight in `G(A')` (the incidence graph of `A'`, that is it exists an
arc from `j` to `i` if `A_ij` is nonzero).

TODO It is much more efficient in time and memory than `mpstar(A) * b`.
"""
mpastarb(A::ArrMP, b::ArrMP)

"""
    B = mpplus(A::ArrMP)

Compute `A * A^* = A + A^2 + ...` of a maxplus matrix A.
See also mpstar.

# Arguments
- A : Max-Plus full square matrix.
- B : Id + A + A^2 + ...

# Examples
```julia-repl
julia> mpplus(MP(1.0))
Max-Plus Inf

julia> mpplus(MP(-1.0))
Max-Plus -1

julia> A = MP([-3.0 -2; -1 0]); B = mpplus(A)
2×2 Max-Plus dense matrix:
  -3   -2
  -1    0

julia> B == (A * mpstar(A))
true
```
"""
mpplus(A::ArrMP)

"""
    mpplus(x::MP)

Make x a 1x1 matrix then call mpplus(A::ArrMP).
See mpplus(A::ArrMP) for more information.
"""
mpplus(x::MP)


text:
@enum BoundaryType begin
    Vaccum
    Reflective
    Periodic
end

struct BoundaryConditions{T<:BoundaryType}
    top::T
    bottom::T
    right::T
    left::T
end

BoundaryConditions(; top=Vaccum, bottom=Vaccum, right=Vaccum, left=Vaccum) =
    BoundaryConditions(top, bottom, right, left)

function boundary_condition(x::Point2D, sides, bcs::BoundaryConditions)

    if x in sides.top
        bc = bcs.top
    elseif x in sides.bottom
        bc = bcs.bottom
    elseif x in sides.right
        bc = bcs.right
    elseif x in sides.left
        bc = bcs.left
    else
        error("Point do not lie in the boundary.")
    end

    return bc
end


text:
file = open(OCTHaGOn.BARON_DIR*"instances.txt", "r")
gams_problem_names = []
while !eof(file)
    # read a new / next line for every iteration
    s = readline(file)
    push!(gams_problem_names, string(s) * ".gms")
 end

 for name in gams_problem_names
    try
        gm = GAMS_to_GlobalModel(OCTHaGOn.GAMS_DIR, name)
        @info "$(name) worked."
    catch
        @info "$(name) didn't work."
    end
 end

 name = "flay05m.gms" # GOOD EXAMPLE
 gm = GAMS_to_GlobalModel(OCTHaGOn.GAMS_DIR, name)
 set_optimizer(gm, SOLVER_SILENT)
 globalsolve_and_time!(gm)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CustomResourceDefinitionVersion describes a version for CRD.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion(;
        additionalPrinterColumns=nothing,
        name=nothing,
        schema=nothing,
        served=nothing,
        storage=nothing,
        subresources=nothing,
    )

    - additionalPrinterColumns::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceColumnDefinition} : additionalPrinterColumns specifies additional columns returned in Table output. See https://kubernetes.io/docs/reference/using-api/api-concepts/#receiving-resources-as-tables for details. Top-level and per-version columns are mutually exclusive. Per-version columns must not all be set to identical values (top-level columns should be used instead). If no top-level or per-version columns are specified, a single column displaying the age of the custom resource is used.
    - name::String : name is the version name, e.g. “v1”, “v2beta1”, etc. The custom resources are served under this version at `/apis/<group>/<version>/...` if `served` is true.
    - schema::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation : schema describes the schema used for validation and pruning of this version of the custom resource. Top-level and per-version schemas are mutually exclusive. Per-version schemas must not all be set to identical values (top-level validation schema should be used instead).
    - served::Bool : served is a flag enabling/disabling this version from being served via REST APIs
    - storage::Bool : storage indicates this version should be used when persisting custom resources to storage. There must be exactly one version with storage=true.
    - subresources::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresources : subresources specify what subresources this version of the defined custom resource have. Top-level and per-version subresources are mutually exclusive. Per-version subresources must not all be set to identical values (top-level subresources should be used instead).
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion <: SwaggerModel
    additionalPrinterColumns::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceColumnDefinition} } # spec name: additionalPrinterColumns
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    schema::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation } # spec name: schema
    served::Any # spec type: Union{ Nothing, Bool } # spec name: served
    storage::Any # spec type: Union{ Nothing, Bool } # spec name: storage
    subresources::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresources } # spec name: subresources

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion(;additionalPrinterColumns=nothing, name=nothing, schema=nothing, served=nothing, storage=nothing, subresources=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        setfield!(o, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("served"), served)
        setfield!(o, Symbol("served"), served)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("storage"), storage)
        setfield!(o, Symbol("storage"), storage)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion, Symbol("subresources"), subresources)
        setfield!(o, Symbol("subresources"), subresources)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion = Dict{Symbol,Symbol}(Symbol("additionalPrinterColumns")=>Symbol("additionalPrinterColumns"), Symbol("name")=>Symbol("name"), Symbol("schema")=>Symbol("schema"), Symbol("served")=>Symbol("served"), Symbol("storage")=>Symbol("storage"), Symbol("subresources")=>Symbol("subresources"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion = Dict{Symbol,String}(Symbol("additionalPrinterColumns")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceColumnDefinition}", Symbol("name")=>"String", Symbol("schema")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceValidation", Symbol("served")=>"Bool", Symbol("storage")=>"Bool", Symbol("subresources")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceSubresources")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("served")) === nothing) && (return false)
    (getproperty(o, Symbol("storage")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceDefinitionVersion }, name::Symbol, val)
end