text:
# some basic functionality testing
x = 1234.56789

@test fmt(x) == "1234.567890"
@test fmt(x;prec=2) == "1234.57"
@test fmt(x,10,3) == "  1234.568"
@test fmt(x,10,3,:left) == "1234.568  "
@test fmt(x,10,3,:ljust) == "1234.568  "
@test fmt(x,10,3,:right) == "  1234.568"
@test fmt(x,10,3,:lrjust) == "  1234.568"
@test fmt(x,10,3,:zpad) == "001234.568"
@test fmt(x,10,3,:zeropad) == "001234.568"
@test fmt(x,:commas) == "1,234.567890"
@test fmt(x,10,3,:left,:commas) == "1,234.568 "
@test fmt(x,:ipre) == "1234.567890"
@test fmt(x,12) == " 1234.567890"

i = 1234567

@test fmt(i) == "1234567"
@test fmt(i,:commas) == "1,234,567"

@test fmt(2 - 3im, 10) == "   2 - 3im"
@test fmt(pi - 3im, 15, 2) == "  3.14 - 3.00im"

@test fmt(3//4, 10) == "      3//4"
@test fmt(1//2 + 6//2 * im, 15) == " 1//2 + 3//1*im"

fmt_default!(Rational, 'f', prec = 2)
fmt_default!(Format.ComplexRational, 'f', prec = 2)

@test fmt(3//4, 10, 2) == "      0.75"
@test fmt(3//4, 10, 1) == "       0.8"
@test fmt(1//2 + 6//2 * im, 23) == "  0.500000 + 3.000000im"
@test fmt(1//2 + 6//2 * im, 15, 2) == "  0.50 + 3.00im"

fmt_default!(Int, :commas, width = 12)
@test fmt(i) == "   1,234,567"
@test fmt(x) == "1234.567890"  # default hasn't changed

fmt_default!(:commas)
@test fmt(i) == "   1,234,567"
@test fmt(x) == "1,234.567890"  # width hasn't changed, but added commas

fmt_default!(Int) # resets Integer defaults
@test fmt(i) == "1234567"
@test fmt(i,:commas) == "1,234,567"

reset!(Int)
fmt_default!(UInt16, 'd', :commas)
@test fmt(0xffff) == "65,535"
fmt_default!(UInt32, UInt16, width=20)
@test fmt(0xfffff) == "           1,048,575"

v = pi
@test fmt(v) == "π"
@test fmt(v; width=10) == "         π"

v = MathConstants.eulergamma
@test fmt(v, 10, 2) == "         γ"


text:
const t=10.0

@inline initial_condition_generator(rng,n,model) = (X=[eps()],φ=[1])

const _first_exit = first_exit_x(0.5*eps(),1.0)
function first_exit_or_t10(model,SFM,SFM0)
    t1 = _first_exit(model,SFM,SFM0)
    i2 = SFM.t>10.0
    i2 && SFM = (t=10.0, X=SFM.X, φ=SFM.φ, n=SFM.m)
    return (Ind=(i2||t1.Ind), SFM=t1.SFM)
end

const tau = first_exit_or_t10

const x = range(0,10,length=6_001)

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to draw rulers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

function _draw_vertical_ruler!(buf::IO, line::Int, num_lines::Int)
    # Get the required spacing.
    num_spaces = floor(Int, log10(abs(num_lines))) + 1

    # Convert to line and apply the padding.
    line_str = lpad(line, num_spaces)

    # Write to the buffer.
    write(buf, " ", line_str, " │")

    return nothing
end

function _get_vertical_ruler_spacing(num_lines::Int)
    # Get the required spacing for the numbers.
    num_spaces = floor(Int, log10(abs(num_lines))) + 1

    # Return the required spacing to draw the vertical ruller.
    return num_spaces + 3
end


text:
__precompile__(false)

module WebApp

using Blink

addListener(win::Window, event::String, callback::Function) = nothing #if active(win) Blink.@dot win on($event,$callback); end

closed = false
isclosed(win::Window) = closed
#Blink.@dot win on("closed", ()->)global closed = true)

const WEB_FILE = abspath(joinpath(@__DIR__,"../web/index.html")) #"C:/Users/Mario/AppData/Local/Julia-1.0.0/JuliaOpenGL/JuliaOpenGL/web/index.html"

webWindow = nothing

function init()
  if Blink.AtomShell.isinstalled() return end
  println("> INIT AtomShell...")
  Blink.AtomShell.install()
  println("\n----------\n")
end

function run()
  init()
  create()
end

function create()
  println("> Create Window...")
  global webWindow = Window()
  addListener(webWindow, "closed", ()->global closed = true)
  println("Read File")
  #f = open(WEB_FILE) do file read(file, String) end
  #println("Write into Window")
  #body!(webWindow,f)
  loadfile(webWindow, WEB_FILE)
  #println("Wait...")
  #while !isclosed(webWindow) yield() end
end

function reload()
  if active(webWindow)
    loadfile(webWindow, WEB_FILE)
  else
    close(webWindow)
    create()
  end
end

function cleanUp()
  #close(webWindow)
  #Blink.AtomShell.uninstall()
end

end #WebApp

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# source modelling from volume image data

using Neuroimaging
data_path = joinpath("..", "..", "..", "test", "data", "test-3d.dat")

t = read_VolumeImage(data_path)

NeuroImaging.plot(t)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add AbstractIntegratorAlgorithm or similar

"""
    bat_integrate(
        posterior::BAT.AnyPosterior,
    )::DensitySampleVector

Calculate the integral (evidence) of `posterior`.

Returns a NamedTuple: (integral = x::Measurement.Measurement, ...)

Result properties not listed here are algorithm-specific and are not part
of the stable BAT API.

`posterior` may be a

* [`BAT.AbstractPosteriorDensity`](@ref)

* [`BAT.DistLikeDensity`](@ref)

* [`BAT.DensitySampleVector`](@ref)

* `Distributions.MultivariateDistribution`

Uses the AHMI algorithm by default.
"""
function bat_integrate end
export bat_integrate

function bat_integrate(posterior::DensitySampleVector)
    hmi_data = HMIData(posterior)
    hm_integrate!(hmi_data)
    Z_signal_v = hmi_data.integralestimates["cov. weighted result"].final.estimate
    result = hmi_data.integralestimates["cov. weighted result"].final
    info = hmi_data.integralestimates

    integral = Measurements.measurement(result.estimate, result.uncertainty)
    (integral = integral, info = info)
end

function bat_integrate(posterior::AnyPosterior)
    npar = totalndof(params_shape(posterior))
    nsamples = 10^5 * npar
    samples = bat_sample(posterior, nsamples).samples::DensitySampleVector
    bat_integrate(samples)
end


text:
# Helper Unions
QuantityOrUnits = Union{Unitful.AbstractQuantity,Unitful.Units}
QuantityOrUnitlike = Union{Unitful.AbstractQuantity, Unitful.Unitlike}

"""
    DimBasis(basis_vectors...) -> DimBasis

Create a dimensional basis for a number of `basis_vectors` (quantities, units or dimensions).
A string identifier can optionally be added to each basis vector.
"""
struct DimBasis{T,N}
    basis_vectors::T
    basis_vector_names::N
    basis_dims
    dim_mat

    function DimBasis(basis_vectors::T, basis_vector_names::N=nothing) where
        T<:AbstractVector{<:QuantityOrUnitlike} where
        N<:Union{Nothing,AbstractVector{<:AbstractString}}
        basis_dims = unique_dims(basis_vectors...)
        dim_mat = dim_matrix(basis_dims, basis_vectors...)
        check_basis(dim_mat)
        return new{T,N}(basis_vectors, basis_vector_names, basis_dims, dim_mat)
    end
end

DimBasis(basis_vectors::Vararg{<:QuantityOrUnitlike}) =
DimBasis([basis_vectors...])

function DimBasis(basis_vectors::Vararg{Pair{<:AbstractString,<:QuantityOrUnitlike}})
    new_basis_vectors = [basis_vector.second for basis_vector in basis_vectors]
    new_basis_vector_names = [basis_vector.first for basis_vector in basis_vectors]
    return DimBasis(new_basis_vectors, new_basis_vector_names)
end

# Do not broadcast DimBasis
Base.Broadcast.broadcastable(basis::DimBasis) = Ref(basis)

# Helper Types
QuantityDimBasis = DimBasis{<:AbstractVector{<:Unitful.AbstractQuantity}}

NamedDimBase = DimBasis{T,<:AbstractVector{<:AbstractString}} where T

"""
    unique_dims(all_values...)

Return a vector of unique dimensions for `all_values`, a set of quantities, units or dimensions.
"""
function unique_dims(all_values::Vararg{Unitful.Dimensions})
    basis_dims = Vector{Type{<:Unitful.Dimension}}()
    for dims in all_values
        union!(basis_dims, typeof.(typeof(dims).parameters[1]))
    end
    return basis_dims
end

unique_dims(all_values::Vararg{QuantityOrUnits}) =
    unique_dims(dimension.(all_values)...)

"""
    dim_matrix(basis_dims, all_values...)

Return the dimensional matrix for a set of basis dimensions `basis_dims` and `all_values`, a set of quantities, units or dimensions.
"""
function dim_matrix(basis_dims::Array{Type{<:Unitful.Dimension}}, all_values::Vararg{Unitful.Dimensions})
    dim_mat = zeros(Rational, length(basis_dims), length(all_values))
    for (dimension_ind,dims) in enumerate(all_values)
        for dimension in typeof(dims).parameters[1]
            basis_dim_ind = findfirst(x -> isa(dimension, x), basis_dims)
            if basis_dim_ind == nothing
                error("$(typeof(dimension)) is no basis dimension!")
            end
            dim_mat[basis_dim_ind,dimension_ind] = dimension.power
        end
    end
    return dim_mat
end

dim_matrix(basis_dims::Array{Type{<:Unitful.Dimension}}, all_values::Vararg{<:QuantityOrUnits}) =
dim_matrix(basis_dims, dimension.(all_values)...)

"""
    check_basis(dim_mat)

Use a dimensional matrix to check if a collection of dimensional vectors is a valid basis.
Throw errors if there are to few basis vectors or if the matrix does not have full rank.
"""
function check_basis(dim_mat)
    if size(dim_mat, 2) < size(dim_mat, 1)
        plr_sgl = (size(dim_mat, 2) == 1) ? "vector" : "vectors"
        error("Invalid basis! There are $(size(dim_mat, 1)) dimensions but only $(size(dim_mat, 2)) basis $(plr_sgl).")
    end

    mat_rank = LinearAlgebra.rank(dim_mat)
    if mat_rank < size(dim_mat, 2)
        if mat_rank == 1
            error("Invalid basis! There are $(size(dim_mat, 2)) basis vectors that are all linear dependent.")
        else
            error("Invalid basis! There are $(size(dim_mat, 2)) basis vectors of which only $(mat_rank) are linear independent.")
        end
    end
    nothing
end


text:
function initialiseMRSresultsDatabase!( mrs::MRS )
    # Clear database.
    if haskey( SQLite.tables( mrs.resultsDB ), :name )
        tablesToDelete = SQLite.tables( mrs.resultsDB ).name
        filter!( tname -> tname ∉ [mrs.mpSim.persDBname, mrs.mpSim.histDBname,
            mrs.mpSim.transDBname], tablesToDelete )
        DBInterface.execute( mrs.resultsDB, "BEGIN TRANSACTION" )
        sqliteCmd = string.( "DROP TABLE `", tablesToDelete, "`" )
        DBInterface.execute.( Ref(mrs.resultsDB), sqliteCmd )
        DBInterface.execute( mrs.resultsDB, "COMMIT" )
    end  # if !isempty( SQLite.tables( mrs.resultsDB ) )

    saveSimulationConfiguration( mrs.mpSim, mrs.resultsDB.file )

    sqliteCmd = string( "CREATE TABLE simokay(",
        "\n    simID TEXT )" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )
end  # initialiseMRSresultsDatabase!( mrs )

function copySnapshots( mrs::MRS, tmpFolder::String )
    tStart = now()

    for ii in 1:mrs.nRuns
        tmpname = joinpath( tmpFolder, string( "tmp", ii, ".sqlite" ) )
        cp( mrs.resultsDB.file, tmpname, force=true )
        newDB = SQLite.DB(tmpname)
        DBInterface.execute( newDB, "DROP TABLE config" )
        DBInterface.execute( newDB, "DROP TABLE simokay" )

        tablenames = SQLite.tables(newDB).name
        DBInterface.execute.( Ref(newDB), string.( "ALTER TABLE `", tablenames,
            "` RENAME TO ", tablenames, ii ) )
    end  # for ii in 1:mrs.nRuns

    tElapsed = ( now() - tStart ).value / 1000.0
    mrs.showInfo && @info string( "Initialising individual databases with start populations took ", tElapsed, " seconds." )
end  # copySnapshots( mrs, tmpFolder )

function copySnapshots( mrs::MRS, tmpFolder::String, nThreads::Int )
    tStart = now()

    for ii in 1:nThreads
        tmpname = joinpath( tmpFolder, string( "tmp", ii, ".sqlite" ) )
        cp( mrs.resultsDB.file, tmpname, force=true )
        newDB = SQLite.DB(tmpname)
        DBInterface.execute( newDB, "DROP TABLE config" )
        DBInterface.execute( newDB, "DROP TABLE simokay" )
    end  # for ii in 1:nThreads

    tElapsed = ( now() - tStart ).value / 1000.0
    mrs.showInfo && @info string( "Initialising temporary databases with start populations took ", tElapsed, " seconds." )
end  # copySnapshots( mrs, tmpFolder, nThreads )

function runsim( mrs::MRS, currentSim::Channel{Int},
    resultsFree::Channel{Bool}, tmpFolder::String )
    n = 0

    while ( ( n = take!( currentSim ) ) <= mrs.nRuns )
        tStart = now()
        put!( currentSim, n + 1 )

        tmpSim = deepcopy( mrs.mpSim )
        GC.gc()
        setSimulationName!( tmpSim, string( tmpSim.simName, n ) )
        setSimulationDatabase!( tmpSim,
            joinpath( tmpFolder, string( "tmp", n ) ) )
        run( tmpSim, saveConfig=false )

        # Copy the results to the main database right away.
        copyRunResult( mrs, tmpSim, resultsFree )
        mrs.nComplete += 1
        tElapsed = ( now() - tStart ).value / 1000.0
        mrs.showInfo && @info string( "Simulation ", mrs.nComplete, " of ",
            mrs.nRuns, " completed after ", tElapsed, " seconds." )
    end  # while ( ( n = take!( currentSim ) ) <= mrs.nRuns )

    put!( currentSim, n )
end  # runsim( mrs, currentSim, resultsFree, tmpFolder )

function runsim( mrs::MRS, threadnum::Int, hasSnapshot::Bool,
    currentSim::Channel{Int}, tmpFolder::String, nThreads::Int )
    tmpSim = deepcopy(mrs.mpSim)
    setSimulationDatabase!( tmpSim,
        joinpath( tmpFolder, string( "tmp", threadnum ) ) )
    DBInterface.execute( tmpSim.simDB,
        "CREATE TABLE inventory( repnr INTEGER )" )

    baseName = tmpSim.simName
    baseNames = (tmpSim.persDBname, tmpSim.transDBname, tmpSim.histDBname)
    n = 0

    while ( ( n = take!(currentSim) ) <= mrs.nRuns )
        tStart = now()
        put!( currentSim, n + 1 )

        setSimulationName!( tmpSim, string( baseName, n ) )
        hasSnapshot && copySnapshot( tmpSim, baseNames )
        run( tmpSim, saveConfig=false, nCommits=nThreads )
        sqliteCmd = string( "INSERT INTO inventory( repnr ) VALUES (", n, ")" )
        DBInterface.execute( tmpSim.simDB, sqliteCmd )

        # Copy the results to the main database right away.
        mrs.nComplete += 1
        tElapsed = ( now() - tStart ).value / 1000.0
        mrs.showInfo && @info string( "Simulation ", mrs.nComplete, " of ",
            mrs.nRuns, " completed after ", tElapsed, " seconds." )
        # GC.gc()
    end  # while ( ( n = take!(currentSim) ) <= mrs.nRuns )

    put!( currentSim, n )
end  # runsim( mrs, threadnum, hasSnapshot, currentSim, tmpFolder, nThreads )

function copySnapshot( mpSim::MPsim, baseNames::NTuple{3,String} )
    newNames = (mpSim.persDBname, mpSim.transDBname, mpSim.histDBname)

    # sqliteCmds = string.( "CREATE TABLE `", newNames, "` LIKE `",
    #     baseNames, "`" )
    DBInterface.execute.( Ref(mpSim.simDB), generateTable.( newNames,
        Ref(newNames[1]), Ref(mpSim), (:pers, :trans, :hist) ) )

    sqliteCmds = string.( "INSERT INTO `", newNames, "` SELECT * FROM `",
        baseNames, "`" )
    DBInterface.execute.( Ref(mpSim.simDB), sqliteCmds )
end  # copySnapshot( sim, baseNames )

function copyResults( mrs::MRS, tmpFolder::String )
    tmpDBs = joinpath.( tmpFolder, readdir(tmpFolder) )

    for tmpDBname in tmpDBs
        tmpDB = SQLite.DB(tmpDBname)

        # Attach temporary database.
        sqliteCmd = string( "ATTACH '", tmpDBname, "' AS origin" )
        DBInterface.execute( mrs.resultsDB, sqliteCmd )
        DBInterface.execute( mrs.resultsDB, "BEGIN TRANSACTION" )

        simNrs = DataFrame(DBInterface.execute( tmpDB,
            "SELECT repnr FROM inventory"))[!, :repnr]
        tableNames = (mrs.mpSim.persDBname, mrs.mpSim.transDBname,
            mrs.mpSim.histDBname)
        copyResults.( Ref(mrs), simNrs, Ref(tableNames) )

        # Detach temporary database.
        DBInterface.execute( mrs.resultsDB, "COMMIT" )
        DBInterface.execute( mrs.resultsDB, "DETACH origin" )

    end  # for tmpDB in tmpDBs
end  # copyResults( mrs, tmpFolder )

function copyResults( mrs::MRS, simNr::Int, tableNames::NTuple{3,String} )
    mpSim = mrs.mpSim
    tmpNames = string.( tableNames, simNr )

    # Remove tables if needed.
    # sqliteCmds = string.( "DROP TABLE IF EXISTS `", tmpNames, "`" )
    # DBInterface.execute.( Ref(mrs.resultsDB), sqliteCmds )

    # Create tables.
    DBInterface.execute.( Ref(mrs.resultsDB), generateTable.( tmpNames,
        Ref(tmpNames[1]), Ref(mpSim), (:pers, :trans, :hist) ) )
#=
    sqliteCmd = string( "CREATE TABLE `", tmpNames[1], "`(",
        "\n    `", mpSim.idKey, "` TEXT NOT NULL PRIMARY KEY,",
        "\n    timeEntered REAL,",
        "\n    timeExited REAL,",
        "\n    ageAtRecruitment REAL,",
        "\n    expectedAttritionTime REAL,",
        "\n    currentNode TEXT,",
        join( string.( "\n    `", collect( keys( mpSim.attributeList ) )
            , "` TEXT," ) ),
        "\n    status TEXT )" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )

    sqliteCmd = string( "CREATE TABLE `", tmpNames[2], "`(",
        "\n    `", mpSim.idKey, "` TEXT,",
        "\n    timeIndex REAL,",
        "\n    transition TEXT,",
        "\n    sourceNode TEXT,",
        "\n    targetNode TEXT,",
        "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
        mpSim.persDBname, "`(`", mpSim.idKey, "`) )" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )

    sqliteCmd = string( "CREATE TABLE `", tmpNames[3], "`(",
        "\n    `", mpSim.idKey, "` TEXT,",
        "\n    timeIndex REAL,",
        "\n    attribute TEXT,",
        "\n    value TEXT,",
        "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
        mpSim.persDBname, "`(`", mpSim.idKey, "`) )" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )
=#
    # Copy results to main database.
    sqliteCmds = string.( "INSERT INTO `", tmpNames, "`",
        "\n    SELECT * FROM origin.`", tmpNames, "`" )
    DBInterface.execute.( Ref(mrs.resultsDB), sqliteCmds )
end  # copyResults( mrs, simNr, tableNames )

function copyRunResult( mrs::MRS, mpSim::MPsim, resultsFree::Channel{Bool} )
    take!( resultsFree )

    # Link simulation result database.
    sqliteCmd = string( "ATTACH '", mpSim.simDB.file, "' AS origin" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )

    # Create the tables.
    tableNames = (mpSim.persDBname, mpSim.transDBname, mpSim.histDBname)
    DBInterface.execute.( Ref(mrs.resultsDB), generateTable.( tableNames,
        Ref(tableNames[1]), Ref(mpSim), (:pers, :trans, :hist) ) )

    # Copy results to main database.
    sqliteCmds = string.( "INSERT INTO `", tableNames, "`",
        "\n    SELECT * FROM origin.`", tableNames, "`" )
    DBInterface.execute.( Ref(mrs.resultsDB), sqliteCmds )

    # Unlink simulation result database.
    DBInterface.execute( mrs.resultsDB, "DETACH origin" )

    # Checking off sim completion.
    sqliteCmd = string( "INSERT INTO simokay VALUES ('", mpSim.simName, "')" )
    DBInterface.execute( mrs.resultsDB, sqliteCmd )
    put!( resultsFree, true )
end  # copyRunResult( mrs, mpSim, resultsFree )

function generateTable( tableName::String, persName::String, mpSim::MPsim,
    type::Symbol )
    if type === :pers
        sqliteCmd = string( "CREATE TABLE `", tableName, "`(",
            "\n    `", mpSim.idKey, "` TEXT NOT NULL PRIMARY KEY,",
            "\n    timeEntered REAL,",
            "\n    timeExited REAL,",
            "\n    ageAtRecruitment REAL,",
            "\n    expectedAttritionTime REAL,",
            "\n    currentNode TEXT,",
            "\n    inNodeSince REAL,",
            join( string.( "\n    `", collect( keys( mpSim.attributeList ) )
                , "` TEXT," ) ),
            "\n    status TEXT )" )
    elseif type === :trans
        sqliteCmd = string( "CREATE TABLE `", tableName, "`(",
            "\n    `", mpSim.idKey, "` TEXT,",
            "\n    timeIndex REAL,",
            "\n    transition TEXT,",
            "\n    sourceNode TEXT,",
            "\n    targetNode TEXT,",
            "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
            persName, "`(`", mpSim.idKey, "`) )" )
    elseif type === :hist
        sqliteCmd = string( "CREATE TABLE `", tableName, "`(",
            "\n    `", mpSim.idKey, "` TEXT,",
            "\n    timeIndex REAL,",
            "\n    attribute TEXT,",
            "\n    value TEXT,",
            "\n    FOREIGN KEY (`", mpSim.idKey, "`) REFERENCES `",
            persName, "`(`", mpSim.idKey, "`) )" )
    end

    sqliteCmd
end  # generateTable( tableName, persName, type, mpSim )


text:
using Jutils.Functions
using Jutils.Integration
using Jutils.Mesh
using Jutils.Topologies

function bench_overhead!(suite)
    func = Constant(1.0)
    suite["cmp"]["constant"] = @benchmarkable compile($func)
end

function bench_1d_lag_mass!(suite)
    key = ("1d", "lag", "mass")
    domain, geom = line(1000)

    suite["lpe"]["1d"] = @benchmarkable collect($domain)

    mass = outer(mkbasis(domain, Lagrange, 1))
    suite["cmp"][key] = @benchmarkable compile($mass; dense=false)

    mass = compile(mass; dense=false)
    suite["mke"][key] = @benchmarkable callable($mass)
    suite["itg"][key] = @benchmarkable integrate($mass, $domain, 5)
end

function bench_1d_lag_lapl!(suite)
    key = ("1d", "lag", "lapl")
    domain, geom = line(1000)

    mass = outer(grad(mkbasis(domain, Lagrange, 1), geom)[:,1])
    suite["cmp"][key] = @benchmarkable compile($mass; dense=false)

    mass = compile(mass; dense=false)
    suite["mke"][key] = @benchmarkable callable($mass)
    suite["itg"][key] = @benchmarkable integrate($mass, $domain, 5)
end

function bench_2d_lag_mass!(suite)
    key = ("2d", "lag", "mass")
    domain, geom = rectilinear(30, 30)

    suite["lpe"]["2d"] = @benchmarkable collect($domain)

    mass = outer(mkbasis(domain, Lagrange, 1))
    suite["cmp"][key] = @benchmarkable compile($mass; dense=false)

    mass = compile(mass; dense=false)
    suite["mke"][key] = @benchmarkable callable($mass)
    suite["itg"][key] = @benchmarkable integrate($mass, $domain, 5)
end

function bench_2d_lag_lapl!(suite)
    key = ("2d", "lag", "lapl")
    domain, geom = rectilinear(30, 30)

    mass = outer(grad(mkbasis(domain, Lagrange, 1), geom))
    mass = dropdims(sum(mass, (3,)), (3,))
    suite["cmp"][key] = @benchmarkable compile($mass; dense=false)

    mass = compile(mass; dense=false)
    suite["mke"][key] = @benchmarkable callable($mass)
    suite["itg"][key] = @benchmarkable integrate($mass, $domain, 5)
end

function runbench()
    suite = BenchmarkGroup()
    suite["cmp"] = BenchmarkGroup()
    suite["itg"] = BenchmarkGroup()
    suite["lpe"] = BenchmarkGroup()
    suite["mke"] = BenchmarkGroup()

    bench_overhead!(suite)
    bench_1d_lag_mass!(suite)
    bench_1d_lag_lapl!(suite)
    bench_2d_lag_mass!(suite)
    bench_2d_lag_lapl!(suite)

    tune!(suite)
    return run(suite)
end

max_stars_repo_name: TheBB/Jutils

max_issues_repo_name: TheBB/Jutils

max_forks_repo_name: TheBB/Jutils


text:
using DDCModelsExample

# Define test_speed function
function test_speed(β,δ,ρ,nperiod,nfirm,supportX,nrepeat)
    nsupportX = length(supportX)

    # Construct Π matrix
    Π = 1 ./ (1 .+ abs.(repeat(supportX',nsupportX,1) .- repeat(supportX,1,nsupportX)))
    Π ./= sum(Π;dims=2)

    m = DDCbasic(;supportX,Π,β,δ,ρ)

    t = 0

    for r in 1:nrepeat
        m = DDCbasic(;supportX,Π,β,δ,ρ)
        m.options.seed = 098471 + r
        X,choices = gen_data(m,nperiod,nfirm)
        t += @elapsed begin
            estimate_model(m,X,choices; startvals = [-1;-0.1;0.5])
            calc_stderr(m,X,choices)
        end
    end
    return t
end

# Define arguments
nperiod = 100
nfirm = 10000
supportX = collect(1:5)
β = [-0.1*length(supportX);0.2]
δ = [0.0;1.0]
ρ = 0.95
nrepeat = 10

# Run first time: includes precompilation time
tpc = test_speed(β,δ,ρ,nperiod,nfirm,supportX,nrepeat)

# Run second time: no precompilation time
tnopc = test_speed(β,δ,ρ,nperiod,nfirm,supportX,nrepeat)

max_stars_repo_name: rgreminger/DDCModelsExample

max_issues_repo_name: rgreminger/DDCModelsExample

max_forks_repo_name: rgreminger/DDCModelsExample


text:
#
# Align Minimum Support Sampler
#
# adds a fixed value to sampled value so that minimum of requested support is aligned to minimum of support of
# the subsampler
# the subsampler must provide a minimumsupport method
# note: this is intended to be used with distributions with finite minimum support, but no error is explicitly raised
# if it is infinite
#

type AlignMinimumSupportSampler <: ModifyingSampler
	subsampler::Sampler
	function AlignMinimumSupportSampler(subsampler::Sampler)
		new(subsampler)
	end
end

function sample(s::AlignMinimumSupportSampler, support, cc::ChoiceContext)
	delta = support[1] - minimumsupport(s.subsampler)
	x, trace = sample(s.subsampler, (support[1]-delta, support[2]-delta), cc)
	x + delta, Dict{Symbol, Any}(:sub=>trace, :delta=>delta)
end

amendtrace(s::AlignMinimumSupportSampler, trace, x) = amendtrace(s.subsampler, trace[:sub], x - trace[:delta])


text:
import LinearAlgebra
import LinearAlgebra: isdiag, ishermitian, isposdef, isposdef!, issuccess, issymmetric, istril, istriu,
     tr, det, dot, cross, adjoint, adjoint!, transpose, transpose!, diag, diagm, diagind,
     svdvals, svdvals!, svd, svd!, eigvals, eigvals!, eigvecs, eigen, eigen!

# Gaius.blocked_mul!(c::Matrix{Robust32}, a::Matrix{Robust32}, b::Matrix{Robust32}) = rewrap(blocked_mul!(rewrap(c), rewrap(a), rewrap(b)))
# LinearAlgebra.mul!(c::Matrix{Robust32}, a::Matrix{Robust32}, b::Matrix{Robust32}) = blocked_mul!(c, a, b)

for F in (:+, :-, :*, :/, :\)
  @eval begin
    $F(x::Vector{Robust32}, y::Vector{Robust32}) = rewrap($F(rewrap(x), rewrap(y)))
    $F(x::Vector{Robust32}, y::Vector{Float64})  = rewrap($F(rewrap(x), y))
    $F(x::Vector{Float64}, y::Vector{Robust32})  = rewrap($F(x, rewrap(y)))
    $F(x::Vector{Robust32}, y::Vector{Float32})  = rewrap($F(rewrap(x), y))
    $F(x::Vector{Float32}, y::Vector{Robust32})  = rewrap($F(x, rewrap(y)))

    $F(x::Matrix{Robust32}, y::Matrix{Robust32}) = rewrap($F(rewrap(x), rewrap(y)))
    $F(x::Matrix{Robust32}, y::Matrix{Float64})  = rewrap($F(rewrap(x), y))
    $F(x::Matrix{Float64}, y::Matrix{Robust32})  = rewrap($F(x, rewrap(y)))
    $F(x::Matrix{Robust32}, y::Matrix{Float32})  = rewrap($F(rewrap(x), y))
    $F(x::Matrix{Float32}, y::Matrix{Robust32})  = rewrap($F(x, rewrap(y)))
  end
end

for F in (:tr, :det)
    @eval $F(x::Matrix{Robust32}) = Rob32($F(rewrap(x)))
end

adjoint(x::Matrix{Robust32}) = Adjoint(x)
transpose(x::Matrix{Robust32}) = Transpose(x)

for F in (:isdiag, :ishermitian, :isposdef, :isposdef!, :issuccess, :issymmetric, :istril, :istriu)
  @eval $F(x::Matrix{Robust32}) = $F(rewrap(x))
end

dot(x::Array{Robust32,N}, y::Array{Robust32,N}) where {N} = Rob32(dot(rewrap(x), rewrap(y)))
cross(x::Array{Robust32,1}, y::Array{Robust32,1}) where {N} = rewrap(cross(rewrap(x), rewrap(y)))

for F in (:inv, :sqrt, :exp, :log,
          :sin, :cos, :tan, :csc, :sec, :cot, :asin, :acos, :atan, :acsc, :asec, :acot,
          :sinh, :cosh, :tanh, :csch, :sech, :coth, :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
    @eval $F(x::Matrix{Robust32}) = reinterpret(Robust32, $F(rewrap(x)))
end

# diag diagm diagind
diag(x::Matrix{Robust32}) = rewrap(diag(rewrap(Float64, x)))

svdvals(A::Matrix{Robust32}; kw...) = rewrap(svdvals(rewrap(A); kw...))
eigvals(A::Matrix{Robust32}; kw...) = rewrap(eigvals(rewrap(A); kw...))
svdvals!(A::Matrix{Robust32}) = rewrap(svdvals!(rewrap(A)))
eigvals!(A::Matrix{Robust32}) = rewrap(eigvals!(rewrap(A)))

function svd(x::Matrix{Robust32})
    u, s, v = svd(rewrap(x))
    U = rewrap(u)
    S = rewrap(s)
    V = adjoint(rewrap(adjoint(v)))
    return SVD{Robust32,Robust32,Matrix{Robust32}}(U, S, V)
end

function vd!(x::Matrix{Robust32}; kw...)
    u, s, v = svd!(rewrap(x); kw...)
    U = rewrap(u)
    S = rewrap(s)
    V = adjoint(rewrap(adjoint(v)))
    return SVD{Robust32,Robust32,Matrix{Robust32}}(U, S, V)
end

function eigen(x::Matrix{Robust32}; kw...)
    v, m = eigen(rewrap(x); kw...)
    V = rewrap(v)
    M = rewrap(m)
    return Eigen{ComplexR32,ComplexR32,Matrix{ComplexR32},Vector{ComplexR32}}(V,M)
end

function eigen!(x::Matrix{Robust32}; kw...)
    v, m = eigen!(rewrap(x); kw...)
    V = rewrap(v)
    M = rewrap(m)
    return Eigen{ComplexR32,ComplexR32,Matrix{ComplexR32},Vector{ComplexR32}}(V,M)
end

function lu(x::Matrix{Robust32}, pivot=Val{true}; check=true)
   xx = rewrap(x)
   res = lu(xx, pivot(); check=check)
   yy = rewrap(res.factors)
   return LU{Robust32, Matrix{Robust32}}(yy, res.ipiv, res.info)
end

function lu!(x::Matrix{Robust32}, pivot=Val{true}; check=true)
   xx = rewrap(x)
   res = lu(xx, pivot(); check=check)
   yy = rewrap(res.factors)
   xx[:] = yy[:]
   return LU{Robust32, Matrix{Robust32}}(yy, res.ipiv, res.info)
end

function cholesky(x::Matrix{Robust32}, pivot=Val{false}; check=true)
   xx = rewrap(x)
   res = cholesky(xx, pivot(); check=check)
   yy = rewrap(res.factors)
   return Cholesky{Robust32, Matrix{Robust32}}(yy, res.uplo, res.info)
end

function cholesky!(x::Matrix{Robust32}, pivot=Val{false}; check=true)
   xx = rewrap(x)
   res = cholesky(xx, pivot(); check=check)
   yy = rewrap(res.factors)
   x[:] = yy[:]
   return Cholesky{Robust32, Matrix{Robust32}}(yy, res.uplo, res.info)
end

function lmul!(x::Matrix{Robust32}, y::Matrix{Robust32})
   xx = rewrap(x)
   yy = rewrap(y)
   res = lmul!(xx, yy)
   return rewrap(res)
end

function lmul!(x::AbstractMatrix{Robust32}, y::AbstractVector{Robust32})
   xx = rewrap(x)
   yy = rewrap(y)
   res = lmul!(xx, yy)
   return rewrap(res)
end

#=
[  "/", "\\",
   # "adjoint", "adjoint!",
   "axpby!", "axpy!", "bunchkaufman", "bunchkaufman!", "cholesky", "cholesky!", "cond", "condskeel",
   "copy_transpose!", "copyto!",
   # "cross",
   # "det", "diag",
   "diagind", "diagm",
   # "dot",
   # "eigen", "eigen!",
   "eigmax", "eigmin",
   # "eigvals", "eigvals!",
    "eigvecs",
   "factorize", "givens", "hessenberg", "hessenberg!",
   # "isdiag", "ishermitian", "isposdef", "isposdef!", "issuccess", "issymmetric", "istril", "istriu",
   "kron", "ldiv!", "ldlt", "ldlt!", "lmul!", "logabsdet", "logdet",
   "lowrankdowndate", "lowrankdowndate!", "lowrankupdate", "lowrankupdate!",
   "lq", "lq!",
   # "lu", "lu!",
   "lyap",
   "mul!", "norm", "normalize", "normalize!", "nullspace", "opnorm",
   "ordschur", "ordschur!", "pinv", "qr", "qr!",
   "rank", "rdiv!", "reflect!", "rmul!", "rotate!", "schur", "schur!",
   # "svd", "svd!", "svdvals", "svdvals!",
   "sylvester",
   # "tr", "transpose", "transpose!",
   "tril", "tril!", "triu", "triu!", "×", "⋅" ]

   (:UnitUpperTriangular, :UnitLowerTriangular, :UpperHessenberg, :LowerHessenbert, :UpperTriangular, :LowerTriangular, :Tridiagonal, :SymTridiagonal, :Symmetric)

   ["Adjoint", "BLAS", "Bidiagonal", "BunchKaufman", "Cholesky", "CholeskyPivoted",
    "Diagonal", "Eigen", "Factorization", "GeneralizedEigen", "GeneralizedSVD", "GeneralizedSchur",
    "Hermitian", "Hessenberg", "I", "LAPACK", "LAPACKException", "LDLt", "LQ", "LU",
    "LowerTriangular", "QR", "QRPivoted", "SVD",
    "Schur", "SymTridiagonal", "Symmetric", "Transpose", "Tridiagonal",
    "UniformScaling", "UnitLowerTriangular", "UnitUpperTriangular", "UpperHessenberg", "UpperTriangular"]

julia> listnames(:eigen)
  list(23)
  –––––––– ––––––– –––––––– –––––––– –––––––– –––––– ––––––– ––––––––– ––––––– –––––––– ––––––– –––––––––
  chebspec circul  dingdong forsythe grcar    invol  minij   oscillate pascal  rosser   tridiag wilkinson
  chow     clement fiedler  frank    hadamard lotkin neumann parter    poisson sampling wathen

julia> listnames(:graph)
  list(3)
  ––––––– ––––––– ––––––––––
  erdrey  gilbert smallworld

julia> listnames(:illcond)
  list(20)
  –––––––– ––––– ––––––– ––––– –––––– ––––––––– –––––– ––––––– ––––––– ––––
  cauchy   frank hilb    invol kms    moler     pascal prolate rosser  triw
  forsythe golub invhilb kahan lotkin oscillate pei    randsvd tridiag vand

julia> listnames(:inverse)
  list(21)
  –––––––– –––––––– –––––––– ––––––– ––––– –––––– ––––– –––––– ––––––– ––––––– ––––
  cauchy   fiedler  hadamard invhilb kahan lehmer magic moler  pei     tridiag vand
  clement  forsythe hilb     invol   kms   lotkin minij pascal poisson triw

julia> listnames(:posdef)
  list(14)
  –––––––– –––––– –––– ––––––– ––– –––––– ––––– ––––– ––––––––– –––––– ––– ––––––– ––––––– ––––––
  cauchy   circul hilb invhilb kms lehmer minij moler oscillate pascal pei poisson tridiag wathen

julia> listnames(:symmetric)
  list(21)
  –––––––– –––––––– ––––––– ––––––– –––––– ––––– ––––––––– ––––––– –––––––– ––––––– –––––––––
  cauchy   clement  fiedler hilb    kms    minij oscillate pei     prolate  tridiag wilkinson
  circul   dingdong hankel  invhilb lehmer moler pascal    poisson randcorr wathen

=#


text:
#fp_growth_performance.jl

include("fp_growth.jl")

#Testing and benchmarking
adult = readtable("/Users/pstey/projects_code/association_rules_jl/data/adult.csv")
T = make_transactions(adult)

find_frequent_patterns(T[1:2], 1)

Profile.clear_malloc_data()

@time find_frequent_patterns(T[1:500], 50);


text:
module Crypto

macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, eval(lib))) ), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

@ctypedef EVP_MD Void
@ctypedef EVP_MD_CTX Void
@ctypedef ENGINE Void

typealias size_t Csize_t

@unix_only const libcrypto = "libcrypto"
@windows_only const libcrypto = Pkg.dir("WinRPM","deps","usr","$(Sys.ARCH)-w64-mingw32","sys-root","mingw","bin","libcrypto-10")

@c Ptr{UInt8} HMAC (Ptr{EVP_MD}, Ptr{Void}, Int32, Ptr{UInt8}, size_t, Ptr{UInt8}, Ptr{UInt32}) libcrypto
@c Ptr{UInt8} MD5 (Ptr{UInt8}, size_t, Ptr{UInt8}) libcrypto
@c Ptr{UInt8} SHA256 (Ptr{UInt8}, size_t, Ptr{UInt8}) libcrypto

@c Ptr{EVP_MD} EVP_md5 () libcrypto
@c Ptr{EVP_MD} EVP_sha1 () libcrypto
@c Ptr{EVP_MD} EVP_sha256 () libcrypto

@c Ptr{EVP_MD_CTX} EVP_MD_CTX_create () libcrypto
@c Int32 EVP_DigestInit_ex (Ptr{EVP_MD_CTX}, Ptr{EVP_MD}, Ptr{ENGINE}) libcrypto
@c Int32 EVP_DigestUpdate (Ptr{EVP_MD_CTX}, Ptr{Void}, size_t) libcrypto
@c Int32 EVP_DigestFinal_ex (Ptr{EVP_MD_CTX}, Ptr{UInt8}, Ptr{UInt32}) libcrypto
@c Union{} EVP_MD_CTX_destroy (Ptr{EVP_MD_CTX},) libcrypto

hmacsha256_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}) =  hmacsha_digest(s, k, EVP_sha256(), 32)
export hmacsha256_digest

hmacsha1_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}) = hmacsha_digest(s, k, EVP_sha1(), 20)
export hmacsha1_digest

function hmacsha_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}, evp, dgst_len)
    if evp == C_NULL error("EVP_sha1() failed!") end

    sig = zeros(UInt8, dgst_len)
    sig_len = zeros(UInt32, 1)

    if isa(k, AbstractString)
        k = convert(Array{UInt8}, k)
    end

    if HMAC(evp, k, length(k), s, length(s), sig, sig_len) == C_NULL error("HMAC() failed!") end
    if (sig_len[1] != dgst_len) error("Wrong length of signature!") end

    return sig
end

function md5_file(s::AbstractString)
    f = open(s)
    md = nothing
    try
        md = md5(f)
    catch e
        rethrow(e)
    finally
        close(f)
    end
    md
end
export md5_file

function md5(s::AbstractString)
    md = zeros(UInt8, 16)
    assert(MD5(s, length(s), md) != C_NULL)
    return md
end

function md5(s::IO)
    evp_md_ctx = EVP_MD_CTX_create()
    assert(evp_md_ctx != C_NULL)

    md = zeros(UInt8, 16)
    try
        evp_md = EVP_md5()
        assert(evp_md != C_NULL)

        rc = EVP_DigestInit_ex(evp_md_ctx, evp_md, C_NULL)
        assert(rc == 1)

        while (!eof(s))
            b = read(s, UInt8, min(nb_available(s), 65536))    # Read in 64 K chunks....

            rc = EVP_DigestUpdate(evp_md_ctx, b, length(b));
            assert(rc == 1)
        end

        rc = EVP_DigestFinal_ex(evp_md_ctx, md, C_NULL)
        assert(rc == 1)

    finally
        EVP_MD_CTX_destroy(evp_md_ctx)
    end

    return md
end
export md5

function sha256(s::AbstractString)
    sha = zeros(UInt8, 32)
    assert(SHA256(s, length(s), sha) != C_NULL)
    return sha
end

end # Module end


text:
using JuMP, Ipopt
using Statistics
using LinearAlgebra

mutable struct MdlParams
    dsafety::Float64
    dt::Float64
    MdlParams(dsafety=10.0, dt=0.250) = new(dsafety, dt)
end

mutable struct MpcParams
	N::Int64
	Q::Array{Float64,2}
	R::Array{Float64,2}
	MpcParams(N=16, Q=Diagonal([1, 50]), R=Diagonal([0.001])) = new(N, Q, R)
end

mutable struct MpcModel
	s0::Array{JuMP.NonlinearParameter, 1}
	sRef::Array{JuMP.NonlinearParameter, 2}
	uRef::Array{JuMP.NonlinearParameter, 2}

	obstacles::Array{JuMP.NonlinearParameter, 2}

	s::Array{JuMP.VariableRef, 2}
	u::Array{JuMP.VariableRef, 2}
	λ::Array{JuMP.VariableRef, 1}

	costS::JuMP.NonlinearExpression
	costU::JuMP.NonlinearExpression
	mdl::JuMP.Model

    dt::Float64
	N::Int64

	function MpcModel(mpcParams::MpcParams, mdlParams::MdlParams)
		println("Building mpc controller ...")

		m = new()

		# get model parameters
		dt = mdlParams.dt
		dsafety = mdlParams.dsafety

		# get mpc parameters
		N = mpcParams.N
		Q = mpcParams.Q
		R = mpcParams.R

		# Create Model
		mdl = Model(Ipopt.Optimizer)
		set_optimizer_attributes(mdl, "max_iter" => 5000)

		@variable(mdl, 0 <= λ[i=1:10] <=1 , start = 1)

		# Create variables (will be optimized)
		smin::Vector{Float64} = [0.0; 0.0] # s sd
		smax::Vector{Float64} = [1000.0; 25.0] # s sd
		@variable(mdl, smin[i] <= s[i=1:2, t=1:(N+1)] <= smax[i], start = [0.0, 20.0][i])

		umin::Vector{Float64} = [-4.0] # decel max
		umax::Vector{Float64} = [2.0] # accel max
		@variable(mdl, umin[i] <= u[i=1:1, t=1:N] <= umax[i], start = 0)

		@NLparameter(mdl, s0[i=1:2] == [0.0, 20.0][i])
		sref::Vector{Float64} = [200.0; 20.0] # s sd
		@NLparameter(mdl, sRef[j=1:2,i=1:N+1] == sref[j])
		@NLparameter(mdl, uRef[j=1:1,i=1:N] == 0)

		# System dynamics
		@NLconstraint(mdl, [i=1:2], s[i,1] == s0[i]) # initial condition

		for i=1:N # CA model
			@constraint(mdl, s[1,i+1] == s[1,i] + dt*s[2,i] + 0.5*dt^2*u[1,i]) # s
			@constraint(mdl, s[2,i+1] == s[2,i] + dt*u[1,i]) # sd
		end

		# Cost definitions
		# Input cost
		# ---------------------------------
		@NLexpression(mdl, costU, 0.5*sum(R[j,j] * sum((u[j,i] - uRef[j,i])^2 for i=1:N) for j = 1:1))

		# State cost
		# ---------------------------------
		@NLexpression(mdl, costS, 0.5*sum(Q[j,j] * sum((s[j,i] - sRef[j,i])^2 for i=2:N+1) for j = 1:2))

		@NLexpression(mdl, binaryλ, 10*sum((λ[i]-0.5)^2 for i =1:10))

		# Objective function
		@NLobjective(mdl, Min, costS + costU - binaryλ)
		#@NLobjective(mdl, Min, costS - binaryλ)
		#@NLobjective(mdl, Min, costS + costU - 1000*(λ-0.5)^2)
		#@NLobjective(mdl, Min, costS - (λ-0.5)^2)

		#for i=1:10
		#	@constraint(mdl, λ[i]*(1-λ[i]) == 0) # Binary ...
		#end

		# First solve
		println("Attempting first solve ...")
		optimize!(mdl)
		sol_stat = termination_status(mdl)
		println("Finished solve 1: $sol_stat")

		#delete(mdl, con)

		m.mdl = mdl
		m.λ = λ
		m.dt = dt
		m.N = N
		m.s0 = s0
		m.s = s
		m.u = u
		m.costS = costS
		m.costU = costU
		m.sRef = sRef
		m.uRef = uRef
		return m
	end
end

#function SolveMpcProblem(mdl::MpcModel, mpcSol::MpcSol, sCurr::Array{Float64}, sRef::Array{Float64,2}, uRef::Array{Float64,2})
function act(mpc::MpcModel, state::Array{Float64}, obstacles::Array{Tuple{Float64, Float64}, 1})::Float64

	# update current initial condition
	for i in 1:length(state)-1
		set_value(mpc.s0[i], state[i])
	end

	#println("state=$state")
	constr = []
	count, i = 1, 1
	while i<=min(10,length(obstacles)) && count <= 10
		tcross, scross = obstacles[i]
		tcrossd = floor(Int, (tcross-state[3])/mpc.dt)
		#println("tcross=$tcross tcrossd=$tcrossd scross=$scross")
		constrained::Bool = (tcrossd > 0) && (tcrossd < mpc.N)
		if constrained
			#con = @constraint(mpc.mdl, (mpc.s[1,tcrossd+1] - scross)^2 >= 12.0^2  )
			#push!(constr, con)

			# It is a hand-engineered OR (Ideally λ should be a binary variable)
			con = @constraint(mpc.mdl, mpc.λ[count] * (mpc.s[1,tcrossd+1]  - (scross + 12)) >= 0)
			push!(constr, con)
			con1 = @constraint(mpc.mdl, (1.0 - mpc.λ[count]) * (mpc.s[1,tcrossd+1] - (scross - 12))  <= 0)
			push!(constr, con1)
			count += 1
		end
		i += 1
	end

	println("Attempting solve ...")
	optimize!(mpc.mdl)
	sol_stat = termination_status(mpc.mdl)
	println("Finished solve: $sol_stat")

	sol_u = value.(mpc.u)
	#println("sol_u $sol_u")
	sol_s = value.(mpc.s)
	#println("sol_s $sol_s")

	sol_λ = value.(mpc.λ)
	println("sol_λ $sol_λ")

	for i in 1:length(constr)
		delete(mpc.mdl, constr[i])
	end
	#exit()

	return sol_u[1]
end

#################### Not Used #######################

# abstract type OCP end
#
# @with_kw mutable struct Mpc <: OCP
# 	T::Float64 = 20.0
# 	dt::Float64 = 0.250
# 	n::Int64 = T/dt
# 	dsafety::Float64 = 10.0
#
# 	smin::Vector{Float64} = [0.0; 0.0]
# 	smax::Vector{Float64} = [1000.0; 25.0]
# 	umin::Float64 = -4.0
# 	umax::Float64 = 2.0
#
# 	# Transtion model
# 	Ts::Array{Float64,2} = [1.0 dt; 0  1]
# 	Ta::Array{Float64,1} = [dt^2/2; dt]
#
# 	s::Array{JuMP.VariableRef, 2} = nothing
#
# 	mdp = nothing
# 	model = nothing
# end


text:
# this needs to be defined outside of the module to make
# Julia print type name without module name when displaying it
struct ⛵⛵⛵⛵⛵
end
Base.show(io::IO, ::⛵⛵⛵⛵⛵) = show(io, "⛵")

module TestShow

using DataFrames, Dates, Random, Test

import Main: ⛵⛵⛵⛵⛵

function capture_stdout(f::Function)
    oldstdout = stdout
    rd, wr = redirect_stdout()
    f()
    redirect_stdout(oldstdout)
    size = displaysize(rd)
    close(wr)
    str = read(rd, String)
    close(rd)
    str, size
end

@testset "Basic show test with allrows and allcols" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0], D = ['\'', '∀', '\$', '\n'])

    refstr = """
    4×4 DataFrame
    │ Row │ A     │ B           │ C       │ D    │
    │     │ Int64 │ String      │ Float32 │ Char │
    ├─────┼───────┼─────────────┼─────────┼──────┤
    │ 1   │ 1     │ x"          │ 1.0     │ '\\'' │
    │ 2   │ 2     │ ∀ε>0: x+ε>x │ 2.0     │ '∀'  │
    │ 3   │ 3     │ z\$          │ 3.0     │ '\$'  │
    │ 4   │ 4     │ A\\nC        │ 4.0     │ '\\n' │"""

    for allrows in [true, false], allcols in [true, false]
        io = IOBuffer()
        show(io, df, allcols=allcols, allrows=allrows)
        str = String(take!(io))
        @test str == refstr
        io = IOBuffer()
        show(io, MIME("text/plain"), df, allcols=allcols, allrows=allrows)
        str = String(take!(io))
        @test str == refstr
    end

    df = DataFrame(A = Vector{String}(undef, 3))
    @test sprint(show, df) == """
    3×1 DataFrame
    │ Row │ A      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ #undef │
    │ 2   │ #undef │
    │ 3   │ #undef │"""
end

@testset "displaysize test" begin
    df_big = DataFrame(reshape(Int64(10000001):Int64(10000000+25*5), 25, 5))

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame. Omitted printing of 2 columns
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 10000001 │ 10000026 │ 10000051 │
    ⋮
    │ 24  │ 10000024 │ 10000049 │ 10000074 │
    │ 25  │ 10000025 │ 10000050 │ 10000075 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allcols=true)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 10000001 │ 10000026 │ 10000051 │
    ⋮
    │ 24  │ 10000024 │ 10000049 │ 10000074 │
    │ 25  │ 10000025 │ 10000050 │ 10000075 │

    │ Row │ x4       │ x5       │
    │     │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┤
    │ 1   │ 10000076 │ 10000101 │
    ⋮
    │ 24  │ 10000099 │ 10000124 │
    │ 25  │ 10000100 │ 10000125 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allrows=true, allcols=true)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 10000001 │ 10000026 │ 10000051 │
    │ 2   │ 10000002 │ 10000027 │ 10000052 │
    │ 3   │ 10000003 │ 10000028 │ 10000053 │
    │ 4   │ 10000004 │ 10000029 │ 10000054 │
    │ 5   │ 10000005 │ 10000030 │ 10000055 │
    │ 6   │ 10000006 │ 10000031 │ 10000056 │
    │ 7   │ 10000007 │ 10000032 │ 10000057 │
    │ 8   │ 10000008 │ 10000033 │ 10000058 │
    │ 9   │ 10000009 │ 10000034 │ 10000059 │
    │ 10  │ 10000010 │ 10000035 │ 10000060 │
    │ 11  │ 10000011 │ 10000036 │ 10000061 │
    │ 12  │ 10000012 │ 10000037 │ 10000062 │
    │ 13  │ 10000013 │ 10000038 │ 10000063 │
    │ 14  │ 10000014 │ 10000039 │ 10000064 │
    │ 15  │ 10000015 │ 10000040 │ 10000065 │
    │ 16  │ 10000016 │ 10000041 │ 10000066 │
    │ 17  │ 10000017 │ 10000042 │ 10000067 │
    │ 18  │ 10000018 │ 10000043 │ 10000068 │
    │ 19  │ 10000019 │ 10000044 │ 10000069 │
    │ 20  │ 10000020 │ 10000045 │ 10000070 │
    │ 21  │ 10000021 │ 10000046 │ 10000071 │
    │ 22  │ 10000022 │ 10000047 │ 10000072 │
    │ 23  │ 10000023 │ 10000048 │ 10000073 │
    │ 24  │ 10000024 │ 10000049 │ 10000074 │
    │ 25  │ 10000025 │ 10000050 │ 10000075 │

    │ Row │ x4       │ x5       │
    │     │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┤
    │ 1   │ 10000076 │ 10000101 │
    │ 2   │ 10000077 │ 10000102 │
    │ 3   │ 10000078 │ 10000103 │
    │ 4   │ 10000079 │ 10000104 │
    │ 5   │ 10000080 │ 10000105 │
    │ 6   │ 10000081 │ 10000106 │
    │ 7   │ 10000082 │ 10000107 │
    │ 8   │ 10000083 │ 10000108 │
    │ 9   │ 10000084 │ 10000109 │
    │ 10  │ 10000085 │ 10000110 │
    │ 11  │ 10000086 │ 10000111 │
    │ 12  │ 10000087 │ 10000112 │
    │ 13  │ 10000088 │ 10000113 │
    │ 14  │ 10000089 │ 10000114 │
    │ 15  │ 10000090 │ 10000115 │
    │ 16  │ 10000091 │ 10000116 │
    │ 17  │ 10000092 │ 10000117 │
    │ 18  │ 10000093 │ 10000118 │
    │ 19  │ 10000094 │ 10000119 │
    │ 20  │ 10000095 │ 10000120 │
    │ 21  │ 10000096 │ 10000121 │
    │ 22  │ 10000097 │ 10000122 │
    │ 23  │ 10000098 │ 10000123 │
    │ 24  │ 10000099 │ 10000124 │
    │ 25  │ 10000100 │ 10000125 │"""

    io = IOContext(IOBuffer(), :displaysize=>(11,40), :limit=>true)
    show(io, df_big, allrows=true, allcols=false)
    str = String(take!(io.io))
    @test str == """
    25×5 DataFrame. Omitted printing of 2 columns
    │ Row │ x1       │ x2       │ x3       │
    │     │ Int64    │ Int64    │ Int64    │
    ├─────┼──────────┼──────────┼──────────┤
    │ 1   │ 10000001 │ 10000026 │ 10000051 │
    │ 2   │ 10000002 │ 10000027 │ 10000052 │
    │ 3   │ 10000003 │ 10000028 │ 10000053 │
    │ 4   │ 10000004 │ 10000029 │ 10000054 │
    │ 5   │ 10000005 │ 10000030 │ 10000055 │
    │ 6   │ 10000006 │ 10000031 │ 10000056 │
    │ 7   │ 10000007 │ 10000032 │ 10000057 │
    │ 8   │ 10000008 │ 10000033 │ 10000058 │
    │ 9   │ 10000009 │ 10000034 │ 10000059 │
    │ 10  │ 10000010 │ 10000035 │ 10000060 │
    │ 11  │ 10000011 │ 10000036 │ 10000061 │
    │ 12  │ 10000012 │ 10000037 │ 10000062 │
    │ 13  │ 10000013 │ 10000038 │ 10000063 │
    │ 14  │ 10000014 │ 10000039 │ 10000064 │
    │ 15  │ 10000015 │ 10000040 │ 10000065 │
    │ 16  │ 10000016 │ 10000041 │ 10000066 │
    │ 17  │ 10000017 │ 10000042 │ 10000067 │
    │ 18  │ 10000018 │ 10000043 │ 10000068 │
    │ 19  │ 10000019 │ 10000044 │ 10000069 │
    │ 20  │ 10000020 │ 10000045 │ 10000070 │
    │ 21  │ 10000021 │ 10000046 │ 10000071 │
    │ 22  │ 10000022 │ 10000047 │ 10000072 │
    │ 23  │ 10000023 │ 10000048 │ 10000073 │
    │ 24  │ 10000024 │ 10000049 │ 10000074 │
    │ 25  │ 10000025 │ 10000050 │ 10000075 │"""
end

@testset "IOContext parameters test" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])
    str1, size = capture_stdout() do
        show(df)
    end
    io = IOContext(IOBuffer(), :limit=>true, :displaysize=>size)
    show(io, df)
    str2 = String(take!(io.io))
    @test str1 == str2

    Random.seed!(1)
    df_big = DataFrame(rand(25,5))
    str1, size = capture_stdout() do
        show(df_big)
    end
    io = IOContext(IOBuffer(), :limit=>true, :displaysize=>size)
    show(io, df_big)
    str2 = String(take!(io.io))
    @test str1 == str2
end

@testset "SubDataFrame show test" begin
    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])
    subdf = view(df, [2, 3], :)
    io = IOBuffer()
    show(io, subdf, allrows=true, allcols=false)
    str = String(take!(io))
    @test str == """
    2×3 SubDataFrame
    │ Row │ A     │ B           │ C       │
    │     │ Int64 │ String      │ Float32 │
    ├─────┼───────┼─────────────┼─────────┤
    │ 1   │ 2     │ ∀ε>0: x+ε>x │ 2.0     │
    │ 2   │ 3     │ z\$          │ 3.0     │"""
    show(io, subdf, allrows=true)
    show(io, subdf, allcols=true)
    show(io, subdf, allcols=true, allrows=true)
end

@testset "Test showing StackedVector and RepeatedVector" begin
    A = DataFrames.StackedVector(Any[[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    @test sprint(show, A) == "[1, 2, 3, 4, 5, 6, 7, 8, 9]"
    A = DataFrames.RepeatedVector([1, 2, 3], 5, 1)
    @test sprint(show, A) == "[1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]"
    A = DataFrames.RepeatedVector([1, 2, 3], 1, 5)
    @test sprint(show, A) == "[1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]"
end

@testset "Test colors and non-standard values: missing and nothing" begin
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, missing])
    @test sprint(show, df, context=:color=>true) == """
        2×2 DataFrame
        │ Row │ Fish   │ Mass     │
        │     │ \e[90mString\e[39m │ \e[90mFloat64?\e[39m │
        ├─────┼────────┼──────────┤
        │ 1   │ Suzy   │ 1.5      │
        │ 2   │ Amir   │ \e[90mmissing\e[39m  │"""

    df = DataFrame(A = [:Symbol, missing, :missing],
                   B = [missing, "String", "missing"],
                   C = [:missing, "missing", missing])
    @test sprint(show, df, context=:color=>true) == """
        3×3 DataFrame
        │ Row │ A       │ B       │ C       │
        │     │ \e[90mSymbol?\e[39m │ \e[90mString?\e[39m │ \e[90mAny\e[39m     │
        ├─────┼─────────┼─────────┼─────────┤
        │ 1   │ Symbol  │ \e[90mmissing\e[39m │ missing │
        │ 2   │ \e[90mmissing\e[39m │ String  │ missing │
        │ 3   │ missing │ missing │ \e[90mmissing\e[39m │"""

    df_nothing = DataFrame(A = [1.0, 2.0, 3.0], B = ["g", "g", nothing])
    @test sprint(show, df_nothing) == """
    3×2 DataFrame
    │ Row │ A       │ B      │
    │     │ Float64 │ Union… │
    ├─────┼─────────┼────────┤
    │ 1   │ 1.0     │ g      │
    │ 2   │ 2.0     │ g      │
    │ 3   │ 3.0     │        │"""
end

@testset "Test correct width computation" begin
    df = DataFrame([["a"]], [:x])
    @test sprint(show, df) == """
    1×1 DataFrame
    │ Row │ x      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ a      │"""
end

@testset "Test showing special types" begin
    # strings with escapes
    df = DataFrame(a = ["1\n1", "2\t2", "3\r3", "4\$4", "5\"5", "6\\6"])
    @test sprint(show, df) == """
    6×1 DataFrame
    │ Row │ a      │
    │     │ String │
    ├─────┼────────┤
    │ 1   │ 1\\n1   │
    │ 2   │ 2\\t2   │
    │ 3   │ 3\\r3   │
    │ 4   │ 4\$4    │
    │ 5   │ 5"5    │
    │ 6   │ 6\\\\6   │"""

    # categorical
    df = DataFrame(a = categorical([1,2,3]), b = categorical(["a", "b", missing]))
    @test sprint(show, df) == """
    3×2 DataFrame
    │ Row │ a    │ b       │
    │     │ Cat… │ Cat…?   │
    ├─────┼──────┼─────────┤
    │ 1   │ 1    │ a       │
    │ 2   │ 2    │ b       │
    │ 3   │ 3    │ missing │"""

    # BigFloat
    df = DataFrame(a = [big(1.0), missing])
    @test sprint(show, df) == """
    2×1 DataFrame
    │ Row │ a         │
    │     │ BigFloat? │
    ├─────┼───────────┤
    │ 1   │ 1.0       │
    │ 2   │ missing   │"""

    # date types
    df = DataFrame(a = Date(2020, 2, 11), b = DateTime(2020, 2, 11, 15), c = Day(1))
    @test sprint(show, df) == """
    1×3 DataFrame
    │ Row │ a          │ b                   │ c     │
    │     │ Date       │ DateTime            │ Day   │
    ├─────┼────────────┼─────────────────────┼───────┤
    │ 1   │ 2020-02-11 │ 2020-02-11T15:00:00 │ 1 day │"""

    # Irrational
    df = DataFrame(a = π)
    if VERSION < v"1.2.0-DEV.276"
        @test sprint(show, df) == """
        1×1 DataFrame
        │ Row │ a                      │
        │     │ Irrational{:π}         │
        ├─────┼────────────────────────┤
        │ 1   │ π = 3.1415926535897... │"""
    else
        @test sprint(show, df) == """
        1×1 DataFrame
        │ Row │ a        │
        │     │ Irratio… │
        ├─────┼──────────┤
        │ 1   │ π        │"""
    end
end

@testset "Test using :compact parameter of IOContext" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, df) == """
        1×1 DataFrame
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint(show, df, context=:compact=>false) == """
        1×1 DataFrame
        │ Row │ x                 │
        │     │ Float64           │
        ├─────┼───────────────────┤
        │ 1   │ 3.141592653589793 │"""
end

@testset "Test of DataFrameRows and DataFrameColumns" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, eachrow(df)) == """
        1×1 DataFrameRows
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint((io, x) -> show(io, x, summary=false), eachrow(df)) == """

        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint(show, eachcol(df)) == """
        1×1 DataFrameColumns
        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""

    @test sprint((io, x) -> show(io, x, summary=false), eachcol(df)) == """

        │ Row │ x       │
        │     │ Float64 │
        ├─────┼─────────┤
        │ 1   │ 3.14159 │"""
end

@testset "Test empty data frame and DataFrameRow" begin
    df = DataFrame(x = [float(pi)])
    @test sprint(show, df[:, 2:1]) == "0×0 DataFrame\n"
    @test sprint(show, @view df[:, 2:1]) == "0×0 SubDataFrame\n"
    @test sprint(show, df[1, 2:1]) == "DataFrameRow"
end

@testset "consistency" begin
    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(df.c, 5)
    @test_throws AssertionError sprint(show, df)

    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(DataFrames._columns(df), df[:, :a])
    @test_throws AssertionError sprint(show, df)
end

@testset "wide type name" begin
    @test sprint(show, DataFrame(a=⛵⛵⛵⛵⛵())) == """
    1×1 DataFrame
    │ Row │ a       │
    │     │ ⛵⛵⛵… │
    ├─────┼─────────┤
    │ 1   │ "⛵"    │"""

    @test sprint(show, DataFrame(a=categorical([Int64(2)^54]))) == """
    1×1 DataFrame
    │ Row │ a                 │
    │     │ CategoricalValue… │
    ├─────┼───────────────────┤
    │ 1   │ 18014398509481984 │"""

    @test sprint(show, DataFrame(a=categorical([Int64(2)^53]))) == """
    1×1 DataFrame
    │ Row │ a                │
    │     │ Categorical…     │
    ├─────┼──────────────────┤
    │ 1   │ 9007199254740992 │"""

    @test sprint(show, DataFrame(a=categorical([Int64(2)^37]))) == """
    1×1 DataFrame
    │ Row │ a            │
    │     │ Categorical… │
    ├─────┼──────────────┤
    │ 1   │ 137438953472 │"""

    @test sprint(show, DataFrame(a=categorical([Int64(2)^36]))) == """
    1×1 DataFrame
    │ Row │ a           │
    │     │ Cat…        │
    ├─────┼─────────────┤
    │ 1   │ 68719476736 │"""

    @test sprint(show, DataFrame(a=Union{Function,Missing}[missing])) == """
    1×1 DataFrame
    │ Row │ a         │
    │     │ Function? │
    ├─────┼───────────┤
    │ 1   │ missing   │"""
end

@testset "wide type name" begin
    df = DataFrame(A = Int32.(1:3), B = ["x", "y", "z"])

    io = IOBuffer()
    show(io, df, eltypes=true)
    str = String(take!(io))
    @test str == """
    3×2 DataFrame
    │ Row │ A     │ B      │
    │     │ Int32 │ String │
    ├─────┼───────┼────────┤
    │ 1   │ 1     │ x      │
    │ 2   │ 2     │ y      │
    │ 3   │ 3     │ z      │"""

    io = IOBuffer()
    show(io, df, eltypes=false)
    str = String(take!(io))
    @test str == """
    3×2 DataFrame
    │ Row │ A │ B │
    ├─────┼───┼───┤
    │ 1   │ 1 │ x │
    │ 2   │ 2 │ y │
    │ 3   │ 3 │ z │"""
end

end # module


text:
using Distributed
using Random

import NODAL

addprocs()

@everywhere begin
    using NODAL

    function tour_cost(x::Configuration, parameters::Dict{Symbol, Any})
        result = parse(Float64, read(`./tour_cost $(x["Tour"].value)`, String))
        result
    end

end

function tsp()
    !isfile("tour_cost") && run(`make`)

    tour = ["1"]
    for i = 2:48
        push!(tour, string(i))
    end
    shuffle!(tour)

    configuration = Configuration([PermutationParameter(tour ,"Tour")],
                                   "TSP Solution")

    tuning_run = Run(cost                = tour_cost,
                     starting_point      = configuration,
                     stopping_criterion  = elapsed_time_criterion,
                     duration            = 60,
                     reporting_criterion = elapsed_time_reporting_criterion,
                     report_after        = 5,
                     methods             = [[:iterative_first_improvement 2];
                                            [:iterative_greedy_construction 2];
                                            [:iterative_probabilistic_improvement 2];
                                            [:randomized_first_improvement 2];
                                            [:simulated_annealing 2];])

    @spawn optimize(tuning_run)
    result = take!(tuning_run.channel)

    print("$(result.current_time) ")
    println(result.cost_minimum)

    while !result.is_final
        result = take!(tuning_run.channel)
        print("$(result.current_time) ")
        println(result.cost_minimum)
    end
end

tsp()


text:
using DataFrames, Dates
import FileIO:save
# Tools for runners

function runner_eachlearner(learnerslist,
                    specs::Union{SpecsMultinomial, SpecsGaussian},
                    foldername::String, filename::String)
    for i in 1:length(learnerslist)
        println("learnersname" * " : " * learnerslist[i].name)
        runner_core([learnerslist[i]], specs, foldername, filename)
    end
    println("Finishing: " *Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS"))
end

function runner_core(learnerslist, specs::Union{SpecsMultinomial, SpecsGaussian},
                    foldername::String, filename::String)

    savepath = makehomesavepath(foldername)
    mkdir(savepath)

    if typeof(specs) == SpecsMultinomial
        settingslearner = collect(Iterators.product(collect(learnerslist),
                            specs.stochasticities, specs.changeprobabilities))
    elseif typeof(specs) == SpecsGaussian
        settingslearner = collect(Iterators.product(collect(learnerslist),
                            specs.sigmas, specs.changeprobabilities))
    end
    settings = collect(Iterators.product(settingslearner, 1:specs.nseeds))
    df = comparison(settings, specs)

    save(joinpath(savepath, filename * learnerslist[1].name * ".jld2"),
            Dict("df" => df, "specs" => specs))
end


text:
export mutate_weight, mutate_connect, mutate_neuron, mutate_enabled, mutate

"Mutate the weight of genes"
function mutate_weight(ind::NEATIndiv, cfg::Dict)
    #TODO Add mutation power
    # https://github.com/FernandoTorres/NEAT/blob/master/src/genome.cpp
    ind_mut = NEATIndiv(ind)
    for g in values(ind_mut.genes)
        g.weight = g.weight + randn()*cfg["weight_factor"]
    end
    ind_mut.fitness.= -Inf
    ind_mut
end

function indexof(a::Array{Float64}, f::Float64)
    findall(x->x==f, a)[1]
end

"Add a connection between 2 random neurons"
function mutate_connect(ind::NEATIndiv, cfg::Dict)
    sort!(ind.neuron_pos)

    ind_mut = NEATIndiv(ind)
    nb_neur = length(ind_mut.neuron_pos)
    # println(nb_neur, " ", ind.neuron_pos)

    if cfg["hyperNEAT"] # HyperNEAT
        n_in = 4
        n_out = 1
    else # NEAT
        n_in = cfg["n_in"]
        n_out = cfg["n_out"]
    end

    # Valid neuron pairs
    valid = trues(nb_neur, nb_neur)

    # Remove existing ones
    for g in values(ind.genes)
        i_origin = indexof(ind.neuron_pos, g.origin)
        i_dest = indexof(ind.neuron_pos, g.destination)
        valid[i_origin, i_dest]=false
    end

    for orig in 1:nb_neur
        # Remove links to self
        valid[orig, orig]=false

        orig_pos = ind.neuron_pos[orig]

        for dest in 1:nb_neur
            dest_pos = ind.neuron_pos[dest]

            # Remove links towards input neurons and bias neuron
            if dest_pos <= 0
                valid[orig, dest]=false
            end

            # Remove links between output neurons (would not support adding a neuron)
            if orig_pos >= 1 && dest_pos >= 1
                valid[orig, dest]=false
            end

            # Remove recurrence if needed
            if !cfg["allow_recurrence"] && orig_pos > dest_pos
                valid[orig, dest]=false
            end
        end
    end

    ind_mut.fitness.= -Inf
    ind_mut.network = Network(n_in, n_out, Dict()) # Reset network

    # Filter invalid ones
    conns = findall(valid)
    if length(conns) > 0
        shuffle!(conns) # Pick random
        cfg["innovation_max"] += 1

        i = ind.neuron_pos[conns[1][1]]
        j = ind.neuron_pos[conns[1][2]]
        if j <= 0 # Catching error where destination is an input
            return ind_mut # Refuse mutation
        end

        g = Gene(i, j, cfg["innovation_max"])

        ind_mut.genes[cfg["innovation_max"]]=g
    end
    ind_mut
end

"Remove a random connection"
function mutate_disconnect(ind::NEATIndiv, cfg::Dict)
    ind_mut = NEATIndiv(ind)
    if length(ind_mut.genes)<=1 # Always keep 1 gene
        return ind_mut
    end
    k = rand(collect(keys(ind_mut.genes))) # pick a random gene
    pop!(ind_mut.genes, k) # remove it
    if cfg["hyperNEAT"] # HyperNEAT
        n_in = 4
        n_out = 1
    else # NEAT
        n_in = cfg["n_in"]
        n_out = cfg["n_out"]
    end
    ind_mut.network = Network(n_in, n_out, Dict()) # Reset network
    ind_mut
end

"Split a connection into 2 connections with a new neuron"
function mutate_neuron(ind::NEATIndiv, cfg::Dict)
    ind_mut = NEATIndiv(ind)

    if length(ind.genes)==0
        return ind_mut
    end

    # Connection to split
    g = rand(collect(values(ind_mut.genes)))
    g.enabled = false

    # Create neuron between origin and destination, in [0; 1]
    n = random_position(g.origin, g.destination)
    push!(ind_mut.neuron_pos, n)
    sort!(ind_mut.neuron_pos)
    ind_mut.activ_functions[n] = rand(cfg["activation_functions"])

    # Create connections
    i = cfg["innovation_max"]
    ind_mut.genes[i+1] = Gene(g.origin, n, i+1)
    ind_mut.genes[i+2] = Gene(n, g.destination, i+2)
    cfg["innovation_max"] +=2
    ind_mut.fitness.= -Inf
    ind_mut
end

"Switch enabled"
function mutate_enabled(ind::NEATIndiv, cfg::Dict)
    ind_mut = NEATIndiv(ind)
    g = rand(collect(values(ind_mut.genes)))
    g.enabled = !g.enabled
    ind_mut
end

function mutate(ind::NEATIndiv, cfg::Dict)
    # println("mutation")
    if rand() < cfg["p_mutate_add_neuron"]
        c =  mutate_neuron(ind, cfg)
    elseif rand() < cfg["p_mutate_add_connection"]
        c = mutate_connect(ind, cfg)
    elseif rand() < cfg["p_mutate_remove_connection"]
        c = mutate_disconnect(ind, cfg)
    elseif rand() < cfg["p_mutate_weights"]
        c = mutate_weight(ind, cfg)
    elseif rand() < cfg["p_mutate_enabled"]
        c = mutate_enabled(ind, cfg)
    else
        # return clone if no mutation occurs
        c = NEATIndiv(ind)
    end
    build!(c)
    c
end


text:
# This file was generated, do not modify it. # hide
@pipeline HotTreeClf(hot = OneHotEncoder(),
                     tree = DecisionTreeClassifier()) is_probabilistic=true

mdl = HotTreeClf()
mach = machine(mdl, X, y)
fit!(mach);


text:
module SupplyChainOptimization

include("Modeling.jl")
include("Querying.jl")
include("Visualization.jl")

using Base: Bool, product
using JuMP
using HiGHS

export SupplyChain,
      Product,
      Customer,
      Storage,
      Supplier,
      Plant,
      Lane,
      Location,
      Demand,
      Node,
      add_lane!,
      add_product!,
      add_customer!,
      add_storage!,
      add_supplier!,
      add_demand!,
      add_plant!,
      optimize_network!,
      get_total_costs,
      get_total_fixed_costs,
      get_total_transportation_costs,
      get_production,
      get_shipments,
      get_receipts,
      get_inventory_at_start,
      get_inventory_at_end,
      is_opened,
      is_opening,
      is_closing,
      haversine,
      plot_flows,
      plot_costs,
      plot_network,
      plot_inventory

function get_demand(supply_chain, customer, product, time)
    for demand in supply_chain.demand
        if demand.customer == customer && demand.product == product
            return demand.demand[time]
        end
    end
    return 0
end

function get_service_level(supply_chain, customer, product)
    for demand in supply_chain.demand
        if demand.customer == customer && demand.product == product
            return demand.service_level
        end
    end
    return 1.0
end

function get_lanes_in(supply_chain, node)
    if(haskey(supply_chain.lanes_in, node))
        return supply_chain.lanes_in[node]
    else
        return Lane[]
    end
end

function get_lanes_out(supply_chain, node)
    if(haskey(supply_chain.lanes_out, node))
        return supply_chain.lanes_out[node]
    end
    return Lane[]
end

function get_bom(production, output, input)
    if(haskey(production.bill_of_material, output))
        if(haskey(production.bill_of_material[output], input))
            return production.bill_of_material[output][input]
        end
    end
    return 0
end

function get_maximum_throughput(node, product)
    if(haskey(node.maximum_throughput, product))
        return node.maximum_throughput[product]
    else
        return 0
    end
end

function get_additional_stock_cover(node, product)
    if(haskey(node.additional_stock_cover, product))
        return node.additional_stock_cover[product]
    else
        return 0
    end
end

function check_model(supply_chain)
    for production in supply_chain.plants
        for product in supply_chain.products
            if (haskey(production.bill_of_material, product) && !haskey(production.unit_cost, product)) ||
            (!haskey(production.bill_of_material, product) && haskey(production.unit_cost, product)) ||
            (haskey(production.bill_of_material, product) && !haskey(production.time, product)) ||
            (!haskey(production.bill_of_material, product) && haskey(production.time, product))
                throw(ArgumentError("Production $production must have the same produts in its bill_of_material, its unit_cost and its time data."))
            end
        end
    end
end

"""
    optimize_network!(supply_chain::SupplyChain, optimizer=HiGHS.Optimizer)

Optimizes the supply chain.
"""
function optimize_network!(supply_chain, optimizer=HiGHS.Optimizer)
    create_network_optimization_model!(supply_chain, optimizer)
    set_optimizer_attribute(supply_chain.optimization_model, "mip_heuristic_effort", 0.35)
    optimize_network_optimization_model!(supply_chain)
end

"""
Creates an optimization model.
"""
function create_network_optimization_model!(supply_chain, optimizer)
    supply_chain.optimization_model = create_network_optimization_model(supply_chain, optimizer)
    set_optimizer_attribute(supply_chain.optimization_model, "primal_feasibility_tolerance", 1e-5)
end

"""
Optimizes an optimization model.
"""
function optimize_network_optimization_model!(supply_chain)
    JuMP.optimize!(supply_chain.optimization_model)
end

"""
Creates an optimization model.
"""
function create_network_optimization_model(supply_chain, optimizer, bigM=100_000)
    check_model(supply_chain)

    times = 1:supply_chain.horizon
    products = supply_chain.products
    customers = supply_chain.customers
    storages = supply_chain.storages
    suppliers = supply_chain.suppliers
    plants = supply_chain.plants
    plants_storages = [x for x in union(plants, storages)]
    lanes = supply_chain.lanes

    m = Model(optimizer)#; bridge_constraints = false)

    @variable(m, total_costs >= 0)
    @variable(m, total_transportation_costs >= 0)
    @variable(m, total_fixed_costs >= 0)
    @variable(m, total_holding_costs >= 0)

    @variable(m, total_costs_per_period[times] >= 0)
    @variable(m, total_transportation_costs_per_period[times] >= 0)
    @variable(m, total_fixed_costs_per_period[times] >= 0)
    @variable(m, total_holding_costs_per_period[times] >= 0)

    @variable(m, opened[plants_storages, times], Bin)
    @variable(m, opening[plants_storages, times], Bin)
    @variable(m, closing[plants_storages, times], Bin)

    #@variable(m, serviced_by[products, storages, customers, times], Bin)

    @variable(m, lost_sales[products, customers, times] >= 0)

    @variable(m, bought[products, suppliers, times] >= 0)

    @variable(m, produced[products, plants, times] >= 0)

    @variable(m, stored_at_start[products, storages, times] >= 0)
    @variable(m, stored_at_end[products, storages, times] >= 0)

    @variable(m, used[lanes, times], Bin)
    @variable(m, sent[products, lanes, times] >= 0)
    @variable(m, received[products, lanes, times] >= 0)

    #single source constraint
    #@constraint(m, [p=products, c=customers, t=times], sum(serviced_by[p, s, c, t] for s in storages) <= 1)
    #@constraint(m, [p=products, s=storages, c=customers, t=times], sum(received[p, l, t] for l in get_lanes_in(supply_chain, c) if l.origin == s) <= bigM * serviced_by[p, s, c, t])

    @constraint(m, [p=products, s=storages; haskey(s.initial_inventory, p)], stored_at_start[p, s, 1] == s.initial_inventory[p])
    @constraint(m, [p=products, s=storages; !haskey(s.initial_inventory, p)], stored_at_start[p, s, 1] <= stored_at_end[p, s, supply_chain.horizon])

    @constraint(m, [p=products, l=lanes, t=times; t > l.time], received[p, l, t] == sent[p, l, t - l.time])
    @constraint(m, [p=products, l=lanes, t=times; t <= l.time], received[p, l, t] == 0)

    @constraint(m, [l=lanes, t=times; l.minimum_quantity > 0], sum(sent[p, l, t] for p in products) <= bigM * used[l, t])
    @constraint(m, [l=lanes, t=times; l.minimum_quantity > 0], sum(sent[p, l, t] for p in products) >= l.minimum_quantity * used[l, t])

    #@constraint(m, [s=storages, t=times], !opened[s, t] => { sum(sent[p, l, t] for p in products, l in get_lanes_out(supply_chain, s)) == 0 })
    @constraint(m, [s=storages, t=times], sum(sent[p, l, t] for p in products, l in get_lanes_out(supply_chain, s)) <= bigM * opened[s, t])
    @constraint(m, [p=products, s=storages, t=times; !isinf(get_maximum_throughput(s, p))], sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s)) <= get_maximum_throughput(s, p))
    #@constraint(m, [s=storages, t=times], !opened[s, t] => { sum(received[p, l, t] for p in products, l in get_lanes_in(supply_chain, s)) == 0 })
    @constraint(m, [s=storages, t=times], sum(received[p, l, t] for p in products, l in get_lanes_in(supply_chain, s)) <= bigM * opened[s, t])

    @constraint(m, [s=plants_storages], opening[s, 1] >= opened[s, 1] + (1 - s.initial_opened) - 1)
    @constraint(m, [s=plants_storages], opening[s, 1] <= opened[s, 1])
    @constraint(m, [s=plants_storages], opening[s, 1] <= 1 - s.initial_opened)

    @constraint(m, [s=plants_storages, t=times; t > 1], opening[s, t] >= opened[s, t] + (1 - opened[s, t-1]) - 1)
    @constraint(m, [s=plants_storages, t=times; t > 1], opening[s, t] <= opened[s, t])
    @constraint(m, [s=plants_storages, t=times; t > 1], opening[s, t] <= 1 - opened[s, t-1])

    @constraint(m, [s=plants_storages], closing[s, 1] >= (1 - opened[s, 1]) + s.initial_opened - 1)
    @constraint(m, [s=plants_storages], closing[s, 1] <= 1 - opened[s, 1])
    @constraint(m, [s=plants_storages], closing[s, 1] <= s.initial_opened)

    @constraint(m, [s=plants_storages, t=times; t > 1], closing[s, t] >= (1 - opened[s, t]) + opened[s, t-1] - 1)
    @constraint(m, [s=plants_storages, t=times; t > 1], closing[s, t] <= 1 - opened[s, t])
    @constraint(m, [s=plants_storages, t=times; t > 1], closing[s, t] <= opened[s, t-1])

    @constraint(m, [s=plants_storages, t=times; isinf(s.opening_cost)], opending[s, t] == 0)
    @constraint(m, [s=plants_storages, t=times; isinf(s.closing_cost)], closing[s, t] == 0)

    @constraint(m, [p=products, s=storages, t=times], stored_at_end[p, s, t] == stored_at_start[p, s, t]
                                                                            + sum(received[p, l, t] for l in get_lanes_in(supply_chain, s))
                                                                            - sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s))
                                                                            )
    @constraint(m, [p=products, s=storages, t=times; t > 1], stored_at_start[p, s, t] == stored_at_end[p, s, t - 1])
    @constraint(m, [p=products, s=storages, t=times], stored_at_end[p, s, t] >= get_additional_stock_cover(s, p) * sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s)))

    @constraint(m, [p=products, s=suppliers, t=times], bought[p, s, t] == sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s)))
    @constraint(m, [p=products, s=suppliers, t=times; !isinf(get_maximum_throughput(s, p))], sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s)) <= get_maximum_throughput(s, p))

    #@constraint(m, [p=products, s=plants, t=times], closed[s, t] => { produced[p, s, t] == 0 })
    for s in plants, p in products
        if haskey(s.time, p)
            @constraint(m, [t=times, ti=t:min(t+s.time[p], supply_chain.horizon)], produced[p, s, t] <= bigM * opened[s, ti])
        else
            @constraint(m, [t=times], produced[p, s, t] == 0)
        end
    end
    @constraint(m, [p=products, s=plants, t=times; haskey(s.time, p) && (t + s.time[p] <= supply_chain.horizon)], produced[p, s, t] == sum(sent[p, l, t + s.time[p]] for l in get_lanes_out(supply_chain, s)))
    @constraint(m, [p=products, s=plants, t=times; !isinf(get_maximum_throughput(s, p))], sum(sent[p, l, t] for l in get_lanes_out(supply_chain, s)) <= get_maximum_throughput(s, p))
    @constraint(m, [p=products, s=plants, t=times], sum(produced[p2, s, t] * get_bom(s, p2, p) for p2 in products) == sum(received[p, l, t] for l in get_lanes_in(supply_chain, s)))

    @constraint(m, [p=products, c=customers, t=times], sum(received[p, l, t] for l in get_lanes_in(supply_chain, c)) == get_demand(supply_chain, c, p, t) - lost_sales[p, c, t])

    @constraint(m, [p=products, c=customers], sum(lost_sales[p, c, t] for t in times) <= (1 - get_service_level(supply_chain, c, p)) * sum(get_demand(supply_chain, c, p, t) for t in times))

    @constraint(m, [t=times], total_transportation_costs_per_period[t] == sum(sent[p, l, t] * l.unit_cost for p in products, l in lanes))
    @constraint(m, total_transportation_costs == sum(total_transportation_costs_per_period[t] for t in times))

    @constraint(m, [t=times], total_fixed_costs_per_period[t] == sum(opened[s, t] * s.fixed_cost for s in plants_storages))
    @constraint(m, total_fixed_costs == sum(total_fixed_costs_per_period[t] for t in times))

    @constraint(m, [t=times], total_holding_costs_per_period[t] == sum(stored_at_end[p, s, t] * p.unit_holding_cost for p in products, s in storages))
    @constraint(m, total_holding_costs == sum(total_holding_costs_per_period[t] for t in times))

    @constraint(m, [t=times], total_costs_per_period[t] == total_transportation_costs_per_period[t] +
                       total_fixed_costs_per_period[t] +
                       sum(opening[s, t] * s.opening_cost for s in plants_storages if !isinf(s.opening_cost)) +
                       sum(closing[s, t] * s.closing_cost for s in plants_storages if !isinf(s.closing_cost)) +
                       sum(sum(received[p, l, t] * s.unit_handling_cost[p] for l in get_lanes_in(supply_chain, s)) for p in products for s in storages if haskey(s.unit_handling_cost, p)) +
                       sum(bought[p, s, t] * s.unit_cost[p] for p in products, s in suppliers if haskey(s.unit_cost, p)) +
                       sum(produced[p, s, t] * s.unit_cost[p] for p in products, s in plants if haskey(s.unit_cost, p)) +
                       total_holding_costs)

    @constraint(m, total_costs == sum(total_costs_per_period[t] for t in times))

    @objective(m, Min, 1.0 * total_costs)

    return m
end

"""
Computes the great circle distance between two locations. The distance is expressed in meter.
"""
function haversine(location1::Location, location2::Location)
    return haversine(location1.latitude, location1.longitude, location2.latitude, location2.longitude)
end

function haversine(lat1, lon1, lat2, lon2)
    R = 6371e3

    Δlat = lat2 - lat1
    Δlon = lon2 - lon1

    return 2 * R * asin(sqrt(sind(Δlat / 2) ^ 2 + cosd(lat1) * cosd(lat2) * sind(Δlon / 2) ^ 2))
end

end # module


text:
abstract type Simplexer end

struct AffineSimplexer <: Simplexer
    a::Float64
    b::Float64
end
AffineSimplexer(;a = 0.025, b = 0.5) = AffineSimplexer(a, b)

function simplexer(S::AffineSimplexer, initial_x::Tx) where Tx
    n = length(initial_x)
    initial_simplex = Tx[copy(initial_x) for i = 1:n+1]
    for j = 1:n
        initial_simplex[j+1][j] = (1+S.b) * initial_simplex[j+1][j] + S.a
    end
    initial_simplex
end

abstract type NMParameters end

struct AdaptiveParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

AdaptiveParameters(;  α = 1.0, β = 1.0, γ = 0.75 , δ = 1.0) = AdaptiveParameters(α, β, γ, δ)
parameters(P::AdaptiveParameters, n::Integer) = (P.α, P.β + 2/n, P.γ - 1/2n, P.δ - 1/n)

struct FixedParameters <: NMParameters
    α::Float64
    β::Float64
    γ::Float64
    δ::Float64
end

FixedParameters(; α = 1.0, β = 2.0, γ = 0.5, δ = 0.5) = FixedParameters(α, β, γ, δ)
parameters(P::FixedParameters, n::Integer) = (P.α, P.β, P.γ, P.δ)

struct NelderMead{Ts <: Simplexer, Tp <: NMParameters} <: ZerothOrderOptimizer
    initial_simplex::Ts
    parameters::Tp
end

"""
# NelderMead
## Constructor
```julia
NelderMead(; parameters = AdaptiveParameters(),
             initial_simplex = AffineSimplexer())
```

The constructor takes 2 keywords:

* `parameters`, an instance of either `AdaptiveParameters` or `FixedParameters`,
and is used to generate parameters for the Nelder-Mead Algorithm
* `initial_simplex`, an instance of `AffineSimplexer`

## Description
Our current implementation of the Nelder-Mead algorithm is based on [1] and [3].
Gradient-free methods can be a bit sensitive to starting values and tuning parameters,
so it is a good idea to be careful with the defaults provided in Optim.jl.

Instead of using gradient information, Nelder-Mead is a direct search method. It keeps
track of the function value at a number of points in the search space. Together, the
points form a simplex. Given a simplex, we can perform one of four actions: reflect,
expand, contract, or shrink. Basically, the goal is to iteratively replace the worst
point with a better point. More information can be found in [1], [2] or [3].

## References
- [1] Nelder, John A. and R. Mead (1965). "A simplex method for function minimization". Computer Journal 7: 308–313. doi:10.1093/comjnl/7.4.308
- [2] Lagarias, Jeffrey C., et al. "Convergence properties of the Nelder–Mead simplex method in low dimensions." SIAM Journal on Optimization 9.1 (1998): 112-147
- [3] Gao, Fuchang and Lixing Han (2010). "Implementing the Nelder-Mead simplex algorithm with adaptive parameters". Computational Optimization and Applications. doi:10.1007/s10589-010-9329-3
"""
function NelderMead(; kwargs...)
    KW = Dict(kwargs)
    if haskey(KW, :initial_simplex) || haskey(KW, :parameters)
        initial_simplex, parameters = AffineSimplexer(), AdaptiveParameters()
        haskey(KW, :initial_simplex) && (initial_simplex = KW[:initial_simplex])
        haskey(KW, :parameters) && (parameters = KW[:parameters])
        return NelderMead(initial_simplex, parameters)
    else
        return NelderMead(AffineSimplexer(), AdaptiveParameters())
    end
end

Base.summary(::NelderMead) = "Nelder-Mead"

# centroid except h-th vertex
function centroid!(c::AbstractArray{T}, simplex, h=0) where T
    n = length(c)
    fill!(c, zero(T))
    @inbounds for i in 1:n+1
        if i != h
            xi = simplex[i]
            c .+= xi
        end
    end
    rmul!(c, T(1)/n)
end

centroid(simplex, h) = centroid!(similar(simplex[1]), simplex, h)

nmobjective(y::Vector, m::Integer, n::Integer) = sqrt(var(y) * (m / n))

function print_header(method::NelderMead)
    @printf "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
    @printf "------   --------------    --------------\n"
end

function Base.show(io::IO, trace::OptimizationTrace{T, NelderMead}) where T
    @printf io "Iter     Function value    √(Σ(yᵢ-ȳ)²)/n \n"
    @printf io "------   --------------    --------------\n"
    for state in trace.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState{NelderMead})
    @printf io "%6d   %14e    %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

mutable struct NelderMeadState{Tx, T, Tfs} <: ZerothOrderState
    x::Tx
    m::Int
    simplex::Vector{Tx}
    x_centroid::Tx
    x_lowest::Tx
    x_second_highest::Tx
    x_highest::Tx
    x_reflect::Tx
    x_cache::Tx
    f_simplex::Tfs
    nm_x::T
    f_lowest::T
    i_order::Vector{Int}
    α::T
    β::T
    γ::T
    δ::T
    step_type::String
end

function initial_state(method::NelderMead, options, d, initial_x)
    T = eltype(initial_x)
    n = length(initial_x)
    m = n + 1
    simplex = simplexer(method.initial_simplex, initial_x)
    f_simplex = zeros(T, m)

    value!!(d, first(simplex))

    @inbounds for i in 1:length(simplex)
        f_simplex[i] = value(d, simplex[i])
    end
    # Get the indices that correspond to the ordering of the f values
    # at the vertices. i_order[1] is the index in the simplex of the vertex
    # with the lowest function value, and i_order[end] is the index in the
    # simplex of the vertex with the highest function value
    i_order = sortperm(f_simplex)

    α, β, γ, δ = parameters(method.parameters, n)

NelderMeadState(similar(initial_x), # Variable to hold final minimizer value for MultivariateOptimizationResults
          m, # Number of vertices in the simplex
          simplex, # Maintain simplex in state.simplex
          centroid(simplex,  i_order[m]), # Maintain centroid in state.centroid
          similar(initial_x), # Store cache in state.x_lowest
          similar(initial_x), # Store cache in state.x_second_highest
          similar(initial_x), # Store cache in state.x_highest
          similar(initial_x), # Store cache in state.x_reflect
          similar(initial_x), # Store cache in state.x_cache
          f_simplex, # Store objective values at the vertices in state.f_simplex
          T(nmobjective(f_simplex, n, m)), # Store nmobjective in state.nm_x
          f_simplex[i_order[1]], # Store lowest f in state.f_lowest
          i_order, # Store a vector of rankings of objective values
          T(α),
          T(β),
          T(γ),
          T(δ),
          "initial")
end

function update_state!(f::F, state::NelderMeadState{T}, method::NelderMead) where {F, T}
    # Augment the iteration counter
    shrink = false
    n, m = length(state.x), state.m

    centroid!(state.x_centroid, state.simplex, state.i_order[m])
    copyto!(state.x_lowest, state.simplex[state.i_order[1]])
    copyto!(state.x_second_highest, state.simplex[state.i_order[n]])
    copyto!(state.x_highest, state.simplex[state.i_order[m]])
    state.f_lowest = state.f_simplex[state.i_order[1]]
    f_second_highest = state.f_simplex[state.i_order[n]]
    f_highest = state.f_simplex[state.i_order[m]]

    # Compute a reflection
    @inbounds for j in 1:n
        state.x_reflect[j] = state.x_centroid[j] + state.α * (state.x_centroid[j]-state.x_highest[j])
    end

    f_reflect = value(f, state.x_reflect)
    if f_reflect < state.f_lowest
        # Compute an expansion
        @inbounds for j in 1:n
            state.x_cache[j] = state.x_centroid[j] + state.β *(state.x_reflect[j] - state.x_centroid[j])
        end
        f_expand = value(f, state.x_cache)

        if f_expand < f_reflect
            copyto!(state.simplex[state.i_order[m]], state.x_cache)
            @inbounds state.f_simplex[state.i_order[m]] = f_expand
            state.step_type = "expansion"
        else
            copyto!(state.simplex[state.i_order[m]], state.x_reflect)
            @inbounds state.f_simplex[state.i_order[m]] = f_reflect
            state.step_type = "reflection"
        end
        # shift all order indeces, and wrap the last one around to the first
        i_highest = state.i_order[m]
        @inbounds for i = m:-1:2
            state.i_order[i] = state.i_order[i-1]
        end
        state.i_order[1] = i_highest
    elseif f_reflect < f_second_highest
        copyto!(state.simplex[state.i_order[m]], state.x_reflect)
        @inbounds state.f_simplex[state.i_order[m]] = f_reflect
        state.step_type = "reflection"
        sortperm!(state.i_order, state.f_simplex)
    else
        if f_reflect < f_highest
            # Outside contraction
            @simd for j in 1:n
                @inbounds state.x_cache[j] = state.x_centroid[j] + state.γ * (state.x_reflect[j]-state.x_centroid[j])
            end
            f_outside_contraction = value(f, state.x_cache)
            if f_outside_contraction < f_reflect
                copyto!(state.simplex[state.i_order[m]], state.x_cache)
                @inbounds state.f_simplex[state.i_order[m]] = f_outside_contraction
                state.step_type = "outside contraction"
                sortperm!(state.i_order, state.f_simplex)

            else
                shrink = true
            end
        else # f_reflect > f_highest
            # Inside constraction
            @simd for j in 1:n
                @inbounds state.x_cache[j] = state.x_centroid[j] - state.γ *(state.x_reflect[j] - state.x_centroid[j])
            end
            f_inside_contraction = value(f, state.x_cache)
            if f_inside_contraction < f_highest
                copyto!(state.simplex[state.i_order[m]], state.x_cache)
                @inbounds state.f_simplex[state.i_order[m]] = f_inside_contraction
                state.step_type = "inside contraction"
                sortperm!(state.i_order, state.f_simplex)
            else
                shrink = true
            end
        end
    end

    if shrink
        for i = 2:m
            ord = state.i_order[i]
            copyto!(state.simplex[ord], state.x_lowest + state.δ*(state.simplex[ord]-state.x_lowest))
            state.f_simplex[ord] = value(f, state.simplex[ord])
        end
        step_type = "shrink"
        sortperm!(state.i_order, state.f_simplex)
    end

    state.nm_x = nmobjective(state.f_simplex, n, m)
    false
end

function after_while!(f, state, method::NelderMead, options)
    sortperm!(state.i_order, state.f_simplex)
    x_centroid_min = centroid(state.simplex, state.i_order[state.m])
    f_centroid_min = value(f, x_centroid_min)
    f_min, i_f_min = findmin(state.f_simplex)
    x_min = state.simplex[i_f_min]
    if f_centroid_min < f_min
        x_min = x_centroid_min
        f_min = f_centroid_min
    end
    f.F = f_min
    state.x .= x_min
end
# We don't have an f_x_previous in NelderMeadState, so we need to special case these
pick_best_x(f_increased, state::NelderMeadState) = state.x
pick_best_f(f_increased, state::NelderMeadState, d) = value(d)

function assess_convergence(state::NelderMeadState, d, options)
    g_converged = state.nm_x <= options.g_tol # Hijact g_converged for NM stopping criterior
    return false, false, g_converged, g_converged, false
end

function initial_convergence(d, state::NelderMeadState, method::NelderMead, initial_x, options)
    nmobjective(state.f_simplex, state.m, length(initial_x)) < options.g_tol
end

function trace!(tr, d, state, iteration, method::NelderMead, options)
    dt = Dict()
    if options.extended_trace
        dt["centroid"] = copy(state.x_centroid)
        dt["step_type"] = state.step_type
    end
    update!(tr,
    iteration,
    state.f_lowest,
    state.nm_x,
    dt,
    options.store_trace,
    options.show_trace,
    options.show_every,
    options.callback)
end


text:
function partial_state_selection_graph!(state::TransformationState)
    find_solvables!(state; allow_symbolic=true)
    var_eq_matching = complete(pantelides!(state))
    complete!(state.structure)
    partial_state_selection_graph!(state.structure, var_eq_matching)
end

function ascend_dg(xs, dg, level)
    while level > 0
        xs = Int[dg[x] for x in xs]
        level -= 1
    end
    return xs
end

function ascend_dg_all(xs, dg, level, maxlevel)
    r = Int[]
    while true
        if level <= 0
            append!(r, xs)
        end
        maxlevel <= 0 && break
        xs = Int[dg[x] for x in xs if dg[x] !== nothing]
        level -= 1
        maxlevel -= 1
    end
    return r
end

function pss_graph_modia!(structure::SystemStructure, var_eq_matching, varlevel, inv_varlevel, inv_eqlevel)
    @unpack eq_to_diff, var_to_diff, graph, solvable_graph = structure

    # var_eq_matching is a maximal matching on the top-differentiated variables.
    # Find Strongly connected components. Note that after pantelides, we expect
    # a balanced system, so a maximal matching should be possible.
    var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
    var_eq_matching = Matching{Union{Unassigned, SelectedState}}(var_eq_matching)
    for vars in var_sccs
        # TODO: We should have a way to not have the scc code look at unassigned vars.
        if length(vars) == 1 && varlevel[vars[1]] != 0
            continue
        end

        # Now proceed level by level from lowest to highest and tear the graph.
        eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
        isempty(eqs) && continue
        maxlevel = level = maximum(map(x->inv_eqlevel[x], eqs))
        old_level_vars = ()
        ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph, complete(Matching(ndsts(graph)))); dir=:in)
        while level >= 0
            to_tear_eqs_toplevel = filter(eq->inv_eqlevel[eq] >= level, eqs)
            to_tear_eqs = ascend_dg(to_tear_eqs_toplevel, invview(eq_to_diff), level)

            to_tear_vars_toplevel = filter(var->inv_varlevel[var] >= level, vars)
            to_tear_vars = ascend_dg_all(to_tear_vars_toplevel, invview(var_to_diff), level, maxlevel)

            if old_level_vars !== ()
                # Inherit constraints from previous level.
                # TODO: Is this actually a good idea or do we want full freedom
                # to tear differently on each level? Does it make a difference
                # whether we're using heuristic or optimal tearing?
                removed_eqs = Int[]
                removed_vars = Int[]
                for var in old_level_vars
                    old_assign = ict.graph.matching[var]
                    if !isa(old_assign, Int) || ict.graph.matching[var_to_diff[var]] !== unassigned
                        continue
                    end
                    # Make sure the ict knows about this edge, so it doesn't accidentally introduce
                    # a cycle.
                    ok = try_assign_eq!(ict, var_to_diff[var], eq_to_diff[old_assign])
                    @assert ok
                    var_eq_matching[var_to_diff[var]] = eq_to_diff[old_assign]
                    push!(removed_eqs, eq_to_diff[ict.graph.matching[var]])
                    push!(removed_vars, var_to_diff[var])
                end
                to_tear_eqs = setdiff(to_tear_eqs, removed_eqs)
                to_tear_vars = setdiff(to_tear_vars, removed_vars)
            end
            filter!(var->ict.graph.matching[var] === unassigned, to_tear_vars)
            filter!(eq->invview(ict.graph.matching)[eq] === unassigned, to_tear_eqs)
            tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, to_tear_vars)
            for var in to_tear_vars
                var_eq_matching[var] = ict.graph.matching[var]
            end
            old_level_vars = to_tear_vars
            level -= 1
        end
    end
    for var in 1:ndsts(graph)
        if varlevel[var] !== 0 && var_eq_matching[var] === unassigned
            var_eq_matching[var] = SelectedState()
        end
    end
    return var_eq_matching
end

struct SelectedState; end
function partial_state_selection_graph!(structure::SystemStructure, var_eq_matching)
    @unpack eq_to_diff, var_to_diff, graph, solvable_graph = structure
    eq_to_diff = complete(eq_to_diff)

    inv_eqlevel = map(1:nsrcs(graph)) do eq
        level = 0
        while invview(eq_to_diff)[eq] !== nothing
            eq = invview(eq_to_diff)[eq]
            level += 1
        end
        level
    end

    varlevel = map(1:ndsts(graph)) do var
        level = 0
        while var_to_diff[var] !== nothing
            var = var_to_diff[var]
            level += 1
        end
        level
    end

    inv_varlevel = map(1:ndsts(graph)) do var
        level = 0
        while invview(var_to_diff)[var] !== nothing
            var = invview(var_to_diff)[var]
            level += 1
        end
        level
    end

    # TODO: Should pantelides just return this?
    for var in 1:ndsts(graph)
        if var_to_diff[var] !== nothing
            var_eq_matching[var] = unassigned
        end
    end

    var_eq_matching = pss_graph_modia!(structure,
        complete(var_eq_matching), varlevel, inv_varlevel, inv_eqlevel)

    var_eq_matching
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
#=
      Driver program to run the c_sw kernel

        This program acheives the goals of:
            * reading a NetCDF file
            * printing statistics on the input & output state in ./test/test_output
            * running the c_sw kernel located in sw_core.jl
            * creating a new NetCDF file with the data after the kernel runs

=#

module program

include("./sw_core.jl")

using TOML, Printf
using NCDatasets, OffsetArrays
using .sw_core_mod

export main

    function main(input::Dict)

        # define vars
        datasize = input["name"]
        nc_input_file = input["input_file"]
        print_output_file = input["test_out_file"]
        nc_output_file = input["output_file"]

    # Define the dataset using NCDatasets pkg from each datafile
    ds = NCDataset(nc_input_file, "a")

    # Open the export file with write/create/truncate privileges ("w")
    io = open(print_output_file,"w")

    # Assign Variables from the NetCDF file to a Julia Struct
    #
    # Could be written
    #
    # current_state = State(ds)
    # different_state = State(different_ds)
    assign_variables(ds)

    # Close the NetCDF dataset
    close(ds)

    # Print input state
    print_state("Input State", current_state, io)

    # Call the Julia kernel
    println("Run the kernel $datasize : ")
    @time for k = 1 : current_state.npz
        c_sw!(current_state, k)
    end

    # Print output state
    print_state("Output State", current_state, io)

    # Write a new NetCDF file
    println("Write NetCDF file $datasize : ")
    @time write_state(nc_output_file, current_state)

    # Close the IO file
    close(io)

    return true

    end # function main

end # module program