text:
"""
    mechanical_energy(mechanism, storage)

    mechanism's total mechanical energy

    mechanism: Mechanism
    storage: Storage
"""
function mechanical_energy(mechanism::Mechanism, storage::Storage)
    ke = kinetic_energy(mechanism, storage)
    pe = potential_energy(mechanism, storage)
    me = ke + pe
    return me
end

"""
    kinetic_energy(mechanism, storage)

    mechanism's kinetic energy from linear and angular velocity

    mechanism: Mechanism
    storage: Storage
"""
function kinetic_energy(mechanism::Mechanism, storage::Storage{T,N}) where {T,N}
    ke = zeros(T,N)
    for i = 1:N
        ke[i] = kinetic_energy(mechanism, storage, i)
    end
    return ke
end

function kinetic_energy(mechanism::Mechanism, storage::Storage{T,N}, t::Int) where {T,N}
    ke = 0.0
    for (i,body) in enumerate(mechanism.bodies)
        vl = storage.vl[i][t]
        ωl = storage.ωl[i][t]
        ke += 0.5 * body.mass * vl' * vl
        ke += 0.5 * ωl' * body.inertia * ωl
    end
    return ke
end

"""
    potential_energy(mechanism, storage)

    mechanism's potential energy from gravity and springs

    mechanism: Mechanism
    storage: Storage
"""
function potential_energy(mechanism::Mechanism, storage::Storage{T,N}) where {T,N}
    pe = zeros(T,N)
    for i = 1:N
        pe[i] = potential_energy(mechanism, storage, i)
    end
    return pe
end

function potential_energy(mechanism::Mechanism{T,Nn,Ne,Nb}, storage::Storage{T,Ns}, t::Int) where {T,Nn,Ne,Nb,Ns}
    pe = 0.0

    # Gravity
    for (i,body) in enumerate(mechanism.bodies)
        pe -= body.mass * dot(mechanism.gravity, storage.x[i][t]) # TODO: confirm this is correct
    end

    # Springs
    for joint in mechanism.joints
        parent_id = joint.parent_id
        child_id = joint.child_id
        if joint.spring
            for element in (joint.translational, joint.rotational)
                if element.spring > 0
                    xb = storage.x[child_id - Ne][t] # TODO this is sketchy way to get the correct index
                    qb = storage.q[child_id - Ne][t] # TODO this is sketchy way to get the correct index

                    if parent_id != 0
                        xa = storage.x[parent_id - Ne][t] # TODO this is sketchy way to get the correct index
                        qa = storage.q[parent_id - Ne][t] # TODO this is sketchy way to get the correct index
                    else
                        xa, qa = current_configuration(mechanism.origin.state)
                    end

                    force = spring_force(:parent, element, xa, qa, xb, qb)
                    pe += 0.5 * dot(force, force) ./ element.spring
                end
            end
        end
    end
    return pe
end


text:
function lsodafun(t::T1, y::T2, yp::T3, userfun::UserFunctionAndData) where {T1, T2, T3}
  y_ = unsafe_wrap(Array,y,userfun.neq)
  ydot_ = unsafe_wrap(Array,yp,userfun.neq)
  userfun.func(t, y_, ydot_,userfun.data)
  return Int32(0)
end

function lsoda_0(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdata::Any=nothing, reltol::Union{Float64,Vector}=1e-4, abstol::Union{Float64,Vector}=1e-10)
  neq = Int32(length(y0))
  userfun = UserFunctionAndData(f, userdata,neq)

  atol = ones(Float64,neq)
  rtol = ones(Float64,neq)

  if typeof(abstol) == Float64
    atol *= abstol
  else
    atol = copy(abstol)
  end

  if typeof(reltol) == Float64
    rtol *= reltol
  else
    rtol = copy(reltol)
  end

  t = Float64[tspan[1]]
  tout = Float64[tspan[2]]
  y = copy(y0)

  #
  opt = lsoda_opt_t()
    opt.ixpr = 0
    opt.rtol = pointer(rtol)
    opt.atol = pointer(atol)
    opt.itask = 1
  #

  ctx = lsoda_context_t()
    ctx.function_ = @cfunction($lsodafun,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ref{UserFunctionAndData})).ptr
    ctx.neq = neq
    ctx.state = 1
    ctx.data = pointer_from_objref(userfun)

  lsoda_prepare(ctx,opt)
  for i=1:12
    lsoda(ctx,y,t,tout[1])
	  @assert (ctx.state >0) string("LSODA error istate = ", ctx.state)
    @printf("at t = %12.4e y= %14.6e %14.6e %14.6e\n",t[1],y[1], y[2], y[3])
    tout[1] *= 10.0E0
  end
end

"""
  lsoda(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdata::Any=nothing, reltol::Union{Float64,Vector}=1e-4, abstol::Union{Float64,Vector}=1e-10)

Solves a set of ordinary differential equations using the LSODA algorithm. The vector field encoded in an inplace f::Function needs to have the self-explanatory arguments f(t, y, ydot, data)
"""
function lsoda(f::Function, y0::Vector{Float64}, tspan::Vector{Float64}; userdata::Any=nothing, reltol::Union{Float64,Vector}=1e-4, abstol::Union{Float64,Vector}=1e-10,nbsteps = 10000)
  neq = Int32(length(y0))
  userfun = UserFunctionAndData(f, userdata, neq)

  atol = ones(Float64,neq)
  rtol = ones(Float64,neq)

  yres = zeros(length(tspan), length(y0))

  if typeof(abstol) == Float64
    atol *= abstol
  else
    atol = copy(abstol)
  end

  if typeof(reltol) == Float64
    rtol *= reltol
  else
    rtol = copy(reltol)
  end

  t    = Float64[tspan[1]]
  tout = Float64[tspan[2]]
  y    = copy(y0)

  opt = lsoda_opt_t()
    opt.mxstep = nbsteps
    opt.ixpr = 0
    opt.rtol = pointer(rtol)
    opt.atol = pointer(atol)
    opt.itask = 1

  ctx_ptr = lsoda_context_t()
    ctx_ptr.function_ = @cfunction($lsodafun,Cint,(Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ref{UserFunctionAndData})).ptr
    ctx_ptr.neq = neq
    ctx_ptr.state = 1
    ctx_ptr.data = pointer_from_objref(userfun)

  lsoda_prepare(ctx_ptr,opt)
  yres[1,:] = y0

  for k in 2:length(tspan)
		tout[1] = tspan[k]
	  lsoda(ctx_ptr,y,t,tout[1])
		@assert (ctx_ptr.state >0) string("LSODA error istate = ", ctx_ptr.state, ", error = ",unsafe_string(ctx_ptr.error))
		yres[k,:] = copy(y)
  end

  lsoda_free(ctx_ptr)
  return yres
end

"""
  lsoda_evolve!(ctx::lsoda_context_t,y::Vector{Float64},tspan::Vector{Float64})

Solves a set of ordinary differential equations using the LSODA algorithm and the context variable ctx. This avoid re-allocating ctx. You have to be carefull to remember the current time or this function will return an error.
"""
function lsoda_evolve!(ctx::lsoda_context_t,y::Vector{Float64},tspan::Vector{Float64})
	@assert ctx.neq == length(y)
# if userdata != nothing
# 		# this functionality is not working yet
# 		# ctx.data.data = userdata
# 		# unsafe_pointer_to_objref(ctx.data).data = userdata
# 	end

	t    = Float64[tspan[1]]
	tout = Float64[tspan[2]]
	lsoda(ctx,y,t,tout[1])
	@assert (ctx.state >0) string("LSODA error istate = ", ctx.state)
end

ext: jl


text:
using Documenter,DiffEqDevTools,DiffEqBase,FiniteElementDiffEq,
      DiffEqProblemLibrary, StokesDiffEq, OrdinaryDiffEq

makedocs(modules=[DiffEqDevTools,DiffEqBase,FiniteElementDiffEq,
                  StokesDiffEq,OrdinaryDiffEq,DiffEqProblemLibrary],
         doctest=false, clean=true,
         format =:html,
         sitename="DifferentialEquations.jl",
         authors="Chris Rackauckas",
         pages = Any[
         "Home" => "index.md",
         "Tutorials" => Any[
           "tutorials/ode_example.md",
           "tutorials/sde_example.md",
           "tutorials/dde_example.md",
           "tutorials/dae_example.md",
           "tutorials/discrete_stochastic_example.md",
           "tutorials/jump_diffusion.md",
           "tutorials/fempoisson_example.md",
           "tutorials/femheat_example.md",
           "tutorials/femstochastic_example.md"
         ],
         "Basics" => Any[
           "basics/overview.md",
           "basics/common_solver_opts.md",
           "basics/solution.md",
           "basics/plot.md",
           "basics/integrator.md",
           "basics/compatibility_chart.md"
         ],
         "Equation Types" => Any[
           "types/discrete_types.md",
           "types/ode_types.md",
           "types/sde_types.md",
           "types/dde_types.md",
           "types/dae_types.md",
           "types/jump_types.md",
           "types/fem_types.md",
           "types/stokes_types.md"
         ],
         "Solver Algorithms" => Any[
           "solvers/discrete_solve.md",
           "solvers/ode_solve.md",
           "solvers/sde_solve.md",
           "solvers/dde_solve.md",
           "solvers/dae_solve.md",
           "solvers/fempoisson_solve.md",
           "solvers/femheat_solve.md",
           "solvers/fdmstokes_solve.md"
         ],
         "Additional Features" => Any[
           "features/performance_overloads.md",
           "features/data_arrays.md",
           "features/linear_nonlinear.md",
           "features/callback_functions.md",
           "features/callback_library.md",
           "features/monte_carlo.md",
           "features/low_dep.md",
           "features/mesh.md",
           "features/output_specification.md",
           "features/progress_bar.md"
         ],
         "Analysis Tools" => Any[
           "analysis/parameterized_functions.md",
           "analysis/parameter_estimation.md",
           "analysis/bifurcation.md",
           "analysis/sensitivity.md",
           "analysis/uncertainty_quantification.md",
           "analysis/dev_and_test.md"
         ],
         "Modeling Tools" => Any[
             "models/multiscale.md",
             "models/financial.md",
             "models/biological.md"
         ]
         ])

deploydocs(
   repo = "github.com/JuliaDiffEq/DiffEqDocs.jl.git",
   target = "build",
   osname = "linux",
   julia = "0.5",
   deps = nothing,
   make = nothing)


text:
"""Wrapper for Python calls."""
function pychain(graph::IndexedGraph,
                 plan::Plan,
                 n_steps::Int,
                 reversible::Bool,
                 min_pop::Int,
                 max_pop::Int,
                 twister::MersenneTwister,
                 params::Dict)::Dict
    if reversible
        proposal = reversible_recom_until_step
    else
        proposal = recom_until_step
    end
    steps = Array{Dict}(undef, 0)
    step_count = 0
    initial_assignment = Dict(node - 1 => assignment for (node, assignment)
                              in enumerate(plan.assignment))
    max_balanced_cuts = 0
    proposal_args = Dict(Symbol(k) => v for (k, v) in params)
    while step_count < n_steps
        flip, self_loops, reasons = proposal(graph, plan, min_pop, max_pop, twister;
                                             (; proposal_args...)...)
        max_balanced_cuts = maximum([max_balanced_cuts, flip.balanced_cuts])
        update!(plan, graph, flip)
        # Use Python indexing
        pyflip = Dict(node - 1 => assignment for (node, assignment)
                      in zip(flip.nodes, flip.new_assignments))
        push!(steps, Dict("flip" => pyflip,
                          "self_loops" => self_loops,
                          "reasons" => reasons,
                          "district_adj" => plan.district_adj))
        step_count += (self_loops + 1)
    end
    return Dict("initial_assignment" => initial_assignment,
                "steps" => steps,
                "max_balanced_cuts" => max_balanced_cuts)
end


text:
module AnalogDataAssimilation

include("utils.jl")
include("models.jl")
include("time_series.jl")
include("state_space.jl")
include("catalog.jl")
include("generate_data.jl")
include("data_assimilation.jl")
include("model_forecasting.jl")
include("local_linear.jl")
include("analog_forecasting.jl")

export forecast

include("ensemble_kalman_filters.jl")
include("ensemble_kalman_smoothers.jl")
include("particle_filters.jl")

end


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("Kline.csv"), DataFrame);

# New col log_pop, set log() for population data
df[!, :log_pop] = map((x) -> log(x), df[!, :population]);
df[!, :society] = 1:10;

stan12_6 = "
data {
    int N;
    int N_societies;
    int total_tools[N];
    real logpop[N];
    int society[N];
}
parameters{
    real a;
    real bp;
    vector[N_societies] a_society;
    real<lower=0> sigma_society;
}
model{
    vector[N_societies] mu;
    sigma_society ~ cauchy( 0 , 1 );
    a_society ~ normal( 0 , sigma_society );
    bp ~ normal( 0 , 1 );
    a ~ normal( 0 , 10 );
    for ( i in 1:N ) {
        mu[i] = a + a_society[society[i]] + bp * logpop[i];
        mu[i] = exp(mu[i]);
    }
    total_tools ~ poisson( mu );
}
";

# Define the SampleModel.

m12_6s = SampleModel("m12.6s",  stan12_6);

# Input data for cmdstan

m12_6_data = Dict("N" => size(df, 1), "N_societies" => 10,
"total_tools" => df[!, :total_tools], "logpop" => df[!, :log_pop],
"society" => df[!, :society]);

# Sample using cmdstan's sample option

rc12_6s = stan_sample(m12_6s, data=m12_6_data);

# Describe the draws

if success(rc12_6s)
  part12_6s = read_samples(m12_6s, :particles)
  part12_6s |> display
end


text:
using Chain
using DataFrames
using DataFramesMeta
using Dates
using StatsBase

function get_mentions_by(twit_data::DataFrame, date_func::Function, manip_func::Function = mean)
    twit_data = add_mentions(twit_data)

    mentions_by = @chain twit_data begin
        @rtransform :temp = date_func(:date)
        @by(:temp, :nmentions = manip_func(:n_mentions))
    end

    rename!(mentions_by, :temp => Symbol(date_func))

    return mentions_by
end

function add_mentions(twit_data::DataFrame)
    twit_data.mentions = get_mentions.(lowercase.(twit_data.text))
    twit_data.n_mentions = length.(twit_data.mentions)

    return twit_data
end

function make_aog_graph(df::DataFrame)
    aog = data(df) *
        visual(Lines) *
        mapping(Symbol(names(df)[1]), Symbol(names(df)[2]))

    return aog
end

function wrap_in_makie(out_of_AoG, annotations, use_mad = false, fig = Figure(), use_log = false,)
    if use_log
        ax1 = Axis(fig[1,1],
            yscale = log10,
            xscale = log10,
            xlabel = annotations[1],
            ylabel = annotations[2],
            title = annotations[3],
        )
    else
        ax1 = Axis(fig[1,1],
            xlabel = annotations[1],
            ylabel = annotations[2],
            title = annotations[3],
        )
    end

    # dates = string.(out_of_AoG.data[1])
    days = length(out_of_AoG.data[1])
    # out_of_AoG.data.x = collect(1:days)
    ax1.xticks = (out_of_AoG.data[1][1:7:days], out_of_AoG.data[3][1:7:days])

    draw!(ax1, out_of_AoG)

    if use_mad
        mean_data = mean(out_of_AoG.data[2])
        mad_data = mad(out_of_AoG.data[2])

        hlines!(ax1, mean_data)
        hspan!(ax1, mean_data-mad_data , mean_data+mad_data, color = (:blue, 0.2))
    end

    return fig
end

function cumulative_sorting(df,column)
    df = @chain df begin
        @subset(:true_type .!= "retweeted")
        groupby(:author_id)
        combine(column => sum => :tot)
        @orderby(-:tot)
        @transform(_,
        :cumulative = cumsum(:tot),
        :rownumber = rownumber.(eachrow(_))) .+ 1 # .+1 to make space for intercept
        @select :rownumber :cumulative
        rename(:rownumber => :x, :cumulative => :y)
    end
    #we add the intercept at 1,1 (that is 0,0 when logged)
    df = [DataFrame(x=1,y=1);df]
    return df
end

function get_column_summary(twit_data::DataFrame, column::Symbol, date_func::Function, manip_func::Function = sum)
    df = select(twit_data, :date, column)
    rename!(df, column => "data")
    df = @chain df begin
        @rtransform :d = date_func(:date)
        @by(:d, :y = manip_func(:data))
        @rtransform :x = dayofyear(:d)
        @select(:x, :y, :d)
    end

    return df
end

function plot_graph(data::DataFrame, annotations, add_mad::Bool = false, f = Figure(), use_log = false)
    aog = make_aog_graph(data)
    fig = wrap_in_makie(aog, annotations, add_mad, f, use_log)
    return fig
end


text:
@testset "RSwM2" begin

using DiffEqNoiseProcess

W = WienerProcess(0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2))

dt = 0.2
calculate_step!(W,dt)
W.dt,W.curt,W.curW,W.dW

reject_step!(W,0.1)
W.dt,W.curt,W.curW,W.dW
DiffEqNoiseProcess.top(W.S₁)

accept_step!(W,0.2)
W.dt,W.curt,W.curW,W.dW

reject_step!(W,0.1)
W.dt,W.curt,W.curW,W.dW
DiffEqNoiseProcess.top(W.S₁)

for i in 1:100
  reject_step!(W,0.05)
  accept_step!(W,0.1)
end

dt = dt/100
for t in dt:dt:W.t[end]-dt
  W(t)
end

end


text:
module TensorBalancing
export nBalancing, qnBalancing

"""
    calcResidual(P::Array{T,2})

Calculate the residual for balancing.
"""
calcResidual(P, scale=1.0) = norm(vcat((scale * sum(P, 1)' .- 1), (scale * sum(P, 2) .- 1)))

"""
    genβ{T<:AbstractArray}(P::Array{T,2})

Collect indices for β.
"""
function genβ{T<:AbstractFloat}(P::Array{T,2})
    ylen, xlen = size(P)
    return vcat([(1,1)], [(1,i) for i = 2:xlen], [(i,1) for i = 2:ylen])
end

"""
    calcη!{T<:AbstractArray}(η::Array{T,2}, P::Array{T,2})

Update η coordinate of a matrix.
"""
function calcη!{T<:AbstractFloat}(η::Array{T,2}, P::Array{T,2})
    shape = size(P)
    idx = [shape...]
    η[idx...] = P[idx...]
    for i in shape[1]-1:-1:1
        η[i,shape[2]] = η[i+1, shape[2]] + P[i, shape[2]]
    end
    for j in shape[2]-1:-1:1
        η[shape[1], j] = η[shape[1], j+1] + P[shape[1], j]
    end
    for i in shape[1]-1:-1:1
        for j in shape[2]-1:-1:1
            η[i,j] = η[i+1, j] + η[i, j+1] - η[i+1, j+1] + P[i,j] # de morgan
        end
    end
    η ./= η[1,1]
end

"""
    calcMarginals!{T<:AbstractArray}(η::Array{T,1}, A::Array{T,2})

Update η, a list containing the limb part of cumulative sums of a matrix.
"""
function calcMarginals!{T<:AbstractFloat}(η::Array{T,1}, A::Array{T,2})
    ylen, xlen = size(A)
    cumcol = zeros(ylen)
    for i = xlen-1:-1:1
        cumcol += A[:,i+1]
        η[i+1] = sum(cumcol)
    end
    cumcol += A[:,1]
    for i = length(cumcol)-1:-1:1
      cumcol[i] += cumcol[i+1]
    end
    η[xlen+1:end] = cumcol[2:end]
    η[1] = cumcol[1]
end

"""
    genTargetη(P::Array{T}, β)

Generate target value of elements in η.
"""
function genTargetη{T<:AbstractFloat}(P::AbstractArray{T}, β)
    shape = size(P)
    dim = length(β[1])
    return [prod(map(x->x[1]/x[2], zip([shape...]-[idx...]+ones(Int32, dim), shape))) for idx in β]
end

"""
    cumθ{I<:Integer, T<:AbstractFloat}(θ::Array{T,1}, β::Array{Tuple{I,I},1}, shape)

Calculate cumulative sum of θ for each row and column of a matrix.
"""
function cumθ{I<:Integer, T<:AbstractFloat}(θ::Array{T,1}, β::Array{Tuple{I,I},1}, shape)
    ylen, xlen = shape
    cumθ_x = zeros(xlen)
    cumθ_y = zeros(ylen)
    cumθ!(cumθ_x, cumθ_y, θ, β, shape)
    return cumθ_x, cumθ_y
end

"""
    cumθ!{I<:Integer, T<:AbstractFloat}(cumθ_x, cumθ_y, θ::Array{T,1}, β::Array{Tuple{I,I},1}, shape)

Calculate cumulative sum of θ for each row and column of a matrix.
"""
function cumθ!{I<:Integer, T<:AbstractFloat}(cumθ_x::Array{T,1}, cumθ_y::Array{T,1}, θ::Array{T,1}, β::Array{Tuple{I,I},1}, shape)
    ylen, xlen = shape
    for i in 1:length(β)
        if β[i][1] == 1
            cumθ_x[β[i][2]] = θ[i]
        else
            cumθ_y[β[i][1]] = θ[i]
        end
    end
    for i = 2:xlen
        cumθ_x[i] += cumθ_x[i-1]
    end
    for i = 2:ylen
        cumθ_y[i] += cumθ_y[i-1]
    end
end

function genSubβ{T<:AbstractFloat}(β, η::Array{T,2}; ϵ=1e-6)
    subβ_ = sort(β, 1; by=pos->η[pos...], rev=true)
    #subβ = [subβ_[1]]
    subβ = [1]
    for i = 2:length(subβ_)
        if abs(η[subβ[end]...] - η[subβ_[i]...]) > ϵ
            #append!(subβ, [subβ_[i]])
            push!(subβ, i)
        else
            common = tuple(map(x->max(x...), zip(subβ[end], subβ_[i]))...)
            if abs(η[subβ[end]...] - η[common...]) > ϵ
                #append!(subβ, [subβ_[i]])
                push!(subβ, i)
            end
        end
    end
    return subβ
end

"""
    nBalancing{T<:AbstractFloat}(P::Array{T,2}; ϵ=1e-9, max_iter=100, show_trace=false, show_every=100)

Matrix balancing by E-projection with Newton's method.

# Parameters
* P - Input matrix
* ϵ - Threshold for the error of η-coordinate
* max_iter - Max number of iterations for optimization
* show_trace - Whether to print trace
* show_every - Trace information is printed every `show_every`'th iteration
"""
function nBalancing{T<:AbstractFloat}(P::Array{T,2}; ϵ=1e-9, max_iter=100, show_trace=false, show_every=100)
    A = copy(P)
    ylen, xlen = size(A)
    β = genβ(A)
    η_target = genTargetη(A, β)
    lengthβ = length(β)

    # allocate memory
    η = zeros(size(A))
    θ, λ = zeros(lengthβ), 1.0 # θ and Lagrange multiplier
    δθλ = zeros(lengthβ+1)
    grad = zeros(lengthβ+1)
    jacobian = zeros(lengthβ+1, lengthβ+1)
    cumθ_y, cumθ_x = zeros(ylen), zeros(xlen)

    if show_trace
        @printf "Iter     Function value    Gradient norm\n"
    end

    # e-projection loop
    for count = 1:max_iter
        sm = sum(A)
        calcη!(η, A)
        for i = 1:lengthβ
            grad[i] .= (1 + λ * sm) * η[β[i]...] - η_target[i]
        end
        grad[end] = sm - 1

        res = calcResidual(A,xlen)
        if res < ϵ
            gnorm = sqrt(sum(grad[1:end-1].^2))
            show_trace && @printf "%6d   %14e    %14e\n" count-1 log(sm)-η_target'θ gnorm
            break
        end
        if show_trace && (count - 1) % show_every == 0
            gnorm = sqrt(sum(grad[1:end-1].^2))
            @printf "%6d   %14e    %14e\n" count-1 log(sm)-η_target'θ gnorm
        end

        # update jacobian
        for i in 1:lengthβ
            for j in i:lengthβ
                common = tuple(map(x->max(x...), zip(β[i], β[j]))...)
                jacobian[i,j] = (1 + λ * sm) * η[common...] - η[β[i]...] * η[β[j]...]
                jacobian[j,i] = jacobian[i,j]
            end
            jacobian[i,end] = jacobian[end,i] = sm * η[β[i]...]
        end

        # update parameters
        try
            δθλ .= jacobian \ grad # solve a linear equation
        catch e
            warn("Jacobian matrix got singular. Projection loop is terminated.")
            break
        end
        θ .-= δθλ[1:end-1]
        λ -= δθλ[end]
        # make cumulative sum of θ for fast update
        cumθ!(cumθ_x, cumθ_y, θ, β, size(A))
        # update matrix
        for j = 1:xlen, i = 1:ylen
            if P[i,j] > 1e-15
                A[i,j] = P[i,j] * exp(cumθ_y[i] + cumθ_x[j])
            end
        end
    end

    return A
end

include("QN.jl")

end # module


text:
using ModelerToolbox, JLD, HDF5
using Base.Test

type Lemons
    n::Int
    color::Symbol
end

ModelerToolbox.jld_fn(::Lemons) = Pkg.dir("ModelerToolbox", "test", "lemons.jld")
ModelerToolbox.group_keys(l::Lemons) = (l.n, l.color)

type Apples
    n::Int
    color::Symbol
end

ModelerToolbox.jld_fn(::Apples) = Pkg.dir("ModelerToolbox", "test", "apples.jld")
ModelerToolbox.group_keys(a::Apples) = a.n

# create some test objects
l1 = Lemons(1, :yellow)
l2 = Lemons(3, :yellow)
l3 = Lemons(1, :greenish)
a1 = Apples(1, :red)
a2 = Apples(3, :yellow)
a3 = Apples(1, :green)  # NOTE: this should have the same key as a1

# delete the files if they already exist
isfile(jld_fn(l1)) && rm(jld_fn(l1))
isfile(jld_fn(a1)) && rm(jld_fn(a1))

# make sure we can init the file
init_jld(l1)
@test isfile(jld_fn(l1))
@test !isfile(jld_fn(a1))

init_jld(a1)
@test isfile(jld_fn(a1))

# make sure the group was created for the first object
@test groupname(l1) == ("group1", true, true)
@test groupname(l2) == ("group2", false, false)
# NOTE: repeat group2 because it was never created
@test groupname(l3) == ("group2", false, false)

# make sure the group was created for the first object
@test groupname(a1) == ("group1", true, true)
@test groupname(a2) == ("group2", false, false)
@test groupname(a3) == ("group1", true, true)

# now make all the groups
@test groupname!(l1) == ("group1", true, true)
@test groupname!(l2) == ("group2", true, true)
@test groupname!(l3) == ("group3", true, true)

# make sure the group was created for the first object
@test groupname!(a1) == ("group1", true, true)
@test groupname!(a2) == ("group2", true, true)
@test groupname!(a3) == ("group1", true, true)

# test with_group and with_file
@test length(with_group(names, l1)) == 1
@test sort(with_file(names, l1)) == ["group1", "group2", "group3", "groups"]

# test delete_group!
delete_group!(l1)
jldopen(jld_fn(l1), "r") do f
    @test !exists(f, "group1")

    # also make sure it was removed from the map
    group_map = read(f, "groups")
    @test !haskey(group_map, group_keys(l1))
end

# test write_jld
write_jld(l1, "foo", "bar")
jldopen(jld_fn(l1), "r") do f
    # groupname for l1 is now 4 because we deleted group1 above
    @test exists(f, "group4")
    g = f["group4"]
    @test exists(g, "foo")
    @test read(g, "foo") == "bar"
end

# test read_jld
@test read_jld(l1, "foo") == "bar"

# test delete_jld!
delete_jld!(l1, "foo")
jldopen(jld_fn(l1), "r") do f
    # groupname for l1 is now 4 because we deleted group1 above
    g = f["group4"]
    @test !exists(g, "foo")
end

@test_throws ErrorException delete_jld!(l1, "foo")

# test haskey_jld
@test !haskey_jld(l1, "foo")
write_jld(l1, "bing", "bong")
@test haskey_jld(l1, "bing")

# test names_jld
@test sort(names_jld(l1)) == ["bing", "visted"]

# test groups
grps = groups(l1)
@test length(grps) == 3
@test sort(collect(values(grps))) == ["group2", "group3", "group4"]
@test sort(collect(values(groups(a1)))) == ["group1", "group2"]

# test write_jld!
write_jld!(l1, bing=42)
@test read_jld(l1, "bing") == 42

write_jld!(l1, hello="world")
@test read_jld(l1, "hello") == "world"


text:
repo_name = "Test41"
copy_package("JuliennedArrays", repo_name)

username = settings("username")
github_token = settings("github_token")
ssh_keygen_file = settings("ssh_keygen_file")
set_up(username, github_token, ssh_keygen_file = ssh_keygen_file)

github, travis = put_online(repo_name)
delete(github)


text:
# see: https://stackoverflow.com/a/3272490/1346276

abstract type Tree{N} end

struct Leave{N} <: Tree{N}
    label::Int
end

struct Branch{N} <: Tree{N}
    left::Tree{N}
    right::Tree{N}
end

abstract type Context{N} end

struct LeftContext{N} <: Context{N}
    right::Tree{N}
    parent::Context{N}
end

function (context::LeftContext{N})(t::Tree{N}) where N
    context.parent(Branch(t, context.right))
end

struct RightContext{N} <: Context{N}
    left::Tree{N}
    parent::Context{N}
end

function (context::RightContext{N})(t::Tree{N}) where N
    context.parent(Branch(context.left, t))
end

struct NoContext{N} <: Context{N} end

function (context::NoContext{N})(t::Tree{N}) where N
    t
end

struct Cont{N}
    context::Context{N}
    chunk::Tree{N}
end

function (continuation::Cont{N})(k::Function) where N
    newchunk = k(continuation.chunk)::Tree{N}
    tree = continuation.context(newchunk)::Tree{N}
    return tree, continuation.chunk
end

function randomsplit_impl{N}(t::Leave{N}, n::Int, context::Context{N}, continuation::Cont{N})
    n += 1
    if rand() ≤ 1/n
        return Cont{N}(context, t), n
    else
        return continuation, n
    end
end

function randomsplit_impl{N}(t::Branch{N}, n::Int, context::Context{N}, continuation::Cont{N})
    c1, n = randomsplit_impl(t.left, n, LeftContext{N}(t.right, context), continuation)
    c2, n = randomsplit_impl(t.right, n, RightContext{N}(t.left, context), c1)

    n += 1
    if rand() ≤ 1/n
        return Cont{N}(context, t), n
    else
        return c2, n
    end
end

function randomsplit{N}(action, t::Tree{N})
    default_context = NoContext{N}()
    cont, _ = randomsplit_impl(t, 0, default_context, Cont{N}(default_context, t))
    cont(action)
end
randomchild{N}(t::Tree{N}) = randomsplit(identity, t)[2]

const t = Branch{1}(Branch{1}(Leave{1}(1), Leave{1}(2)), Leave{1}(4))
const t2 = Branch{1}(Leave{1}(1), Leave{1}(2))

# to test for a uniform distribution of splits, we use some hackish hashing by summation:
sumtree(t::Leave) = t.label
sumtree(t::Branch) = sumtree(t.left) + sumtree(t.right)

# function testvalues(N)
#     zd = [sumtree(randomsplit(t, x -> (Leave(64), nothing))[1]) for _ in 1:N]
#     [(i, sum(zd .== i)) for i in unique(zd)]
# end

# k -> Branch(k(Leave()), Leave())
# k -> Branch(Leave(), k(Leave()))
# k -> k(Branch(Leave(), Leave()))

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
using CubicEoSDatabase
using Test

#####
##### Tests
#####

include("types.jl")
include("errors.jl")
include("requests.jl")


text:
using Test, SpecialFunctions, Statistics

function test_digitize(sep; C=1000, K=1000)
    cal0 = randn(C)
    cal1 = randn(C) .+ sep

    # flip 1000 coins, and add noise
    rbits  = round.(Int, rand(K))
    rvolts = sep*rbits + randn(K)

    # estimate of prob of being correct
    dm, pe0m, pe1m = digitize(rvolts, cal0, cal1, mode=:max)
    de, pe0e, pe1e = digitize(rvolts, cal0, cal1, mode=:equal)
    # compute theoretical prob of being correct
    pe = 1/2 - 1/2*erf(sep/(2*sqrt(2)))

    return pe,
           sum(abs.(dm-rbits))/length(rbits),
           sum(abs.(de-rbits))/length(rbits),
           pe0m,
           pe1m,
           pe0e,
           pe1e

end

# run a bunch of tests, store each as a column of a matrix, take the mean
# the rows are:
#   - predicted error rate (asymptotic)
#   - observed error rate for :max mode
#   - observed error rate for :equal mode
#   - estimated error rate based on cals for :max mode (0 preps)
#   - estimated error rate based on cals for :max mode (1 preps)
#   - estimated error rate based on cals for :equal mode (0 preps)
#   - estimated error rate based on cals for :equal mode (1 preps)
test1 = mean(reduce(hcat,[[test_digitize(0.75)...] for _ in 1:1000]))
test2 = mean(reduce(hcat,[[test_digitize(0.50)...] for _ in 1:1000]))
# make sure these are all close to the prediction
@test all(map(x->isapprox(x,0.35,atol=5e-2),test1))
@test all(map(x->isapprox(x,0.40,atol=5e-2),test2))


text:
##############################################################################
##
## TODO
##
##############################################################################

# - Test different WIF types
# - Add negative WIF tests (e.g., invalid compression bits)
# - Test different key types (e.g., compressed)
# - Add negative key tests
# - add test for leading 0
# - to_varstring

using Coin
using Base.Test

##############################################################################
##
## Base58 tests
##
##############################################################################

base58data = parseint(BigInt, "800c28fca386c7a227600b2fe50b7cae11ec86d3bf1fbe471be89827e19d72aa1d507a5b8d", 16)

# Base 58 encoding
@test encode58(base58data) == "5HueCGU8rMjxEXxiPuD5BDku4MkFqeZyd4dZ1jvhTVqvbTLvyTJ"

# Base 58 decoding
@test decode58("5HueCGU8rMjxEXxiPuD5BDku4MkFqeZyd4dZ1jvhTVqvbTLvyTJ") == base58data

##############################################################################
##
## Key generation tests
##
##############################################################################

secret_key = "18E14A7B6A307F426A94F8114701E7C8E774E7F9A47E2C2035DB29A206321725"
public_key = get_pub_key(secret_key)
@test public_key == "16UwLL9Risc3QfPqBUvKofHmBQ7wMtjvM"

##############################################################################
##
## Wallet Interchange Format tests
##
##############################################################################

# Private key to WIF
@test private2wif("0c28fca386c7a227600b2fe50b7cae11ec86d3bf1fbe471be89827e19d72aa1d") == "5HueCGU8rMjxEXxiPuD5BDku4MkFqeZyd4dZ1jvhTVqvbTLvyTJ"

# WIF to private key
@test wif2private("5HueCGU8rMjxEXxiPuD5BDku4MkFqeZyd4dZ1jvhTVqvbTLvyTJ") == "0c28fca386c7a227600b2fe50b7cae11ec86d3bf1fbe471be89827e19d72aa1d"

# WIF checksum
@test wif_check_sum("5HueCGU8rMjxEXxiPuD5BDku4MkFqeZyd4dZ1jvhTVqvbTLvyTJ")

@test pub2base58("0450863AD64A87AE8A2FE83C1AF1A8403CB53F53E486D8511DAD8A04887E5B23522CD470243453A299FA9E77237716103ABC11A1DF38855ED6F2EE187E9C582BA6") == "16UwLL9Risc3QfPqBUvKofHmBQ7wMtjvM"

##############################################################################
##
## Transaction tests
##
##############################################################################

# Test header generation
payload = "0100000001484d40d45b9ea0d652fca8258ab7caa42541eb52975857f96fb50cd732c8b481000000008a47304402202cb265bf10707bf49346c3515dd3d16fc454618c58ec0a0ff448a676c54ff71302206c6624d762a1fcef4618284ead8f08678ac05b13c84235f1654e6ad168233e8201410414e301b2328f17442c0b8310d787bf3d8a404cfbd0704f135b6ad4b2d3ee751310f981926e53a6e8c39bd7d3fefd576c543cce493cbac06388f2651d1aacbfcdffffffff0162640100000000001976a914c8e90996c7c6080ee06284600c684ed904d14c5c88ac00000000"
# @test create_header(magic_mainnet, "tx", payload) == "f9beb4d9747800000000000000000000df000000ea0f5494"

# Transaction output
tx_out = Tx_Output(123, "abc")
@test bytearray(tx_out) == [0x7b, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x0a, 0xbc]
tx_out = Tx_Output(5000000, "76A9141AA0CD1CBEA6E7458A7ABAD512A9D9EA1AFB225E88AC")
@test bytearray(tx_out) == [0x40, 0x4b, 0x4c, 0x00, 0x00, 0x00, 0x00, 0x00, 0x19, 0x76, 0xa9, 0x14, 0x1a, 0xa0, 0xcd, 0x1c, 0xbe, 0xa6, 0xe7, 0x45, 0x8a, 0x7a, 0xba, 0xd5, 0x12, 0xa9, 0xd9, 0xea, 0x1a, 0xfb, 0x22, 0x5e, 0x88, 0xac]
tx_out = Tx_Output(91234, [0x76, 0xa9, 0x14, 0xc8, 0xe9, 0x09, 0x96, 0xc7, 0xc6, 0x08, 0x0e, 0xe0, 0x62, 0x84, 0x60, 0x0c, 0x68, 0x4e, 0xd9, 0x04, 0xd1, 0x4c, 0x5c, 0x88, 0xac])

# Transaction outpoint (i.e., reference to previous transaction)
@test_throws ErrorException OutPoint("abc", 1)
outpoint = OutPoint("abcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabcaabca", 1)
@test bytearray(outpoint) == [0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0xca, 0xab, 0x01, 0x00, 0x00, 0x00]
outpoint = OutPoint("81b4c832d70cb56ff957589752eb4125a4cab78a25a8fc52d6a09e5bd4404d48", 0)
@test bytearray(outpoint) == [0x48, 0x4d, 0x40, 0xd4, 0x5b, 0x9e, 0xa0, 0xd6, 0x52, 0xfc, 0xa8, 0x25, 0x8a, 0xb7, 0xca, 0xa4, 0x25, 0x41, 0xeb, 0x52, 0x97, 0x58, 0x57, 0xf9, 0x6f, 0xb5, 0x0c, 0xd7, 0x32, 0xc8, 0xb4, 0x81, 0x00, 0x00, 0x00, 0x00]

# Transaction input
scriptSig = "47304402202cb265bf10707bf49346c3515dd3d16fc454618c58ec0a0ff448a676c54ff71302206c6624d762a1fcef4618284ead8f08678ac05b13c84235f1654e6ad168233e8201410414e301b2328f17442c0b8310d787bf3d8a404cfbd0704f135b6ad4b2d3ee751310f981926e53a6e8c39bd7d3fefd576c543cce493cbac06388f2651d1aacbfcd"
tx_in = Tx_Input(outpoint, scriptSig)
@test bytearray(tx_in) == [0x48, 0x4d, 0x40, 0xd4, 0x5b, 0x9e, 0xa0, 0xd6, 0x52, 0xfc, 0xa8, 0x25, 0x8a, 0xb7, 0xca, 0xa4, 0x25, 0x41, 0xeb, 0x52, 0x97, 0x58, 0x57, 0xf9, 0x6f, 0xb5, 0x0c, 0xd7, 0x32, 0xc8, 0xb4, 0x81, 0x00, 0x00, 0x00, 0x00, 0x8a, 0x47, 0x30, 0x44, 0x02, 0x20, 0x2c, 0xb2, 0x65, 0xbf, 0x10, 0x70, 0x7b, 0xf4, 0x93, 0x46, 0xc3, 0x51, 0x5d, 0xd3, 0xd1, 0x6f, 0xc4, 0x54, 0x61, 0x8c, 0x58, 0xec, 0x0a, 0x0f, 0xf4, 0x48, 0xa6, 0x76, 0xc5, 0x4f, 0xf7, 0x13, 0x02, 0x20, 0x6c, 0x66, 0x24, 0xd7, 0x62, 0xa1, 0xfc, 0xef, 0x46, 0x18, 0x28, 0x4e, 0xad, 0x8f, 0x08, 0x67, 0x8a, 0xc0, 0x5b, 0x13, 0xc8, 0x42, 0x35, 0xf1, 0x65, 0x4e, 0x6a, 0xd1, 0x68, 0x23, 0x3e, 0x82, 0x01, 0x41, 0x04, 0x14, 0xe3, 0x01, 0xb2, 0x32, 0x8f, 0x17, 0x44, 0x2c, 0x0b, 0x83, 0x10, 0xd7, 0x87, 0xbf, 0x3d, 0x8a, 0x40, 0x4c, 0xfb, 0xd0, 0x70, 0x4f, 0x13, 0x5b, 0x6a, 0xd4, 0xb2, 0xd3, 0xee, 0x75, 0x13, 0x10, 0xf9, 0x81, 0x92, 0x6e, 0x53, 0xa6, 0xe8, 0xc3, 0x9b, 0xd7, 0xd3, 0xfe, 0xfd, 0x57, 0x6c, 0x54, 0x3c, 0xce, 0x49, 0x3c, 0xba, 0xc0, 0x63, 0x88, 0xf2, 0x65, 0x1d, 0x1a, 0xac, 0xbf, 0xcd, 0xff, 0xff, 0xff, 0xff]

# Transaction
tx = Tx([tx_in], [tx_out])
@test bytearray(tx) == [0x01, 0x00, 0x00, 0x00, 0x01, 0x48, 0x4d, 0x40, 0xd4, 0x5b, 0x9e, 0xa0, 0xd6, 0x52, 0xfc, 0xa8, 0x25, 0x8a, 0xb7, 0xca, 0xa4, 0x25, 0x41, 0xeb, 0x52, 0x97, 0x58, 0x57, 0xf9, 0x6f, 0xb5, 0x0c, 0xd7, 0x32, 0xc8, 0xb4, 0x81, 0x00, 0x00, 0x00, 0x00, 0x8a, 0x47, 0x30, 0x44, 0x02, 0x20, 0x2c, 0xb2, 0x65, 0xbf, 0x10, 0x70, 0x7b, 0xf4, 0x93, 0x46, 0xc3, 0x51, 0x5d, 0xd3, 0xd1, 0x6f, 0xc4, 0x54, 0x61, 0x8c, 0x58, 0xec, 0x0a, 0x0f, 0xf4, 0x48, 0xa6, 0x76, 0xc5, 0x4f, 0xf7, 0x13, 0x02, 0x20, 0x6c, 0x66, 0x24, 0xd7, 0x62, 0xa1, 0xfc, 0xef, 0x46, 0x18, 0x28, 0x4e, 0xad, 0x8f, 0x08, 0x67, 0x8a, 0xc0, 0x5b, 0x13, 0xc8, 0x42, 0x35, 0xf1, 0x65, 0x4e, 0x6a, 0xd1, 0x68, 0x23, 0x3e, 0x82, 0x01, 0x41, 0x04, 0x14, 0xe3, 0x01, 0xb2, 0x32, 0x8f, 0x17, 0x44, 0x2c, 0x0b, 0x83, 0x10, 0xd7, 0x87, 0xbf, 0x3d, 0x8a, 0x40, 0x4c, 0xfb, 0xd0, 0x70, 0x4f, 0x13, 0x5b, 0x6a, 0xd4, 0xb2, 0xd3, 0xee, 0x75, 0x13, 0x10, 0xf9, 0x81, 0x92, 0x6e, 0x53, 0xa6, 0xe8, 0xc3, 0x9b, 0xd7, 0xd3, 0xfe, 0xfd, 0x57, 0x6c, 0x54, 0x3c, 0xce, 0x49, 0x3c, 0xba, 0xc0, 0x63, 0x88, 0xf2, 0x65, 0x1d, 0x1a, 0xac, 0xbf, 0xcd, 0xff, 0xff, 0xff, 0xff, 0x01, 0x62, 0x64, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x19, 0x76, 0xa9, 0x14, 0xc8, 0xe9, 0x09, 0x96, 0xc7, 0xc6, 0x08, 0x0e, 0xe0, 0x62, 0x84, 0x60, 0x0c, 0x68, 0x4e, 0xd9, 0x04, 0xd1, 0x4c, 0x5c, 0x88, 0xac, 0x00, 0x00, 0x00, 0x00]

priv_key = Crypto.random(256)
pub_key = Crypto.ec_pub_key(priv_key)
keys = [Keys(priv_key, pub_key)]

outpoints = [OutPoint("81b4c832d70cb56ff957589752eb4125a4cab78a25a8fc52d6a09e5bd4404d48", 1)]
addresses = ["1KKKK6N21XKo48zWKuQKXdvSsCf95ibHFa"]
amounts = [91234]
tx = create_tx(keys, outpoints, addresses, amounts)

##############################################################################
##
## Utility tests
##
##############################################################################

# Reverse endian of hex string
@test reverse_endian("") == ""
@test reverse_endian("0c28fca386c7a227600b2fe50b7cae11ec86d3bf1fbe471be89827e19d72aa1d") == "1daa729de12798e81b47be1fbfd386ec11ae7c0be52f0b6027a2c786a3fc280c"
@test reverse_endian(0) == 0
@test reverse_endian(0x22) == 0x22
@test reverse_endian(0x333) == 0x3303
@test reverse_endian(0x2234) == 0x3422
@test reverse_endian(0xf9beb4d9) == 0xd9b4bef9
@test reverse_endian(BigInt(19238471923847192837419283749128374912837491823742198374)) == 9830438508025927557749821611753759704769006686345681864
# reverse_endian assumes that BigInt is at least 16 bytes and Int is system-dependent
@test reverse_endian(BigInt(24)) == 31901471898837980949691369446728269824
@test reverse_endian(23) == 1657324662872342528
@test reverse_endian(-23) == -1585267068834414593

# Test VarInt conversion
@test to_varint(1) == [0x01]
@test to_varint(252) == [0xfc]
@test to_varint(253) == [0xfd, 0x00, 0xfd]
@test to_varint(0xffff) == [0xfd, 0xff, 0xff]
@test to_varint(0x10000) == [0xfe, 0x00, 0x01, 0x00, 0x00]
@test to_varint(0x100000000) == [0xff, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00]


text:
function generateall(input)
    cubes = Set{Tuple{Int,Int,Int}}()
    hypercubes = Set{Tuple{Int,Int,Int,Int}}()
    for (i, r) in enumerate(input)
        for (j, c) in enumerate(r)
            if c == '#'
                push!(cubes, (i, j, 0))
                push!(hypercubes, (i, j, 0, 0))
            end
        end
    end
    return cubes, hypercubes
end

function generatecuberange()
    r = (-1:1)
    return collect((x, y, z) for x in r for y in r for z in r if x != 0 || y != 0 || z != 0)
end

function generatehypercuberange()
    r = (-1:1)
    return collect((x, y, z, w) for x in r for y in r for z in r for w in r if x != 0 || y != 0 || z != 0 || w != 0)
end

function nextcubes(cubes)
    newcubes = deepcopy(cubes)
    marked = Dict()
    for (x, y, z) in cubes
        neighbors = 0
        for (dx, dy, dz) in generatecuberange()
            cube = (x + dx, y + dy, z + dz)
            if ∈(cube, cubes)
                neighbors += 1
            else
                push!(marked, cube => (get!(marked, cube, 0) + 1))
            end
        end
        if ∉(neighbors, [2, 3])
            pop!(newcubes, (x, y, z))
        end
    end
    for (cube, count) in marked
        if count == 3
            push!(newcubes, cube)
        end
    end
    return newcubes
end

function nexthypercubes(hypercubes)
    newhypercubes = deepcopy(hypercubes)
    marked = Dict()
    for (x, y, z, w) in hypercubes
        neighbors = 0
        for (dx, dy, dz, dw) in generatehypercuberange()
            cube = (x + dx, y + dy, z + dz, w + dw)
            if ∈(cube, hypercubes)
                neighbors += 1
            else
                push!(marked, cube => (get!(marked, cube, 0) + 1))
            end
        end
        if ∉(neighbors, [2, 3])
            pop!(newhypercubes, (x, y, z, w))
        end
    end
    for (cube, count) in marked
        if count == 3
            push!(newhypercubes, cube)
        end
    end
    return newhypercubes
end

function cycle(something, next, cycle = 6)
    i = 0
    while i < cycle
        something = next(something)
        i += 1
    end
    return something
end

input = readlines("d17.in")
cubes, hypercubes = generateall(input)
cubes = cycle(cubes, nextcubes)
hypercubes = cycle(hypercubes, nexthypercubes)

println("P1: ", length(cubes))
println("P2: ", length(hypercubes))

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
"""
    struct BlockCollection{TGP,TGPl,TW,TWn,TX}
        blocks::Vector{BiBlock{_L,TGP,TGPl,TW,TWn,TX} where _L}
    end

Gathers all blocks relevant for an entire single recording.

    BlockCollection(
        sp::SamplingPair{TGP,TGPl,TW,TWn,TX},
        ranges,
        ρρ=0.0,
        ll_hist_len=0
    ) where {TGP,TGPl,TW,TWn,TX}

Base constructor.
"""
struct BlockCollection{TGP,TGPl,TW,TWn,TX}
    blocks::Vector{BiBlock{_L,TGP,TGPl,TW,TWn,TX} where _L}

    function BlockCollection(
            sp::SamplingPair{TGP,TGPl,TW,TWn,TX},
            ranges,
            ρρ=0.0,
            ll_hist_len=0
        ) where {TGP,TGPl,TW,TWn,TX}
        N = length(ranges)
        ρρ = (typeof(ρρ) <: Number ? fill(ρρ, N) : ρρ)
        blocks = map(i->BiBlock(sp, ranges[i], ρρ[i], i==N, ll_hist_len), 1:N)
        new{TGP,TGPl,TW,TWn,TX}(blocks)
    end
end

#===============================================================================

        IMPUTATION

===============================================================================#
"""
    draw_proposal_path!(bc::BlockCollection)

Sample proposal paths on each block, compute log-likelihoods along the way.
Assumes `bc.blocks[i].b.XX[1].x[1]` are starting points for each block. Uses
the preconditioned Crank-Nicolson scheme.
"""
draw_proposal_path!(bc::BlockCollection) = draw_proposal_path!.(bc.blocks)

#===============================================================================

        ACCEPT/REJECT DECISION

===============================================================================#

"""
    accept_reject_proposal_path!(bc::BlockCollection, mcmciter)

Accept/reject decision of the Metropolis-Hastings algorithm for the step of path
imputation, done separately for each block.
"""
function accept_reject_proposal_path!(bc::BlockCollection, mcmciter)
    for biblock in bc.blocks
        accept_reject_proposal_path!(biblock, mcmciter)
    end
end

#===============================================================================

        SWAPS

===============================================================================#

"""
    swap_paths!(bc::BlockCollection)

For each block in the collection swap `XX` and `WW` containers between
proposal-acceptance pair.
"""
swap_paths!(bc::BlockCollection) = swap_paths!.(bc.blocks)

"""
    swap_XX!(bc::BlockCollection)

For each block in the collection swap `XX` containers between
proposal-acceptance pair.
"""
swap_XX!(bc::BlockCollection) = swap_XX!.(bc.blocks)

"""
    swap_WW!(bc::BlockCollection)

For each block in the collection swap `WW` containers between
proposal-acceptance pair.
"""
swap_WW!(bc::BlockCollection) = swap_WW!.(bc.blocks)

"""
    swap_PP!(bc::BlockCollection)

For each block in the collection swap `PP` containers between
proposal-acceptance pair.
"""
swap_PP!(bc::BlockCollection) = swap_PP!.(bc.blocks)

"""
    swap_ll!(bc::BlockCollection)

For each block in the collection swap `ll` fields between
proposal-acceptance pair.
"""
swap_ll!(bc::BlockCollection) = swap_ll!.(bc.blocks)

#===============================================================================

        UTILITY

===============================================================================#

"""
    loglikhd!(bc::BlockCollection)

For each `BiBlock` in a collection compute the log-likelihood for the accepted
block, evaluated at a sampled path and store the result in internal fields `ll`.
"""
loglikhd!(bc::BlockCollection) = loglikhd!.(bc.blocks)

"""
    loglikhd°!(bc::BlockCollection)

For each `BiBlock` in a collection compute the log-likelihood for the proposal
block, evaluated at a sampled path and store the result in internal fields `ll`.
"""
loglikhd°!(bc::BlockCollection) = loglikhd°!.(bc.blocks)

@doc raw"""
    fetch_ll(bc::BlockCollection)

Retreive the log-likelihood for the entire accepted path.
!!! warning
    The function uses only internal fields `ll` for this computation, which
    means that for the call to this function to make sense the log-likelihood
    must have been previously computed and stored in the field `ll`. If it
    hasn't been done, then you must first call `loglikhd!(bc)`.
"""
fetch_ll(bc::BlockCollection) = mapreduce(bb->bb.b.ll, +, bc.blocks)

@doc raw"""
    fetch_ll°(bc::BlockCollection)

Retreive the log-likelihood for the entire proposed path.
!!! warning
    The function uses only internal fields `ll` for this computation, which
    means that for the call to this function to make sense the log-likelihood
    must have been previously computed and stored in the field `ll`. If it
    hasn't been done, then you must first call `loglikhd°!(bc)`.
"""
fetch_ll°(bc::BlockCollection) = mapreduce(bb->bb.b°.ll, +, bc.blocks)

"""
    save_ll!(bc::BlockCollection, i::Int)

For each block in the collection commit the current proposal and accepted
log-likelihood fields `ll` to history, at index `i`.
"""
save_ll!(bc::BlockCollection, i::Int) = save_ll!.(bc.blocks, i)

"""
    ll_of_accepted(bb::BlockCollection, i)

Return an array of log-likelihoods (one for each block) of the paths that were
accepted at the `i`th iteration.
"""
ll_of_accepted(bc::BlockCollection, i) = ll_of_accepted.(bc.blocks, i)

"""
    accpt_rate(bb::BlockCollection, range)

Compute the acceptance rate over the `range` of MCMC accept/reject history for
each block in the collection.
"""
function accpt_rate(bc::BlockCollection, range)
    map(bc.blocks) do biblock
        accpt_rate(biblock, range)
    end
end

#===============================================================================

        SETTING ARTIFICIAL OBSERVATIONS

===============================================================================#

"""
    GP.set_obs!(bc::BlockCollection)

Freeze an artificial observation at the terminal point of each block. For a
terminal block nothing is done.
"""
GP.set_obs!(bc::BlockCollection) = set_obs!.(bc.blocks)

"""
    GP.recompute_guiding_term!(bc::BlockCollection, [::Val{:_only}])

For each block in the collection recompute the guiding terms of both the
proposal and the accepted laws. If an additional flag `Val(:P_only)` is passed,
then recomputes the guiding term on accepted law only. If `Val(:P°_only)`, then
recomputes for the proposal law only.
"""
function GP.recompute_guiding_term!(bc::BlockCollection)
    recompute_guiding_term!.(bc.blocks)
end

function GP.recompute_guiding_term!(bc::BlockCollection, ::Val{:P_only})
    _rgt_only_P!.(bc.blocks)
end

function GP.recompute_guiding_term!(bc::BlockCollection, ::Val{:P°_only})
    _rgt_only_P°!.(bc.blocks)
end

_rgt_only_P!(bb::BiBlock) = recompute_guiding_term!(bb.b)
_rgt_only_P°!(bb::BiBlock) = recompute_guiding_term!(bb.b°)

"""
    find_W_for_X!(bc::BlockCollection)

For each block find the Wiener process `bb.b.WW` that reconstructs the path
`bb.b.XX` under the accepted law `bb.b.PP` (possibly including `bb.b.P_last`).
"""
find_W_for_X!(bc::BlockCollection) = find_W_for_X!.(bc.blocks)

#===============================================================================

        SETTING PARAMETERS

===============================================================================#

"""
    is_critical_update(bb::BlockCollection, pnames)

For each block in the collection verify whether the update characterized by a
list of parameter names stored in `pnames` is `critical` in a sense of prompting
for recomputation of the guiding term.
"""
GP.is_critical_update(bc::BlockCollection, pnames) = map(
    bb->GP.is_critical_update(bb.b.PP, pnames.θ_local_aux, pnames.θ_local_obs),
    bc.blocks
)

"""
    set_proposal_law!(
        bc::BlockCollection,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bc, pnames);
        skip=0
    )

For each block in the collection set the parameters in `bb.b°.PP` and
`bb.b°.P_last` to `θ°` and make sure that all other parameters are shared with
`bb.b.PP` and `bb.b.P_last`. Recompute the guiding term if needed, and then,
compute the proposal trajectory `bb.b°.XX` for the proposal point `θ°`.
"""
function set_proposal_law!(
        bc::BlockCollection,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bc, pnames);
        skip=0
    )
    for i in eachindex(bc.blocks)
        set_proposal_law!(
            bc.blocks[i], θ°, pnames.blocks[i], critical_change[i]; skip=skip
        )
    end
end


text:
module SolverInterface

using Base.Meta
using Compat
using Compat.SparseArrays
using Compat.LinearAlgebra

const methods_by_tag = Dict{Symbol,Vector{Symbol}}()

macro define_interface(args)
    @assert args.head == :block
    code = quote end
    for line in args.args
        (isexpr(line, :line) || isa(line, LineNumberNode)) && continue
        (isa(line,Symbol) || isexpr(line,:tuple)) || error("Unexpected code in block")
        if isa(line,Symbol)
            fname = line
            tags = []
        else
            @assert isexpr(line,:tuple)
            fname = line.args[1]
            tags = line.args[2:end]
        end
        func = esc(fname)
        code = quote
            $code
            $(func)() = throw(MethodError($(func),()))
            function $(func)(::Int)
                q = Any[1]
                return q[1] # infer Any return type, workaround for julia issue #9479
            end
        end
        for t in tags
            t = quot(t)
            code = quote
                $code
                if haskey(methods_by_tag,$t)
                    push!(methods_by_tag[$t],$(quot(fname)))
                else
                    methods_by_tag[$t] = [$(quot(fname))]
                end
            end
        end
        push!(code.args, Expr(:export, func))
    end
    return code
end

abstract type AbstractMathProgModel end
export AbstractMathProgModel

# immutable type which we dispatch solvers on
abstract type AbstractMathProgSolver end
export AbstractMathProgSolver

# create dummy method to define function so that we can attach methods in other modules
# these are the common methods for AbstractMathProgModel
@define_interface begin
    getsolution
    getobjval
    optimize!
    status
    getobjbound, rewrap
    getobjgap, rewrap
    getrawsolver, rewrap
    getsolvetime, rewrap
    setsense!
    getsense
    numvar
    numconstr
    freemodel!
    setvartype!
    getvartype
    loadproblem!
    setwarmstart!
end

# solver parameters, may be implemented by AbstractMathProgModel or AbstractMathProgSolver
@define_interface begin
    setparameters!
end

include("LinearQuadratic.jl")
include("callbacks.jl")
include("Nonlinear.jl")
include("Conic.jl")

# Solver conversion routines
include("conic_to_lpqp.jl")
include("lpqp_to_conic.jl")
include("nonlinear_to_lpqp.jl")

end


text:
module UnitfulOptim

using Unitful
import Unitful: Quantity, Units
import Optim: OptimizationResults
import Optim: optimize, minimizer, minimum

struct UnitfulOptimizationResults
    opt::OptimizationResults
    xunit::Units
    funit::Units
end

function optimize(f, a::Quantity, b::Quantity; kws...)
    c, d = promote(a, b)
    optimize(f, c, d; kws...)
end

function optimize{T<:Real,D,U}(f, a::Quantity{T,D,U}, b::Quantity{T,D,U}; kws...)
    g(x::T) = ustrip(f(x*U()))
    opt = optimize(g, ustrip(a), ustrip(b); kws...)
    UnitfulOptimizationResults(opt, U(), unit(f(a)))
end

for f in (:minimizer,)
    @eval function $f(uopt::UnitfulOptimizationResults)
        $f(uopt.opt) * uopt.xunit
    end
end

for f in (:minimum,)
    @eval function $f(uopt::UnitfulOptimizationResults)
        $f(uopt.opt) * uopt.funit
    end
end

end # module


text:
### Latent models ####
include("prior.jl")
include("posterior.jl")

## Exact Gaussian Process
struct LatentGP{T,Tpr<:GPPrior,Tpo<:Posterior{T},O} <: AbstractLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function LatentGP(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return LatentGP(
        GPPrior(deepcopy(kernel), deepcopy(mean), cholesky(Matrix{T}(I, dim, dim))),
        Posterior(dim, zeros(T, dim), cholesky(Matrix{T}(I(dim)))),
        deepcopy(opt),
    )
end

@traitimpl IsFull{LatentGP}

## AbstractVarLatent

abstract type AbstractVarLatent{T,Tpr,Tpo} <: AbstractLatent{T,Tpr,Tpo} end

## Variational Gaussian Process
mutable struct VarLatent{T,Tpr<:GPPrior,Tpo<:VarPosterior{T},O} <:
               AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function VarLatent(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return VarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean), cholesky(Matrix{T}(I, dim, dim))),
        VarPosterior{T}(dim),
        deepcopy(opt),
    )
end

@traitimpl IsFull{VarLatent}

## Sparse Variational Gaussian Process

mutable struct SparseVarLatent{
    T,Tpr<:GPPrior,Tpo<:VarPosterior{T},Topt,TZ<:AbstractVector,TZopt
} <: AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    Z::TZ
    Knm::Matrix{T}
    κ::Matrix{T}
    K̃::Vector{T}
    opt::Topt
    Zopt::TZopt
end

function SparseVarLatent(
    T::DataType,
    dim::Int,
    S::Int,
    Z::AbstractVector,
    kernel::Kernel,
    mean::PriorMean,
    opt=nothing,
    Zopt=nothing
)
    return SparseVarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean), cholesky(Matrix{T}(I(dim)))),
        VarPosterior{T}(dim),
        deepcopy(Z),
        Matrix{T}(undef, S, dim),
        Matrix{T}(undef, S, dim),
        Vector{T}(undef, S),
        deepcopy(opt),
        deepcopy(Zopt),
    )
end

@traitimpl IsSparse{SparseVarLatent}

## Monte-Carlo Gaussian Process

struct SampledLatent{T,Tpr<:GPPrior,Tpo<:SampledPosterior{T}} <: AbstractLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
end

function SampledLatent(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean)
    return SampledLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean), cholesky(Matrix{T}(I, dim, dim))),
        SampledPosterior(dim, zeros(T, dim), Symmetric(Matrix{T}(I(dim)))),
    )
end

@traitimpl IsFull{SampledLatent}

## Online Sparse Variational Process

mutable struct OnlineVarLatent{T,Tpr<:GPPrior,Tpo<:AbstractVarPosterior{T},Topt,TZ<:AbstractVector,
    TZalg<:InducingPoints.OnIPSA,TZopt} <:
               AbstractVarLatent{T,Tpo,Tpr}
    prior::Tpr
    post::Tpo
    Z::TZ
    Zalg::TZalg
    Knm::Matrix{T}
    κ::Matrix{T}
    K̃::Vector{T}
    Zupdated::Bool
    opt::Topt
    Zopt::TZopt
    Zₐ::AbstractVector
    Kab::Matrix{T}
    κₐ::Matrix{T}
    K̃ₐ::Matrix{T}
    invDₐ::Symmetric{T,Matrix{T}}
    prev𝓛ₐ::T
    prevη₁::Vector{T}
end

function OnlineVarLatent(
    T::DataType,
    dim::Int,
    nSamplesUsed::Int,
    Z::AbstractVector,
    Zalg::InducingPoints.OnIPSA,
    kernel::Kernel,
    mean::PriorMean,
    opt=nothing,
    Zopt=nothing
)
    return OnlineVarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean), cholesky(Matrix{T}(I, dim, dim))),
        OnlineVarPosterior{T}(dim),
        Z,
        Zalg,
        Matrix{T}(undef, nSamplesUsed, dim),
        Matrix{T}(undef, nSamplesUsed, dim),
        Vector{T}(undef, nSamplesUsed),
        false,
        deepcopy(opt),
        deepcopy(Zopt),
        deepcopy(Z),
        Matrix{T}(I, dim, dim),
        Matrix{T}(I, dim, dim),
        Matrix{T}(I, dim, dim),
        Symmetric(Matrix{T}(I, dim, dim)),
        zero(T),
        Vector{T}(undef, dim),
    )
end

@traitimpl IsSparse{OnlineVarLatent}

## Variational Student-T Process

mutable struct TVarLatent{T<:Real,Tpr<:TPrior,Tpo<:VarPosterior{T},O} <:
               AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function TVarLatent(T::DataType, ν::Real, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return TVarLatent(
        TPrior(
            deepcopy(kernel),
            deepcopy(mean),
            cholesky(Matrix{T}(I, dim, dim)),
            ν,
            rand(T),
            rand(T),
        ),
        VarPosterior{T}(dim),
        deepcopy(opt),
    )
end

@traitimpl IsFull{TVarLatent}

### Functions

prior(gp::AbstractLatent) = gp.prior
kernel(gp::AbstractLatent) = kernel(prior(gp))
setkernel!(gp::AbstractLatent, kernel::Kernel) = setkernel!(prior(gp), kernel)
pr_mean(gp::AbstractLatent) = mean(prior(gp))
pr_mean(gp::AbstractLatent, X::AbstractVector) = mean(prior(gp), X)
setpr_mean!(gp::AbstractLatent, μ₀::PriorMean) = setmean!(prior(gp), μ₀)
pr_cov(gp::AbstractLatent) = cov(prior(gp))
pr_cov(gp::TVarLatent) = prior(gp).χ * cov(prior(gp))
pr_cov!(gp::AbstractLatent, K::Cholesky) = gp.prior.K = K

posterior(gp::AbstractLatent) = gp.post
Distributions.dim(gp::AbstractLatent) = dim(posterior(gp))
Distributions.mean(gp::AbstractLatent) = mean(posterior(gp))
Distributions.cov(gp::AbstractLatent) = cov(posterior(gp))
Distributions.var(gp::AbstractLatent) = var(posterior(gp))
nat1(gp::AbstractVarLatent) = nat1(posterior(gp))
nat2(gp::AbstractVarLatent) = nat2(posterior(gp))

mean_f(model::AbstractGP) = mean_f.(model.f)

@traitfn mean_f(gp::T) where {T <: AbstractLatent; IsFull{T}} = mean_f(mean(gp))
@traitfn mean_f(gp::T) where {T <: AbstractLatent; !IsFull{T}} = mean_f(mean(gp), gp.κ)

mean_f(μ::AbstractVector) = μ
mean_f(μ::AbstractVector, κ::AbstractMatrix) = κ * μ

var_f(model::AbstractGP) = var_f.(model.f)

@traitfn var_f(gp::T) where {T <: AbstractLatent; IsFull{T}} = var_f(cov(gp))
@traitfn var_f(gp::T) where {T <: AbstractLatent; !IsFull{T}} = var_f(cov(gp), gp.κ, gp.K̃)

var_f(Σ::AbstractMatrix) = diag(Σ)
var_f(Σ::AbstractMatrix, κ::AbstractMatrix, K̃::AbstractVector) = diag_ABt(κ * Σ, κ) + K̃

Zview(gp::SparseVarLatent) = gp.Z
Zview(gp::OnlineVarLatent) = gp.Z

setZ!(gp::AbstractLatent, Z::AbstractVector) = gp.Z = Z

opt(gp::AbstractLatent) = gp.opt
Zopt(::AbstractLatent) = nothing
Zopt(gp::SparseVarLatent) = gp.Zopt
Zopt(gp::OnlineVarLatent) = gp.Zopt

@traitfn function compute_K!(
    gp::TGP, X::AbstractVector, jitt::Real
) where {TGP <: AbstractLatent; IsFull{TGP}}
    return pr_cov!(gp, cholesky(kernelmatrix(kernel(gp), X) + jitt * I))
end

@traitfn function compute_K!(gp::T, jitt::Real) where {T <: AbstractLatent; !IsFull{T}}
    return pr_cov!(gp, cholesky(kernelmatrix(kernel(gp), gp.Z) + jitt * I))
end

function compute_κ!(gp::SparseVarLatent, X::AbstractVector, jitt::Real)
    gp.Knm = kernelmatrix(kernel(gp), X, gp.Z)
    gp.κ = copy(gp.Knm / pr_cov(gp))
    gp.K̃ = kernelmatrix_diag(kernel(gp), X) .+ jitt - diag_ABt(gp.κ, gp.Knm)
    return all(gp.K̃ .> 0) || error("K̃ has negative values")
end

function compute_κ!(gp::OnlineVarLatent, X::AbstractVector, jitt::Real)
    # Covariance with the model at t-1
    gp.Kab = kernelmatrix(kernel(gp), gp.Zₐ, gp.Z)
    gp.κₐ = gp.Kab / pr_cov(gp)
    Kₐ = Symmetric(kernelmatrix(kernel(gp), gp.Zₐ) + jitt * I)
    gp.K̃ₐ = Kₐ - gp.κₐ * transpose(gp.Kab)

    # Covariance with a new batch
    gp.Knm = kernelmatrix(kernel(gp), X, gp.Z)
    gp.κ = gp.Knm / pr_cov(gp)
    gp.K̃ = kernelmatrix_diag(kernel(gp), X) .+ jitt - diag_ABt(gp.κ, gp.Knm)
    all(gp.K̃ .> 0) || error("K̃ has negative values")
end


text:
export dualcgp_config, monocgp_config

using Dates

function dict_to_namedtuple(d::Dict)
    NamedTuple{Tuple(keys(d))}(values(d))
end

function dualcgp_config(dualcgp_cfg_filename::String, game_name::String)
    cfg = YAML.load_file(dualcgp_cfg_filename)
    dualcgp_config(cfg, game_name)
end

function dualcgp_config(cfg::NamedTuple, game_name::String)
    cfg_dict = Dict(pairs(cfg))
    dualcgp_config(cfg_dict, game_name)
end

"""
    dualcgp_config(cfg::Dict, game_name::String)

Retrieve the encoder and controller configuration files from the main dual CGP
configuration file.
Return both config dictionaries along with the corresponding reducer.
"""
function dualcgp_config(cfg::Dict, game_name::String)
    # Temporarily open a game to retrieve parameters
    seed = cfg["seed"]
    game = Game(game_name, seed)
    grayscale = cfg["grayscale"]
    downscale = cfg["downscale"]
    s = get_state(game, grayscale, downscale)
    img_size = size(s[1])
    n_in = length(s)
    n_out = length(getMinimalActionSet(game.ale))  # One output per legal action
    close!(game)

    # Main config and  initialize sub-configs
    bootstrap = cfg["bootstrap"]
    d_fitness = cfg["d_fitness"]
    n_gen = cfg["n_gen"]
    log_gen = cfg["log_gen"]
    save_gen = cfg["save_gen"]
    encoder_cfg = cfg["encoder"]
    reducer_cfg = cfg["reducer"]
    controller_cfg = cfg["controller"]
    reducer_type = reducer_cfg["type"]
    logid = string(Dates.now(), "_", game_name)
    for k in ["seed", "d_fitness", "n_gen", "log_gen", "save_gen"]
        encoder_cfg[k] = cfg[k]
        controller_cfg[k] = cfg[k]
    end

    # Main config
    main_cfg = Dict()
    for k in keys(cfg)
        if k ∉ ["encoder", "reducer", "controller"]
            main_cfg[k] = cfg[k]
        end
    end
    main_cfg["id"] = logid
    # main_cfg = dict_to_namedtuple(main_cfg)

    # Encoder config
    encoder_cfg["function_module"] = IICGP.CGPFunctions
    encoder_cfg["n_in"] = n_in
    encoder_cfg["img_size"] = img_size
    encoder_cfg["id"] = logid
    encoder_cfg["type"] = "encoder"
    encoder_cfg = get_config(encoder_cfg) # dict to named tuple

    # Reducer
    reducer = Reducer(reducer_cfg, n_in=encoder_cfg.n_out, img_size=img_size)

    # Forward pass to retrieve the number of input of the controller
    enco = IPCGPInd(encoder_cfg)
    enco_out = CartesianGeneticProgramming.process(enco, s)
    features = reducer.reduct(enco_out, reducer.parameters)
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    if !("n_cst_inputs" in keys(controller_cfg))
        controller_cfg["n_cst_inputs"] = 0
    end
    cont_n_in = length(features_flatten) + controller_cfg["n_cst_inputs"]

    # Controller config
    controller_cfg["function_module"] = IICGP.CGPFunctions
    controller_cfg["n_in"] = cont_n_in
    controller_cfg["n_out"] = n_out
    controller_cfg["id"] = logid
    controller_cfg["type"] = "controller"
    controller_cfg = get_config(controller_cfg) # dict to named tuple

    main_cfg, encoder_cfg, controller_cfg, reducer, bootstrap
end

function monocgp_config(monocgp_cfg_filename::String, game_name::String)
    cfg = YAML.load_file(monocgp_cfg_filename)
    monocgp_config(cfg, game_name)
end

"""
    monocgp_config(cfg::Dict, game_name::String)

Retrieve the controller configuration file from the main mono CGP configuration
file.
Return the config dictionary along with the corresponding reducer.
"""
function monocgp_config(cfg::Dict, game_name::String)
    # Temporarily open a game to retrieve parameters
    seed = cfg["seed"]
    game = Game(game_name, seed)
    grayscale = cfg["grayscale"]
    downscale = cfg["downscale"]
    s = get_state(game, grayscale, downscale)
    img_size = size(s[1])
    n_in = length(s)
    n_out = length(getMinimalActionSet(game.ale))  # One output per legal action
    close!(game)

    # Main config and  initialize sub-configs
    bootstrap = cfg["bootstrap"]
    d_fitness = cfg["d_fitness"]
    n_gen = cfg["n_gen"]
    log_gen = cfg["log_gen"]
    save_gen = cfg["save_gen"]
    reducer_cfg = cfg["reducer"]
    controller_cfg = cfg["controller"]
    reducer_type = reducer_cfg["type"]
    logid = string(Dates.now(), "_", game_name)
    for k in ["seed", "d_fitness", "n_gen", "log_gen", "save_gen"]
        controller_cfg[k] = cfg[k]
    end

    # Main config
    main_cfg = Dict()
    for k in keys(cfg)
        if k ∉ ["encoder", "reducer", "controller"]
            main_cfg[k] = cfg[k]
        end
    end
    # main_cfg = dict_to_namedtuple(main_cfg)

    # Reducer
    reducer = Reducer(reducer_cfg, n_in=n_in, img_size=img_size)

    # Forward pass to retrieve the number of input of the controller
    features = reducer.reduct(s, reducer.parameters)
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    cont_n_in = length(features_flatten)

    # Controller config
    controller_cfg["function_module"] = IICGP.CGPFunctions
    controller_cfg["n_in"] = cont_n_in
    controller_cfg["n_out"] = n_out
    controller_cfg["id"] = logid
    controller_cfg = get_config(controller_cfg) # dict to named tuple

    main_cfg, controller_cfg, reducer, bootstrap
end


text:
include("DeterministicInternal.jl")
include("Deterministic.jl")
include("ProbabilisticInternal.jl")
include("Probabilistic.jl")
include("ScenarioBasedInternal.jl")
include("ScenarioBased.jl")
include("PiecewiseFunction.jl")
include("PiecewiseFunctionInternal.jl")

export get_break_points
export get_data
export get_horizon
export get_initial_time
export get_internal
export get_label
export get_percentiles
export get_resolution
export get_scenario_count
export get_time_series_uuid
export set_break_points!
export set_data!
export set_horizon!
export set_initial_time!
export set_internal!
export set_label!
export set_percentiles!
export set_resolution!
export set_scenario_count!
export set_time_series_uuid!


text:
using Plots, BlackBoxOptim, Random, Dates, LinearAlgebra, Statistics, Printf, Optim, OptimTestProblems
plot(rand(10),xticks=1:1:10)

curr_path=pwd()
include(string(curr_path,"\\functions_opti.jl"))
include(string(curr_path,"\\call_bbox.jl"))
# ibb=[4 9 11] #
ibb=[1 2 4 5 6 7 9 10 11 12 13 14 15 16] #

funcs1=[:x_i, :elliptic, :cigar, :cigtab, :x_5_sq, :sin_01_x_2, :griewank, :quartic, :schwefel1_2, :rastrigin, :sphere, :ellipsoid, :alpine_1]
funcs1=[:sin_01_x_2]
# funcs1=[funcs1;funcs1;funcs1;funcs1;funcs1]

# include(string(curr_path,"\\dependencies_pdf_cdf.jl")) # για cdf
include(string(curr_path,"\\dependencies.jl")) # merged
vars1=Int64(20);evals1=Int64(10000)
history_all,times_all_step,x_evals,opti_evals,cdf_y,inp_yy,inp_x,pdf_y,iter,i_true=run_OHT(vars1,evals1,ibb,funcs1);
sizex=700; sizey=400; sf1=2.7
include(string(curr_path,"\\plots_all.jl"))

history_all_all=history_all.*0.0
nof_runs=10
for i=1:nof_runs
  global ibb,funcs1,history_all,vars1,evals1
  println("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<",i,">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>")
  history_all,times_all_step,x_evals,opti_evals,cdf_y,inp_yy,inp_x,pdf_y,iter,i_true=run_OHT(vars1,evals1,ibb,funcs1);
  include(string(curr_path,"\\plots_all.jl"))
  history_all_all.+=history_all
end
history_all=history_all_all./nof_runs
sizex=500; sizey=400; sf1=2.7
include(string(curr_path,"\\plots_all.jl"))

savefig(string(curr_path,"\\FIGURES\\history\\all_",nof_runs,"times_",vars1,"vars_",evals1,"evals.pdf"))


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -42.99754755332367 <= q <= 40.38858433100043)

                     add_NL_constraint(m, :(log(1 + exp(-0.8925450194703983 + 0.6110890427800215*log(1 + exp(-0.16626426792027038 + 0.44845618413236066*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.6799320945085388*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7514544713734312*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.795720190818876*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + 0.1740395626916258*log(1 + exp(0.8063063231762757 + 0.4905142959916642*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.7092497630992938*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.8897979792487347*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.11224890387090136*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.8017072131372296*log(1 + exp(0.9210510999614945 + -0.4935968211971775*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.33547745483400737*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.1703969445706579*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.2742027660464488*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.7966067754659663*log(1 + exp(-0.7328670881809516 + 0.6102429954796835*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.6353392761755483*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7830650903096812*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.1627621376749886*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + log(1 + exp(-0.7008872510227278 + -0.258173144464545*log(1 + exp(-0.16626426792027038 + 0.44845618413236066*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.6799320945085388*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7514544713734312*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.795720190818876*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + 0.49550787111793193*log(1 + exp(0.8063063231762757 + 0.4905142959916642*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.7092497630992938*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.8897979792487347*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.11224890387090136*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.8462943069272191*log(1 + exp(0.9210510999614945 + -0.4935968211971775*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.33547745483400737*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.1703969445706579*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.2742027660464488*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.05869131996649157*log(1 + exp(-0.7328670881809516 + 0.6102429954796835*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.6353392761755483*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7830650903096812*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.1627621376749886*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + log(1 + exp(-0.2643911691309153 + -0.31660511055293883*log(1 + exp(-0.16626426792027038 + 0.44845618413236066*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.6799320945085388*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7514544713734312*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.795720190818876*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.09136517532701127*log(1 + exp(0.8063063231762757 + 0.4905142959916642*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.7092497630992938*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.8897979792487347*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.11224890387090136*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.36400896662330684*log(1 + exp(0.9210510999614945 + -0.4935968211971775*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.33547745483400737*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.1703969445706579*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.2742027660464488*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + 0.4876638787119827*log(1 + exp(-0.7328670881809516 + 0.6102429954796835*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.6353392761755483*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7830650903096812*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.1627621376749886*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) + log(1 + exp(0.6679332747384503 + -0.6061375406574494*log(1 + exp(-0.16626426792027038 + 0.44845618413236066*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.6799320945085388*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7514544713734312*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.795720190818876*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.7532640510772808*log(1 + exp(0.8063063231762757 + 0.4905142959916642*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + 0.7092497630992938*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.8897979792487347*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.11224890387090136*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + 0.9959397069462161*log(1 + exp(0.9210510999614945 + -0.4935968211971775*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.33547745483400737*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + 0.1703969445706579*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + -0.2742027660464488*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))) + -0.8544693505755525*log(1 + exp(-0.7328670881809516 + 0.6102429954796835*log(1 + exp(-0.17340316897790675 + 0.7538424045603822*$(x[1]) + -0.010659538195954*$(x[2]) + 0.3269307983949328*$(x[3]) + 0.4439680182837229*$(x[4]) + 0.638115725320251*$(x[5]))) + -0.6353392761755483*log(1 + exp(-0.2603616500556529 + 0.5967362687072955*$(x[1]) + -0.6976319849871899*$(x[2]) + -0.3505130257519107*$(x[3]) + 0.7751645636807667*$(x[4]) + 0.9801389530891269*$(x[5]))) + -0.7830650903096812*log(1 + exp(-0.6109406919333629 + 0.2096051265369132*$(x[1]) + -0.41861538796679865*$(x[2]) + -0.21379191209843995*$(x[3]) + 0.8880999966792658*$(x[4]) + 0.29816559910815243*$(x[5]))) + 0.1627621376749886*log(1 + exp(0.4223539181621372 + 0.05773940633058361*$(x[1]) + -0.7624642273283109*$(x[2]) + -0.13780284117445118*$(x[3]) + -0.6723420679434677*$(x[4]) + 0.47304736250148327*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
######### Stan program example  ###########

using Mamba, Stan

old = pwd()
ProjDir = normpath(homedir(), "Projects", "Julia", "Rob", "Support_Stan", "ODE")
cd(ProjDir)

odemodel = "
functions {
       real[] sho(real t,
                  real[] y,
                  real[] theta,
                  real[] x_r,
                  int[] x_i) {
         real dydt[2];
         dydt[1] <- y[2];
         dydt[2] <- -y[1] - theta[1] * y[2];
         return dydt;
        }
}

data {
  int<lower=1> T;
  int<lower=1> M;
  real t0;
  #real y0[M];
  real ts[T];
  real y[T,M];
}

transformed data {
  real x_r[0];
  int x_i[0];
}

parameters {
  real y0[2];
  vector<lower=0>[2] sigma;
  real theta[1];
}

model {
  real y_hat[T,2];
  sigma ~ cauchy(0,2.5);
  theta ~ normal(0,1);
  y0 ~ normal(0,1);
  y_hat <- integrate_ode(sho, y0, t0, ts, theta, x_r, x_i);
  for (t in 1:T)
   y[t] ~ normal(y_hat[t], sigma);
}
"

odedict =
  Dict(
    "T" => 100,
    "M" => 2,
    "t0" => 0,
    "y0" =>[1, 0],
    "ts" => range(0.1, stop=10, length=100),
    "y" => reshape(
    [1.18447214273212, 0.931328995894153, 1.15546788943838,
    1.11347639821441, 0.942770215168461, 0.599435493714444, 0.633508944792969,
    0.664145932981416, 0.636983924087869, 0.922765008328693, 0.596380019171069,
    0.278177612651595, 0.139369448614547, 0.179787763946766, 0.0891127323154026,
    -0.0169426249192601, -0.00570968277658015, -0.263912169092685,
    -0.148772505096642, -0.480382763289494, -0.404896182791746, -0.386328653102745,
    -0.489160545877774, -0.448720573538264, -0.631560016157571, -0.594120804164941,
    -0.546128595678147, -0.844353075954525, -0.957296903080828, -0.773557074171747,
    -0.823902983465925, -0.870022238251538, -0.845805547618517, -0.862822747019435,
    -0.633621387892182, -0.54005323688692, -0.52759594518313, -0.698260966423782,
    -0.399851384871993, -0.57417194227814, -0.210773536528394, -0.328691086366792,
    -0.415888595882974, -0.405460494417675, -0.386540582062706, -0.301225587058244,
    -0.305429707660651, 0.172795987572527, 0.193262969285858, 0.101787696286287,
    0.240887604250816, 0.206580592212456, 0.687874694334188, 0.256891132411694,
    0.417940378943263, 0.508702198428123, 0.603583494034937, 0.518963474603121,
    0.239370958898509, 0.405413842254413, 0.664537334395703, 0.61855117823015,
    0.482298625607209, 0.602981201676457, 0.488960693669512, 0.632472522128953,
    0.361607240790125, 0.604105168457383, 0.555260897330889, 0.491905227250743,
    0.444851458749041, 0.43665875657333, 0.253492488343665, 0.282978685637294,
    0.348245354984117, 0.359853365329113, 0.161227722113711, 0.0663877602751073,
    -0.108375655042251, -0.242273836124497, -0.200186428760795, 0.0551064659674957,
    -0.0674686721293205, -0.263935344013016, -0.339446542506593, -0.379888014913294,
    -0.205050970863634, -0.429446704363204, -0.228811177404085, -0.344775841425698,
    -0.266696154615485, -0.607526316761665, -0.485276356000595, -0.623683758487882,
    -0.402186147508401, -0.468862353806569, -0.519585081598368, -0.244708739840049,
    -0.410221061395796, -0.272425410718381, 0.0181904776900387, -0.312235493264689,
    -0.381368486156188, -0.37097995457925, -0.631447866640104, -0.453475182715527,
    -0.803672414768249, -0.432186904480258, -0.807890489400406, -0.529710418827556,
    -0.88352798312236, -0.998941025325834, -0.871810569861829, -0.884996699488221,
    -1.14562740835939, -0.731091370545707, -1.04389586568583, -0.803924418504678,
    -0.750210634933026, -0.76615209369324, -0.759865803866824, -0.660889360648131,
    -0.793258906688218, -0.690430982604711, -0.66929377044793, -0.193522247444425,
    -0.130968185844238, -0.238937181200586, -0.38435217993569, -0.120257530068827,
    0.186771342421634, -0.0324051426644301, 0.236004465629662, 0.045278615833121,
    0.366831904664821, 0.18841113450848, 0.211083941237886, 0.311481217801157,
    0.378327568947026, 0.420127897626007, 0.709620445789031, 0.643849729991572,
    0.613348386811806, 0.520286542950945, 0.545129657971766, 0.589478472290371,
    0.536722990830953, 0.630499864509491, 0.754715456167985, 0.486532397567358,
    0.857517549562448, 0.407717911256268, 0.551895261585553, 0.878724986192702,
    0.609205971465613, 0.386122525791358, 0.238987436773361, 0.398304899982967,
    0.248882597950643, 0.132786785682916, 0.173995860263871, -0.0100293452254933,
    0.39937911504977, 0.190839430703432, -0.004520398638622, -0.0726269219830769,
    -0.31655229056679, -0.15460224391486, -0.388478202676613, -0.532056186725431,
    -0.138588618642611, -0.595854661880176, -0.531027468041037, -0.604247540619453,
    -0.673722368761153, -0.455905727528453, -0.782762232356633, -0.587995350317593,
    -0.640490312145211, -0.754692905495262, -0.543258607376073, -0.578773929440613,
    -0.579598805944238, -0.315701576014008, -0.350936596499424, -0.408500418636149,
    -0.332210426245586, -0.424089753625995, -0.299022180880827, -0.0568205107389459,
    -0.176479238422781, -0.100509704966414, 0.22863110851166, -0.130938638084811,
    0.167944957592579, 0.340641058030672, -0.0407477745000139, 0.16645548196369,
    0.147962705604123, 0.172782087145017], 100, 2))

odedata = [
  odedict,
  odedict,
  odedict,
  odedict
]

stanmodel = Stanmodel(name="ode", model=odemodel, nchains=4);

println("\nStanmodel that will be used:")
stanmodel |> display
println("Input observed data dictionary:")
odedata |> display
println()

sim1 = stan(stanmodel, odedata, ProjDir, diagnostics=false, CmdStanDir=CMDSTAN_HOME);

## Subset Sampler Output to variables suitable for describe().
monitor = ["theta.1", "lp__", "accept_stat__"]
sim = sim1[:, monitor, :]
describe(sim)
println()

p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

## Subset Sampler Output to variables suitable for describe().
monitor = ["y0.1", "y0.2", "sigma.1", "sigma.2"]
sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
describe(sim)
println()

p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-sigmaplot", fmt=:svg)
draw(p, nrow=4, ncol=4, filename="$(stanmodel.name)-sigmaplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static Sys.isapple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-summaryplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-summaryplot-$(i).svg"`)
          isfile("$(stanmodel.name)-sigmaplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-sigmaplot-$(i).svg"`)
        end
      end : println()

cd(old)


text:
"""
Subroutine insort

Given an integer array keys of length n, this subroutine uses
an insertion sort to sort the keys in increasing order.

MINPACK-2 Project. March 1998.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function insort(n,keys)
    @inbounds for j=2:n
        ind = keys[j]
        i = j - 1
        while (i > 0) && (keys[i] > ind)
            keys[i+1] = keys[i]
            i = i - 1
        end
        keys[i+1] = ind
    end

    return
end


text:
using RungeKutta
using Test

include("test_utils.jl")
include("test_tableau.jl")
include("test_tableau_partitioned.jl")
include("test_order_conditions.jl")
include("test_symmetry.jl")
include("test_symplecticity.jl")
include("test_gauss.jl")
include("test_lobatto.jl")
include("test_radau.jl")
include("test_tableaus.jl")
include("test_tableaus_prk.jl")
include("test_list.jl")


text:
# AbelianSpace.jl
#
# Defines the immutable AbelianSpace for a complex vector space that is graded by
# the unitary irreps of an Abelian group.

# AbelianSpace
#--------------
immutable AbelianSpace{G<:Abelian} <: UnitaryRepresentationSpace{G}
    dims::Dict{G,Int}
    dual::Bool
end
AbelianSpace{G<:Abelian}(dims::Dict{G,Int})=AbelianSpace{G}(dims, false)

# Corresponding methods:
sectors(V::AbelianSpace) = keys(V.dims)

dim(V::AbelianSpace) = sum(values(V.dims))::Int
dim{G<:Abelian}(V::AbelianSpace{G},c::G) = get(V.dims,c,0)
dual(V::AbelianSpace) = AbelianSpace([conj(c)=>dim(V,c) for c in sectors(V)], !V.dual)
cnumber{G}(V::AbelianSpace{G}) = AbelianSpace(Dict{G,Int}(one(G)=>1),V.dual)
cnumber{G}(::Type{AbelianSpace{G}}) = AbelianSpace(Dict{G,Int}(one(G)=>1))
iscnumber{G}(V::AbelianSpace{G}) = V.dims == Dict{G,Int}(one(G)=>1)

# Comparison
==(V1::AbelianSpace,V2::AbelianSpace)=(V1.dims==V2.dims && V1.dual==V2.dual)

# Show methods
Base.show(io::IO, V::AbelianSpace) = (print(io,"AbelianSpace");showcompact(io,V.dims);print(io, V.dual ? "*" : ""))

# direct sum of AbelianSpaces
function directsum{G}(V1::AbelianSpace{G}, V2::AbelianSpace{G})
    V1.dual == V2.dual || throw(SpaceError("Direct sum of a vector space and its dual do not exist"))
    dims=Dict{G,Int}
    for c in union(sectors(V1),sectors(V2))
        dims[c]=dim(V1,c)+dim(V2,c)
    end
    return AbelianSpace(dims,V1.dual)
end

# fusing and splitting AbelianSpaces
function fuse{G}(V1::AbelianSpace{G}, V2::AbelianSpace{G}, V::AbelianSpace{G})
    dims=[c1*c2=>0 for c1 in sectors(V1), c2 in sectors(V2)]
    for c1 in sectors(V1), c2 in sectors(V2)
        dims[c1*c2]+= dim(V1,c1)*dim(V2,c2)
    end
    return V.dims == dims
end

# indexing using sectors
function to_range{G}(s::G,V::AbelianSpace{G})
    offset=0
    for c in sort!(collect(sectors(V)),rev=V.dual,by=s->s.charge)
        if c!=s
            offset+=dim(V,c)
        else
            break
        end
    end
    return offset+(1:dim(V,s))
end


text:
# # Lovász Theta Function
#
# The Lovász theta function is an important concept in graph theory. Consider an undirected graph $G = (V,E)$ with vertex set $V = \{1,\dots,n\}$ and edge set $E$ with  $E_{ij} = 1$ if there exists an edge between the vertices $i$ and $j$. A stable set (or independent set) is a subset of $V$ such that the induced subgraph does not contain edges. The *stability number* $\alpha(G)$ of the graph is equal to the cardinality of the largest stable set. It is closely related to the *Shannon capacity* $\Theta(G)$ that models the amount of information that a noisy communication channel can carry if certain signal values can be confused with each other. Unfortunately, the determination of $\alpha(G)$ is an NP-hard problem and the computational complexity to determine $\Theta(G)$ remains unknown. The Lovász theta function $\vartheta(G)$ was introduced in \[1\], can be computed in polynomial time, and represents an upper bound on both the stability number and the Shannon capacity:
# $$
# \alpha(G) \leq \Theta(G) \leq \vartheta(G).
# $$
# The value of $\vartheta(G)$ can be determined by computing the optimal value $p^*$ of the following SDP:
# $$
# \begin{array}{ll} \text{maximize} &   \text{Tr}(JX)\\
# \text{subject to} &  \text{Tr}(X) = 1, \\
#                   &  X_{ij} = 0, \quad (i, j) \in E, \\
#                   &  X \succeq 0,
# \end{array}
# $$
# with matrix variable $X$, matrix $J \in \mathbf{R}^{n \times n}$ of all ones and $\text{Tr}()$ denoting the matrix trace.
# In this simple example we will compute the value of the Lovász theta function for the *Petersen Graph* (which is known to be $4$).

#-
using LinearAlgebra, COSMO, JuMP, Plots, GraphRecipes

## let's define the Petersen graph
n = 10
E = zeros(n, n)
E[1, 2] = E[1, 5] = E[1, 6] = 1.
E[2, 3] = E[2, 7]  = 1.
E[3, 4] = E[3, 8]  = 1.
E[4, 5] = E[4, 9] = 1.
E[5, 10] = 1.
E[6, 8] = E[6, 9] = 1.
E[7, 9] = E[7, 10] = 1.
E[8, 10] = 1.

## plot the graph
ri = 1.
ro = 2.
coordinates = []
for θ = 90:-72:-198
  push!(coordinates, [ro * cosd(θ), ro * sind(θ)])
end
for θ = 90:-72:-198
  push!(coordinates, [ri * cosd(θ), ri * sind(θ)])
end
graphplot(E, names = 1:n,  x = getindex.(coordinates, 1), y = getindex.(coordinates, 2), fontsize = 10, nodesize = 1, nodeshape =:circle, curvature = 0.)

# Let's solve the SDP with COSMO and JuMP:
#-
model = JuMP.Model(COSMO.Optimizer);

@variable(model, X[1:n, 1:n], PSD)
x = vec(X)
@objective(model, Max, sum(x))
@constraint(model, tr(X)== 1.)
for j = 1:n
  for i = 1:j-1
    if E[i, j] == 1.
      @constraint(model, X[i, j] == 0)
    end
  end
end
status = JuMP.optimize!(model)

# The optimal objective is given by:
#-
JuMP.objective_value(model)
# Which is the correct known value for the Petersen Graph.
#
# ## References
# [1] Lovász - On the Shannon Capacity of a Graph, IEEE Transactions on Information Theory (1979)


text:
__precompile__()
module ObjectTracker

export
    Blob,
    form_blobs,
    is_transient,
    match!,
    Object,
    predict,
    State,
    update

include("Kalman.jl")
include("Blobs.jl")
include("Tracker.jl")
end # module


text:
# Author: Lukasz Bratos

# Funcja oblicza wartość wyrażenia podanego w zadaniu
# n - liczba iteracji, p0 - wartość początkowa p0, r - dana stała
function func(n, p0, r)
    pn = p0
    for i in 1:n
        pn = pn + r * pn * (1 - pn)
    end
    return pn
end

# podpunkt 1

println("Funkcja obliczana normalnie: ", func(40, Float32(0.01), Float32(3.0)))
println("Funkcja obliczana z zaokrągleniem: ", func(30, Float32(0.722), Float32(3.0)))

# podpunkt 2

println("Wynik dla Float32: ", func(40, Float32(0.01), Float32(3.0)))
println("Wynik dla Float64: ", func(40, Float64(0.01), Float64(3.0)))


text:
export make_grid_min_max, BMtransition_from_data

function make_grid_min_max(logx::Vector; n::Int=31, scal::Real=5.0)
    a = minimum(logx)-log(scal)
    b = maximum(logx)+log(scal)
    step = (b-a)/(n-1)
    return a : step : b
end

markov_transition_match_option(μ::Function, Σ::Function, ::Type{Val{true}}, minp::Real,  grid::AbstractVector)     = markov_transition_moment_matching_serial(  μ, Σ, minp, grid)
markov_transition_match_option(μ::Function, Σ::Function, ::Type{Val{true}}, minp::Real,  grids::AbstractVector...) = markov_transition_moment_matching_parallel(μ, Σ, minp, grids...)
markov_transition_match_option(μ::Function, Σ::Function, ::Type{Val{false}}, minp::Real, grids::AbstractVector...) = markov_transition(μ, Σ, minp, grids...)

function BMtransition_from_data(logxy::Array{T,2}, matching::Type; zero_cov::Bool=true, minp::Real=1e-9, kwargs...) where {T<:Real}
    0.0 <= minp < 1.0 || throw(DomainError())
    grids = ((make_grid_min_max(logxy[:,j]; kwargs...) for j in 1:size(logxy, 2))...,)
    Σ = cov(diff(logxy,1))
    zero_cov  &&  Σ .= Diagonal(Σ)
    transition = markov_transition_match_option((xy::NTuple) -> [xy...], (xy::NTuple) -> Σ, matching, minp, grids...)
    return (grids, Σ, transition, )
end

function BMtransition_from_data(logx::Array{T,1}, matching::Type; minp::Real=1e-9, kwargs...) where {T<:Real}
    0.0 <= minp < 1.0 || throw(DomainError())
    grid = make_grid_min_max(logx; kwargs...)
    σ = sqrt(var(diff(logx)))
    transition = markov_transition_match_option( (x) -> x[1], (x) -> σ, matching, minp, grid)
    return (grid, σ^2., transition, )
end


text:
function  kerker_mix(res, mixOpts::kerkerMixOptions)
#The Kerker mixing, following the parameter setup of elliptic
#preconditioner.
#
#Lin Lin
#Last revision: 5/6/2012

return kerker_mix(res, mixOpts.KerkerB, mixOpts.kmult,
                       mixOpts.YukawaK, mixOpts.epsil0)

end

function  kerker_mix(res, KerkerB::Float64, kmult, YukawaK::Float64, epsil0::Float64)
#The Kerker mixing, following the parameter setup of elliptic
#preconditioner.
#
#Lin Lin
#Last revision: 5/6/2012

resfft = fft(res);

resfft = ( epsil0 / (4*pi) * (kmult .+ YukawaK^2) )./
  ( KerkerB .+ epsil0 / (4*pi) * (kmult .+ YukawaK^2) ) .* resfft;

resprec = real(ifft(resfft));
resprec = resprec .- mean(resprec);

return resprec

end

max_stars_repo_name: Forgotten/modelscf

max_issues_repo_name: Forgotten/modelscf

max_forks_repo_name: Forgotten/modelscf


text:
using Replay

instructions = """
using FileIO, Sixel, Plots
gr()
buf = IOBuffer()
show(buf, MIME("image/png"), plot(sin, size=(1000, 750)))
buf |> load |> sixel_encode
"""

replay(instructions, julia_project=@__DIR__)


text:
#include("Misc/mag.jl")
#include("Misc/arf.jl")
#include("Misc/arb.jl")
#include("Misc/acb.jl")
#include("Misc/acb_poly2.jl")
include("Misc/acb_root_ctx.jl")
include("Misc/arb_mat2.jl")
include("Misc/fmpr.jl")