text:
module ChemometricsData
    using MD5, JSON3, MAT, CSV, DataFrames, Crayons,
            StringDistances, DataStructures, HTTP,
            ProgressBars, DataDeps

    global DATA_PATH = Base.joinpath( @__DIR__ , "..", "data" )

    include("manifest.jl")
    include("online_manifest.jl")

    include("Validation.jl")
    export flatten_dir
    include("PostTreatment.jl")
    export numeric_columns, nonnumeric_columns

    """
        suggest_a_dataset(dataset_name_lc)

    if a dataset isn't found try to guess what the user tried to type.
    Note: this is an internal function.
    """
    function suggest_a_dataset( dataset_name )
        dataset_name_lc = lowercase( dataset_name )
        println(Crayon(foreground = :red), "Sorry no such dataset found.")
        db_keys = collect(keys(data_manifest))
        if length( dataset_name ) > 1
            scores = [ Levenshtein()(dataset_name_lc, lowercase(k)) for k in db_keys ]
            println(Crayon(foreground = :yellow), "By any chance did you mean to type?")
            best_guess = first( sortperm(scores) )
            println(Crayon(foreground = :yellow, italics = true), "\t " * db_keys[best_guess])
        end
    end

    """
        search(keyword::String)::Vector{String}

    Enter a keyword, such as "Food", "Forensics", measurement technique, etc, and see what datasets which are available are related too it.
    """
    function search(keyword::String)::Vector{String}
        result::Vector{String} = []
        keyword_lc = lowercase(keyword)
        for ( dataset_name, meta_data ) in data_manifest
            search_regions = vcat(  get( meta_data, "usage", ""),
                                    get( meta_data, "keywords", ""),
                                    get( meta_data, "property values", ""),
                                    get( meta_data, "measurements", "")
                                    )
            if keyword_lc in lowercase.( search_regions )
                push!( result, dataset_name )
            end
        end
        return result
    end
    export search

    """
        describe( dataset_name::String )

    Prints a high level overview of the requested dataset if it exists to the console. If it does not exist it looks for a typo via the Levenshtein distance and informs the user.
    """
    function describe( dataset_name::String )
        dataset_name_lc = lowercase( dataset_name )
        if haskey( data_manifest, dataset_name )
            println(Crayon(foreground = :blue), "Dataset found.")
            sug_path = Base.joinpath( DATA_PATH , dataset_name )
            avail = isdir( sug_path ) ? "Offline" : "Online Only"
            println(Crayon(foreground = :green), "> Available $avail " )
            if avail == "Offline"
                println(Crayon(foreground = :green), "> Contains " * string(length(readdir(sug_path))) * " files." )
            else
                if isdir(sug_path)
                    println(Crayon(foreground = :green), "> Contains " * string(length(readdir(sug_path))) * " files." )
                else
                    println(Crayon(foreground = :green), "> Dataset must be downloaded with \'fetchdata()\' to use." )
                end
            end
            for (k,v) in data_manifest[dataset_name]
                if (v != "") && (v != [""])
                    println(Crayon(foreground = :green), "> $k : $v")
                end
            end
        else
            suggest_a_dataset(dataset_name)
        end
    end
    export describe

    function isCSV(x)
        return isfile(x) && (lowercase(x[(end-2):end]) == "csv")
    end

    """
        load( dataset_name::String )

    Loads the requested dataset as a dataframe into memory, or downloads it from the internet if it exists.
    """
    function load( dataset_name::String )
        dataset_name_lc = lowercase( dataset_name )
        if haskey( data_manifest, dataset_name )
            println(Crayon(foreground = :blue), "Dataset found.")
            sug_path = Base.joinpath( DATA_PATH , dataset_name )
            avail = isdir( sug_path ) ? "Offline" : "Online Only"
            if avail == "Offline"
                assets = readdir( sug_path )
                if length( readdir( sug_path ) ) > 1
                    data = Dict([ f => DataFrame!(CSV.File(Base.joinpath(sug_path, f); threaded=false )) for f in assets if isCSV(Base.joinpath(sug_path, f)) ])
                else
                    data = DataFrame!(CSV.File(Base.joinpath(sug_path, first(assets) ) ))
                end
                println(Crayon(foreground = :blue), "Dataset loaded.")
                if haskey( data_manifest[dataset_name], "Usage Statement" )
                    #did the author specify a special statement for its usage? if so display it!
                    println(Crayon(foreground = :red), data_manifest[dataset_name]["Usage Statement"])
                else
                    tmp = "Please honor the dataset author's contribution to our field by acknowledging their work: "
                    tmp = tmp .* ("\n" .* data_manifest[dataset_name]["references"])
                    println(Crayon(foreground = :green), first( tmp ) )
                end
                return data
            elseif haskey( online_manifest, dataset_name )
                println(Crayon(foreground = :red), "Dataset must be downloaded with the \"fetchdata()\" function.")
            else
                suggest_a_dataset(dataset_name)
            end
        end
        return nothing
    end
    export load

    """
        fetchdata( dataset_name::String )

    Will download a given dataset using information from `online_manifest.jl`.
    """
    function fetchdata( dataset_name::String )
        dataset_name_lc = lowercase( dataset_name )
        if haskey( data_manifest, dataset_name )
            println(Crayon(foreground = :blue), "Dataset found.")
            sug_path = Base.joinpath( DATA_PATH , dataset_name )
            avail = isdir( sug_path ) ? "Offline" : "Online Only"
            avail = "Online Only"
            if avail == "Offline"
                println(Crayon(foreground = :green), "Directory already exists, please delete it to redownload the dataset or use the \"load()\" function to bring it into memory.")
                println(Crayon(foreground = :green, italics = true), sug_path)
            elseif haskey( online_manifest, dataset_name )
                println(Crayon(foreground = :red), "> Creating Directory...")
                mkdir(sug_path) #works
                println(Crayon(foreground = :red), "> Downloading...\n" *
                        "from: ", online_manifest[dataset_name]["URL"], "\n" *
                        "to: ", sug_path)
                HTTP.download( online_manifest[dataset_name]["URL"], sug_path )
                compr_ext = [".zip", ".tar"]
                files = [ f for f in readdir(sug_path) if f[(end-3):end] in compr_ext ]
                (length(files) == 0) && @warn "Although a file was downloaded the file does not match the stored MD5 checksum. \n Please notify ChemometricsData.jl!"
                md5chk = [ f for f in files if check_MD5( Base.joinpath(sug_path, f), online_manifest[dataset_name]["MD5"] ) ]
                while length(files) > 0
                    file_of_interest = first( ( length( md5chk ) > 0 ) ? md5chk : files )#first( files )#
                    cd(sug_path) do #thanks Lyndon!
                        DataDeps.unpack( Base.joinpath( sug_path, file_of_interest ) )
                    end
                    flatten_dir(sug_path) #eeek this changes the filesystem don't screw up!
                    md5chk = []
                    files = [ f for f in readdir(sug_path) if f[(end-3):end] in compr_ext ]
                end
            else
                suggest_a_dataset(dataset_name)
            end
        end
        return nothing
    end
    export fetchdata

    """
        meta( dataset_name::String )

    prints the metadata associated with a dataset to the console. If the dataset does not exist it will print suggestions for similarly named datasets or keywords.
    """
    function meta( dataset_name::String )
        if haskey( data_manifest, dataset_name )
            return data_manifest[ dataset_name ]
        else
            suggest_a_dataset(dataset_name)
        end
        return nothing
    end
    export meta

end


text:
abstract MatchType
type Prefix <: MatchType end
type Glob <: MatchType end

abstract Verbosity
type Quiet   <: Verbosity end
type Normal  <: Verbosity end
type Verbose <: Verbosity end

abstract PrintExit
type Nothing <: PrintExit end
type List    <: PrintExit end
type Version <: PrintExit end

function default_config()
   {:ci => 0.95,
    :gc => true,
    :resamples  => 1000,
    :results => nothing,
    :report => nothing,
    :samples => 100,
    :summary => nothing,
    :template => "report.tpl",
    :quiet => false,
    :no_measurements => false,
    :junit => false}
end

function config(options)
   merge(default_config(), options)
end

ext: jl


text:
function iterate(fn::Function, initval::Number, k::Integer)
  x = sizehint([initval], k)
  for _ in 1:k-1
    push!(x, fn(x[end]))
  end
  x
end

function fixedpoint(x::Array{Float64, 1}; ϵ::Float64=1e-4)
  findfirst(d -> d < ϵ, abs(diff(x)))
end


text:
using CompScienceMeshes
using Test

fn = joinpath(dirname(@__FILE__),"FPLZX.msh")
m = CompScienceMeshes.load_gid_mesh(fn)

@test numcells(m) == 896
@test numvertices(m) == 450
@test vertices(m,1) ≈ [0,0,7]


text:
__precompile__()

module Shoco
    if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
        include("../deps/deps.jl")
    else
        error("Shoco not properly installed. Please run Pkg.build(\"Shoco\")",
              " and restart Julia.")
    end

    export compress, decompress

    if !isdefined(Base, :unsafe_string)
        const unsafe_string = bytestring
    end

    function compress(s::AbstractString)
        isempty(s) && return ""

        # The output should be no longer than the input
        compressed = Array(Cchar, sizeof(s))

        # The function modifies `compressed` and returns the number of bytes written
        nbytes = ccall((:shoco_compress, shoco), Int,
                       (Ptr{Cchar}, Csize_t, Ptr{Cchar}, Csize_t),
                       s, 0, compressed, sizeof(s))

        # Failed compression is shoco's fault, not ours >_>
        nbytes > 0 || throw(ErrorException("Compression failed for input $s"))

        compressed = compressed[1:nbytes]
        compressed[end] == Cchar(0) || push!(compressed, Cchar(0))

        return unsafe_string(pointer(compressed))
    end

    function decompress(s::AbstractString)
        isempty(s) && return ""

        # The decompressed string will be at most twice as long as the input
        decompressed = Array(Cchar, 2 * sizeof(s))

        nbytes = ccall((:shoco_decompress, shoco), Int,
                       (Ptr{Cchar}, Csize_t, Ptr{Cchar}, Csize_t),
                       s, sizeof(s), decompressed, 2 * sizeof(s))

        nbytes > 0 || throw(ErrorException("Decompression failed for input $s"))

        decompressed = decompressed[1:nbytes]
        decompressed[end] == Cchar(0) || push!(decompressed, Cchar(0))

        return unsafe_string(pointer(decompressed))
    end
end


text:
include(joinpath("Atom", "Atom.jl"))
include(joinpath("Molecule", "Molecule.jl"))


text:
export BIFMHelper, functional_dependencies

@doc raw"""
The BIFMHelper node is a node required to perform efficient message passing inconjuction with
the BIFM node. It is required to switch from the backward pass with messages to the forward pass
with marginals.

```julia
out ~ BIFMHelper(in)
```

Interfaces:
1. out - output of the BIFMHelper node, should be connected to the state space model.
2. in - input of the BIFMHelper node, should be connected to the prior for the latent state.

*Note: When performing inference, first subscribe to the marginals (in the order: z, out, in) and then to the free energy score function.*

"""
struct BIFMHelper <: AbstractFactorNode end

@node BIFMHelper Stochastic [ out, in ]

# specify custom functional dependencies for BIFMHelper node
function functional_dependencies(dependencies, factornode::FactorNode{ <: Type{BIFMHelper} }, iindex::Int)
    cindex             = clusterindex(factornode, iindex)

    nodeinterfaces     = interfaces(factornode)
    nodeclusters       = factorisation(factornode)
    nodelocalmarginals = localmarginals(factornode)

    varcluster = @inbounds nodeclusters[ cindex ]

    # output
    if iindex === 2
        mdependencies = (nodeinterfaces[1], )
        return tuple(mdependencies...), ()
    elseif iindex === 1
        qdependencies = TupleTools.deleteat(nodelocalmarginals, cindex)
        return (), tuple(qdependencies...)
    end

end

@average_energy BIFMHelper (q_out::Any, q_in::Any) = begin
    return entropy(q_in)
end


text:
@testset "operators" begin

pbar = Progress(4, 0.2, "operators")

# TODO this should be in a file test/operators.jl
@testset "map_edgevals! undirected" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        g = SimpleValueGraph(CycleGraph(V(5)), E_VAL)
        u = rand(vertices(g))
        w1 = rand_sample(E_VAL)
        w2 = rand_sample(E_VAL)
        map_edgevals!(g) do s, d, value
            return (u ∈ (s, d)) ? w1 : w2
        end
        @test all(edges(g)) do e
            return val(e) == (u ∈ (src(e), dst(e)) ? w1 : w2)
        end
    end
end

next!(pbar)

# TODO this should be in a file test/operators.jl
@testset "map_edgevals! directed" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        g = SimpleValueDiGraph(CycleDiGraph(V(5)), E_VAL)
        u = rand(vertices(g))
        w1 = rand_sample(E_VAL)
        w2 = rand_sample(E_VAL)
        map_edgevals!(g) do s, d, value
            return (s == u) ? w1 : w2
        end
        @test all(edges(g)) do e
            return val(e) == (src(e) == u ? w1 : w2)
        end
    end
end

next!(pbar)

# TODO this should be in a file test/operators.jl
@testset "map_edgevals! with key undirected" begin
    for (V, E_VAL) in product(test_vertex_types, filter(T -> T <: TupleOrNamedTuple, test_edgeval_types))
        g = SimpleValueGraph(CycleGraph(V(5)), E_VAL)
        u = rand(vertices(g))
        key = rand(1:length(E_VAL.parameters))
        w1 = rand_sample(E_VAL)[key]
        w2 = rand_sample(E_VAL)[key]
        map_edgevals!(g, key) do s, d, value
            return (u ∈ (s, d)) ? w1 : w2
        end
        @test all(edges(g)) do e
            return val(e)[key] == (u ∈ (src(e), dst(e)) ? w1 : w2)
        end
    end
end

next!(pbar)

# TODO this should be in a file test/operators.jl
@testset "map_edgevals! with key directed" begin
    for (V, E_VAL) in product(test_vertex_types, filter(T -> T <: TupleOrNamedTuple, test_edgeval_types))
        g = SimpleValueDiGraph(CycleDiGraph(V(5)), E_VAL)
        u = rand(vertices(g))
        key = rand(1:length(E_VAL.parameters))
        w1 = rand_sample(E_VAL)[key]
        w2 = rand_sample(E_VAL)[key]
        map_edgevals!(g, key) do s, d, value
            return (s == u) ? w1 : w2
        end
        @test all(edges(g)) do e
            return val(e)[key] == (u == src(e) ? w1 : w2)
        end
    end
end

next!(pbar)

end # testset


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "COINBLASBuilder"
version = v"1.4.6"

# Collection of sources required to build COINBLASBuilder
sources = [
    "https://github.com/coin-or-tools/ThirdParty-Blas/archive/releases/1.4.6.tar.gz" =>
    "f9601efb98f04fdba220d49d5bda98d2a5a5e2ed7564df339bc7149b0c303f0c",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ThirdParty-Blas-releases-1.4.6/
update_configure_scripts
./get.Blas
mkdir build
cd build/
## STATIC BUILD START
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
## STATIC BUILD END

## DYNAMIC BUILD START
#../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
## DYNAMIC BUILD START

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcoinblas", Symbol(""))
]

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

]

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


text:
##
using Pkg
Pkg.activate(".")
Pkg.instantiate()
include("src/EasyStocks.jl")

## Setup
using CairoMakie # For plotting
using Colors # To make it easier to select colors

mp = ModPar(q=0.0,ψ=0.0,xstar=2.0,η=0.)
np = NumPar(mp,nx=400,nsav=120,nq=1001,nα=101,nh=1,ygrd=[0.0,0.0],xmax=4.0,p1cons=false)
np = NumPar(mp,nx=150,nsav=120,nq=501,nα=101,nh=1,ygrd=[0.0,0.0],xmax=5.0,p1cons=false)
np.hgrd.=0.0
MS1 = ModelSolution(mp,np)
@time SolveModel!(MS1)
ψ = 0.015
mp = ModPar(q=mp.q,ψ=ψ,xstar=mp.xstar,η =mp.η)
MS2 = ModelSolution(mp,np)
@time SolveModel!(MS2)
np.πrs[1] = 0.002
np.πrs[:] = np.πrs/sum(np.πrs)

mp.ψ=0.0
mp.q=0.01
MS1c = ModelSolution(mp,np)
@time SolveModel!(MS1c)
mp.ψ=ψ
MS2c = ModelSolution(mp,np)
@time SolveModel!(MS2c)
mp.q=0.0

mp.ψ=0.0
MS1d = ModelSolution(mp,np)
@time SolveModel!(MS1d)
mp.ψ=ψ
MS2d = ModelSolution(mp,np)
@time SolveModel!(MS2d)
##
# Plot settings
using Colors # To make it easier to select colors
blue = colorant"#2B6EB2"
fig = Figure(resolution = (700, 300))
ax1  = Axis(fig[1,2],xlabel = L"Wealth $x$", ylabel = L"\alpha")
lines!(ax1,MS1.np.xgrd, MS1.α[:,1,1], linestyle=:dash,color=blue,label="Benchmark")
ixstar = findfirst(x->x==1,(np.xgrd .-mp.xstar).>0)
lines!(ax1,MS2.np.xgrd[1:ixstar-1], MS2.α[1:ixstar-1,1,1], color=:orange,label="Housing")
lines!(ax1,MS2.np.xgrd[ixstar:end], MS2.α[ixstar:end,1,1], color=:orange)
ylims!(ax1,(-0.05,1.05))
xlims!(ax1,(1,3))
fig

ax2  = Axis(fig[1,1],xlabel = L"Wealth $x$", ylabel = L"Utility $u(x')$")
lines!(ax2,MS1.np.xgrd, MS1.V[:,1,2], linestyle=:dash,color=blue)
ixstar = findfirst(x->x==1,(np.xgrd .-mp.xstar).>0)
lines!(ax2,MS2.np.xgrd[1:ixstar-1], MS2.V[1:ixstar-1,1,2], color=:orange)
lines!(ax2,MS2.np.xgrd[ixstar:end], MS2.V[ixstar:end,1,2], color=:orange)
ylims!(ax2,(-1.035,-0.165))
xlims!(ax2,(1,3))

leg = Legend(fig[2, 1:2], ax1,orientation = :horizontal,framevisible = false)
Label(fig[1, 1, Top()], L"Terminal Value $U(x')$", valign = :bottom,padding=[0 0 5 0])
Label(fig[1, 2, Top()], L"Portfolio Weight $\alpha$", valign = :bottom,padding=[0 0 5 0])
rowsize!(fig.layout, 2, Fixed(0.1))
hidespines!(ax2,:t,:r)
hidespines!(ax1,:t,:r)
save("figs/value_weight.pdf", fig)
fig

## Plots with disaster
fig = Figure(resolution = (700, 300))
ax1  = Axis(fig[1,1],xlabel = L"Wealth $x$", ylabel = L"\alpha")
lines!(ax1,MS1d.np.xgrd, MS1d.α[:,1,1], linestyle=:dash,color=blue,label="Benchmark")
lines!(ax1,MS2d.np.xgrd[1:ixstar-1], MS2d.α[1:ixstar-1,1,1], color=:orange,label="Housing")
lines!(ax1,MS2d.np.xgrd[ixstar:end], MS2d.α[ixstar:end,1,1], color=:orange)
ylims!(ax1,(-0.05,1.05))
xlims!(ax1,(1,3))

ax2  = Axis(fig[1,2],xlabel = L"Wealth $x$", ylabel = L"\alpha")
discont=sort!([1;ixstar;findfirst(x->x>0,MS2c.α[:,1,1]);findlast(x->x==0,MS2c.α[:,1,1])+1;length(MS2c.α[:,1,1])])
lines!(ax2,MS1d.np.xgrd[1:discont[2]-1], MS1c.α[1:discont[2]-1,1,1], linestyle=:dash,color=blue,label="Benchmark")
lines!(ax2,MS1d.np.xgrd[discont[2]:end], MS1c.α[discont[2]:end,1,1], linestyle=:dash,color=blue)
for i in 1:length(discont)-1
    lines!(ax2,MS2d.np.xgrd[discont[i]:discont[i+1]-1], MS2c.α[discont[i]:discont[i+1]-1], color=:orange,label="Housing")
end
ylims!(ax2,(-0.05,1.05))
xlims!(ax2,(1,3.))
leg = Legend(fig[2, 1:2], ax1,orientation = :horizontal,framevisible = false)
Label(fig[1, 1, Top()], L"Disaster Risk $p_{tail}>0$", valign = :bottom,padding=[0 0 5 0])
Label(fig[1, 2, Top()], L"Participation cost $q>0$", valign = :bottom,padding=[0 0 5 0])
rowsize!(fig.layout, 2, Fixed(0.1))
hidespines!(ax2,:t,:r)
hidespines!(ax1,:t,:r)
save("figs/disaster_cost.pdf", fig)
fig


text:
using Test
using CbiomesProcessing

println("Testing...")

max_issues_repo_name: darwinproject/CBIOMES-Processing-legacy

max_forks_repo_name: darwinproject/CBIOMES-Processing-legacy


text:
using TightlyBound
using Test
using Suppressor

function test1()

    @testset "testing ewald" begin

        @suppress begin
            cscl = makecrys([1.0 0 0; 0 1.0 0; 0 0 1.0]*2/sqrt(3), [0 0 0; 0.5 0.5 0.5], ["Na", "Cl"], units="Bohr");

            gamma_ij_tot_k1 = TightlyBound.Ewald.electrostatics_getgamma(cscl, kappa = 1.0, noU=true)
            gamma_ij_tot_k2 = TightlyBound.Ewald.electrostatics_getgamma(cscl, kappa = 2.0, noU=true)

            rs = makecrys([1.0 1.0 0; 1.0 0.0 1.0; 0 1.0 1.0], [0 0 0; 0.5 0.5 0.5], ["Na", "Cl"], units="Bohr");

            kappa = 3.0
            gamma_ij_tot_rs = TightlyBound.Ewald.electrostatics_getgamma(rs, kappa = kappa, noU=true)

            @test  isapprox(gamma_ij_tot_k1[1,2] - gamma_ij_tot_k1[1,1], gamma_ij_tot_k2[1,2] - gamma_ij_tot_k2[1,1], rtol=1e-5) #do different kappa values agree?
            @test  isapprox(sum( (gamma_ij_tot_k1[1,2] - gamma_ij_tot_k1[1,1])), 2*1.76268, rtol = 1e-4) #madelung constant cscl
            @test  isapprox(sum( (gamma_ij_tot_rs[1,2] - gamma_ij_tot_rs[1,1])), 2*1.74757, rtol = 1e-4) #madelung constant nacl

        end

#        rs_real = CrystalMod.makecrys([10.0 10.0 0; 10.0 0.0 10.0; 0 10.0 10.0], [0 0 0; 0.5 0.5 0.5], ["Na", "Cl"], units="Bohr");
#        gamma_ij_tot_rs = Ewald.electrostatics_getgamma(rs_real)
#        rs_real = CrystalMod.makecrys([3.0 3.0 0; 3.0 0.0 3.0; 0 3.0 3.0], [0 0 0; 0.5 0.5 0.5], ["Na", "Cl"], units="Bohr");
#        gamma_ij_tot_rs = Ewald.electrostatics_getgamma(rs_real)

#        for kappa = exp10.(range(log10(0.05), stop=log10(10), length=30))
#            println("KAPPA $kappa--------------------------------")
#            gamma_ij_tot_rs = Ewald.electrostatics_getgamma(rs_real, kappa = kappa)
#        end

    end
end

function test2()

    @testset "testing ewald units - dipole" begin

        @suppress begin

            c = makecrys([1.0 0 0; 0 1.0 0; 0 0 1.0]*100, [0 0 0; 0.0 0.0 0.01], ["Na", "Cl"], units="Bohr");

            gamma = TightlyBound.Ewald.electrostatics_getgamma(c, noU=true)

            energy, pot = TightlyBound.TB.ewald_energy(c, gamma, [1.0, -1.0])

#        println("energy $energy")

            @test  isapprox(energy, -2.0, rtol=1e-5) #energy in rydberg e2=2; r=1;q=1; e2*q*-q/r  == -2, other atoms far away
        end

    end
end

test1()
test2()


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

JLLWrappers.@generate_main_file_header("Libmount")
JLLWrappers.@generate_main_file("Libmount", UUID("4b2f31a3-9ecc-558c-b454-b3730dcb73e9"))
end  # module Libmount_jll


text:
module ResponseManneville

import IntervalArithmetic: Interval
export Interval

abstract type NormKind end

struct L1 <: NormKind end
struct Linf <: NormKind end
struct Lipschitz <: NormKind end
struct TotalVariation <: NormKind end
struct ℓ1 <: NormKind end
struct ℓinf <: NormKind end

include("DynamicDefinition.jl")
include("BasisDefinition.jl")
include("Contractors.jl")
include("GenericAssembler.jl")
include("GenericEstimate.jl")
include("PwDynamicDefinition.jl")
include("Norms.jl")
include("NormsOfPowers.jl")
include("UlamBasis.jl")
include("C2Basis.jl")
include("InducedLSV.jl")
include("ContractionC1.jl")

using .C2BasisDefinition, .InducedLSVMapDefinition
export ApproxInducedLSV, C2Basis, C1, Ulam, assemble, DiscretizedOperator

end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
@testset "composite" begin

PetscInitialize()

# Create two matrices
A = create_matrix(2, 2; auto_setup = true)
B = create_matrix(2, 2; auto_setup = true)

# Fill

# A
# 1 -
# 3 4
A[1,1] = 1
A[2,1] = 3
A[2,2] = 4

# B
# 4 3
# - 1
B[1,1] = 4
B[1,2] = 3
B[2,2] = 1

# Assemble matrices
assemble!.((A,B))

# Create composite mat from A and B
# C
#  5 3
#  3 5
C = create_composite_add([A, B])

# Create vectors to check C (see below)
x1 = create_vector(2; auto_setup = true)
x2 = create_vector(2; auto_setup = true)
y  = create_vector(2; auto_setup = true)
x1[1] = 1.; x1[2] = 0.
x2[1] = 0.; x2[2] = 1.
assemble!.((x1,x2,y))

# Check result by multiplying with test vectors
MatMult(C, x1, y)
@test vec2array(y) == [5., 3.]

MatMult(C, x2, y)
@test vec2array(y) == [3., 5.]

# Free memory
destroy!.((A, B, C, x1, x2, y))
PetscFinalize()
end

max_stars_repo_name: bmxam/JuliaPetsc

max_issues_repo_name: bmxam/JuliaPetsc

max_forks_repo_name: bmxam/JuliaPetsc


text:
# module NonLazyScalar
#     import Base.Broadcast: broadcasted
#     import ExBroadcast.AbstractScalar
#     const Scalar = AbstractScalar{<:Number}

#     ## broadcasted expand
#     for op in (:+, :*, :&, :|, :xor, :min, :max, :kron)
#         @eval begin
#             @inline broadcasted(::typeof($op),x,y,z,args...) =
#                 broadcasted($op, broadcasted($op,x,y), z, args...)
#         end
#     end

#     for op in (:muladd, :fma,)
#         @eval broadcasted(f::typeof($op), x::Scalar, y::Scalar, z::Scalar) = f(only(x),only(y),only(z))
#     end

#     for op in (
#         :+, :-, :*, :/, :\, :÷, :%, :^, :fld, :cld, :mod,
#         :&, :|, :⊻, :>>>, :>>, :<<,
#         :atan, :hypot,
#         :min, :max,
#         :binomial, :gcd, :lcm,
#         )
#         @eval broadcasted(f::typeof($op), x::Scalar, y::Scalar) = f(only(x),only(y))
#     end

#     for op in (
#         :+, :-, :!, :inv, :~,
#         :sin  ,:cos  ,:tan  ,:sind ,:cosd ,:tand ,:sinpi,:cospi,
#         :asin ,:acos ,:atan ,:asind,:acosd,:atand,
#         :sec  ,:csc  ,:cot  ,:secd ,:cscd ,:cotd ,
#         :asec ,:acsc ,:acot ,:asecd,:acscd,:acotd,
#         :sinh ,:cosh ,:tanh ,:sech ,:csch ,:coth ,
#         :asinh,:acosh,:atanh,:asech,:acsch,:acoth,
#         :sinc,:cosc,
#         :deg2rad,:rad2deg,
#         :log,:log2,:log10,:log1p,
#         :exp,:exp2,:exp10,:ldexp,:expm1,
#         :round,:ceil,:floor,:trunc,
#         :abs,:abs2,:sign,:sqrt,:cbrt,:real,:imag,:conj,:angle,:cis,
#         :factorial,:big,
#         :identity
#         )
#         @eval broadcasted(f::typeof($op), x::Scalar) = f(only(x))
#     end
#     broadcasted(::Type{T}, x::Scalar) where {T} = T(only(x))
# end
# using .NonLazyScalar


text:
using StorageMirrorServer: get_hash, RegistryMeta, query_latest_hash, url_exists
using StorageMirrorServer: download_and_verify
using StorageMirrorServer: timeout_call, TimeoutException

@testset "normalize_upstream" begin
    upstreams = [
        "https://pkg.julialang.org",
        "https://pkg.julialang.org/",
        "pkg.julialang.org",
        "pkg.julialang.org/",
    ]
    for server in upstreams
        @test "https://pkg.julialang.org" == StorageMirrorServer.normalize_upstream(server)
    end
end

registry = RegistryMeta(
    "General",
    "23338594-aafe-5451-b93e-139f81909106",
    "https://github.com/JuliaRegistries/General"
)

upstreams = [
    # "https://mirrors.sjtug.sjtu.edu.cn/julia",
    # "https://mirrors.bfsu.edu.cn/julia"
    "https://pkg.julialang.org",
    "https://us-east.storage.juliahub.com",
]

@testset "get_hash" begin
    treehash = "e1e1667fcb629f2c9ec4ac43464c2ef6177ed696"
    @test treehash == get_hash("/registry/23338594-aafe-5451-b93e-139f81909106/$treehash\n", registry.uuid)
    @test treehash == get_hash(IOBuffer("/registry/23338594-aafe-5451-b93e-139f81909106/$treehash\n"), registry.uuid)

    # uuid not found
    @test nothing === get_hash("/registry/00000000-1111-2222-3333-444444444444/$treehash\n", registry.uuid)
    # invalid lines
    @test nothing === get_hash("hello world", registry.uuid)
    # whitespaces are ignored
    @test treehash == get_hash("hello world\n /registry/23338594-aafe-5451-b93e-139f81909106/$treehash \n", registry.uuid)
end

@testset "url_exists" begin
    # @test url_exists("https://pkg.julialang.org/registries"; timeout=0)
    @test !url_exists("https://pkg.julialang.org/registries_1234")

    @test_throws ArgumentError url_exists("abc")
    err_msg = @capture_err url_exists("https://abc")
    @test occursin("failed to send HEAD request", err_msg)
end

@testset "query_latest_hash" begin
    @testset "single upstream" begin
        @test !isnothing(match(r"[0-9a-f]{40}", query_latest_hash(registry, upstreams[1])))
    end
end

@testset "download_and_verify" begin
    tmp_testdir = mktempdir()

    resource_list = [
        "/artifact/ff8ad169326afd41d46f507a960c1717f4ed1a47",
        "/artifact/f8eb2f3aa430c1ea80c46779ba89580372f0f0db",
        "/package/a603d957-0e48-4f86-8fbd-0b7bc66df689/e4581e3fadda3824e0df04396c85258a2107035d",
        "/package/22bb73d7-edb2-5785-ba1e-7d60d6824784/41731998eb760f9c5e4acc911c9fb33c9643365f",
    ]
    for resource in resource_list
        tarball = joinpath(tmp_testdir, resource[2:end])
        @test download_and_verify(upstreams, resource, tarball)
        @test isfile(tarball)
    end

    server = "https://mirrors.bfsu.edu.cn/julia"
    # invalid hash
    resource = "/artifact/0000"
    tarball = joinpath(tmp_testdir, resource[2:end])
    @test @suppress_err !download_and_verify(server, resource, tarball)
    err_msg = @capture_err download_and_verify(server, resource, tarball)
    @test occursin("bad resource: valid hash not found", err_msg)

    # 404 test
    resource = "/artifact/0000000000000000000000000000000000000000"
    tarball = joinpath(tmp_testdir, resource[2:end])
    @test @suppress_err !download_and_verify(server, resource, tarball)
    err_msg = @capture_err download_and_verify(server, resource, tarball)
    @test occursin("failed to fetch resource", err_msg)

    # timeout test
    # server = "https://us-east.storage.juliahub.com"
    # resource = resource_list[1]
    # tarball = joinpath(tmp_testdir, resource[2:end])
    # config = Dict{Symbol, Any}(:timeout => 0.001)
    # @test @suppress_err !download_and_verify(server, resource, tarball; http_parameters=config)
    # err_msg = @capture_err download_and_verify(server, resource, tarball; http_parameters=config)
    # @test occursin("failed to fetch resource", err_msg)
    # @test occursin("TimeoutException(0.001)", err_msg)
    # @test !isfile(tarball)

    # TODO: test hash mismatch
end


text:
# Different ways to simulate molecules

export
    SteepestDescentMinimizer,
    simulate!,
    VelocityVerlet,
    Verlet,
    StormerVerlet,
    Langevin

"""
    SteepestDescentMinimizer(; <keyword arguments>)

Steepest descent energy minimization.

# Arguments
- `step_size::D=0.01u"nm"`: the initial maximum displacement.
- `max_steps::Int=1000`: the maximum number of steps.
- `tol::F=1000.0u"kJ * mol^-1 * nm^-1"`: the maximum force below which to
    finish minimization.
- `run_loggers::Bool=false`: whether to run the loggers during minimization.
- `log_stream::L=devnull`: stream to print minimization progress to.
"""
struct SteepestDescentMinimizer{D, F, L}
    step_size::D
    max_steps::Int
    tol::F
    run_loggers::Bool
    log_stream::L
end

function SteepestDescentMinimizer(;
                                    step_size=0.01u"nm",
                                    max_steps=1_000,
                                    tol=1000.0u"kJ * mol^-1 * nm^-1",
                                    run_loggers=false,
                                    log_stream=devnull)
    return SteepestDescentMinimizer(step_size, max_steps, tol,
                                    run_loggers, log_stream)
end

"""
    simulate!(system, simulator, n_steps; parallel=true)
    simulate!(system, simulator; parallel=true)

Run a simulation on a system according to the rules of the given simulator.
Custom simulators should implement this function.
"""
function simulate!(sys,
                    sim::SteepestDescentMinimizer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    sim.run_loggers && run_loggers!(sys, neighbors, 0; parallel=parallel)
    E = potential_energy(sys, neighbors)
    println(sim.log_stream, "Step 0 - potential energy ",
            E, " - max force N/A - N/A")
    hn = sim.step_size

    for step_n in 1:sim.max_steps
        F = forces(sys, neighbors; parallel=parallel)
        max_force = maximum(norm.(F))

        coords_copy = sys.coords
        sys.coords += hn * F ./ max_force
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        neighbors_copy = neighbors
        neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                    parallel=parallel)
        E_trial = potential_energy(sys, neighbors)
        if E_trial < E
            hn = 6 * hn / 5
            E = E_trial
            println(sim.log_stream, "Step ", step_n, " - potential energy ",
                    E_trial, " - max force ", max_force, " - accepted")
        else
            sys.coords = coords_copy
            neighbors = neighbors_copy
            hn = hn / 5
            println(sim.log_stream, "Step ", step_n, " - potential energy ",
                    E_trial, " - max force ", max_force, " - rejected")
        end

        sim.run_loggers && run_loggers!(sys, neighbors, step_n;
                                        parallel=parallel)

        if max_force < sim.tol
            break
        end
    end
    return sys
end

# Forces are often expressed per mol but this dimension needs removing for use in the integrator
function remove_molar(x)
    fx = first(x)
    if dimension(fx) == u"𝐋 * 𝐍^-1 * 𝐓^-2"
        T = typeof(ustrip(fx))
        return x / T(Unitful.Na)
    else
        return x
    end
end

"""
    VelocityVerlet(; <keyword arguments>)

The velocity Verlet integrator.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct VelocityVerlet{T, C}
    dt::T
    coupling::C
    remove_CM_motion::Bool
end

function VelocityVerlet(; dt, coupling=NoCoupling(), remove_CM_motion=true)
    return VelocityVerlet(dt, coupling, remove_CM_motion)
end

function simulate!(sys,
                    sim::VelocityVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    run_loggers!(sys, neighbors, 0; parallel=parallel)
    accels_t = accelerations(sys, neighbors; parallel=parallel)
    accels_t_dt = zero(accels_t)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        sys.coords += sys.velocities .* sim.dt .+ (remove_molar.(accels_t) .* sim.dt ^ 2) ./ 2
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        accels_t_dt = accelerations(sys, neighbors; parallel=parallel)

        sys.velocities += remove_molar.(accels_t .+ accels_t_dt) .* sim.dt / 2

        sim.remove_CM_motion && remove_CM_motion!(sys)
        apply_coupling!(sys, sim, sim.coupling)

        run_loggers!(sys, neighbors, step_n; parallel=parallel)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
            accels_t = accels_t_dt
        end
    end
    return sys
end

"""
    Verlet(; <keyword arguments>)

The leapfrog Verlet integrator.
This is a leapfrog integrator, so the velocities are offset by half a time step
behind the positions.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct Verlet{T, C}
    dt::T
    coupling::C
    remove_CM_motion::Bool
end

function Verlet(; dt, coupling=NoCoupling(), remove_CM_motion=true)
    return Verlet(dt, coupling, remove_CM_motion)
end

function simulate!(sys,
                    sim::Verlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    run_loggers!(sys, neighbors, 0; parallel=parallel)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        accels_t = accelerations(sys, neighbors; parallel=parallel)

        sys.velocities += remove_molar.(accels_t) .* sim.dt

        sys.coords += sys.velocities .* sim.dt
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        sim.remove_CM_motion && remove_CM_motion!(sys)
        apply_coupling!(sys, sim, sim.coupling)

        run_loggers!(sys, neighbors, step_n; parallel=parallel)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
        end
    end
    return sys
end

"""
    StormerVerlet(; <keyword arguments>)

The Störmer-Verlet integrator.
Does not currently work with coupling methods that alter the velocity.

# Arguments
- `dt::T`: the time step of the simulation.
- `coupling::C=NoCoupling()`: the coupling which applies during the simulation.
"""
struct StormerVerlet{T, C}
    dt::T
    coupling::C
end

StormerVerlet(; dt, coupling=NoCoupling()) = StormerVerlet(dt, coupling)

function simulate!(sys,
                    sim::StormerVerlet,
                    n_steps::Integer;
                    parallel::Bool=true)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    run_loggers!(sys, neighbors, 0; parallel=parallel)
    coords_last = sys.coords

    for step_n in 1:n_steps
        accels_t = accelerations(sys, neighbors; parallel=parallel)

        coords_copy = sys.coords
        if step_n == 1
            # Use the velocities at the first step since there is only one set of coordinates
            sys.coords += sys.velocities .* sim.dt .+ (remove_molar.(accels_t) .* sim.dt ^ 2) ./ 2
        else
            sys.coords += vector.(coords_last, sys.coords, (sys.box_size,)) .+ remove_molar.(accels_t) .* sim.dt ^ 2
        end
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))

        # This is accurate to O(dt)
        sys.velocities = vector.(coords_copy, sys.coords, (sys.box_size,)) ./ sim.dt

        apply_coupling!(sys, sim, sim.coupling)

        run_loggers!(sys, neighbors, step_n; parallel=parallel)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
            coords_last = coords_copy
        end
    end
    return sys
end

"""
    Langevin(; <keyword arguments>)

The Langevin integrator, based on the Langevin Middle Integrator in OpenMM.
This is a leapfrog integrator, so the velocities are offset by half a time step
behind the positions.

# Arguments
- `dt::T`: the time step of the simulation.
- `temperature::K`: the temperature of the simulation.
- `friction::F`: the friction coefficient of the simulation.
- `remove_CM_motion::Bool=true`: whether to remove the centre of mass motion
    every time step.
"""
struct Langevin{S, K, F, T}
    dt::S
    temperature::K
    friction::F
    remove_CM_motion::Bool
    vel_scale::T
    noise_scale::T
end

function Langevin(; dt, temperature, friction, remove_CM_motion=true)
    vel_scale = exp(-dt * friction)
    noise_scale = sqrt(1 - vel_scale^2)
    return Langevin(dt, temperature, friction, remove_CM_motion, vel_scale, noise_scale)
end

function simulate!(sys,
                    sim::Langevin,
                    n_steps::Integer;
                    parallel::Bool=true,
                    rng=Random.GLOBAL_RNG)
    neighbors = find_neighbors(sys, sys.neighbor_finder; parallel=parallel)
    run_loggers!(sys, neighbors, 0; parallel=parallel)
    sim.remove_CM_motion && remove_CM_motion!(sys)

    for step_n in 1:n_steps
        accels_t = accelerations(sys, neighbors; parallel=parallel)

        sys.velocities += remove_molar.(accels_t) .* sim.dt

        sys.coords += sys.velocities .* sim.dt / 2
        noise = random_velocities(sys, sim.temperature; rng=rng)
        sys.velocities = sys.velocities .* sim.vel_scale .+ noise .* sim.noise_scale

        sys.coords += sys.velocities .* sim.dt / 2
        sys.coords = wrap_coords_vec.(sys.coords, (sys.box_size,))
        sim.remove_CM_motion && remove_CM_motion!(sys)

        run_loggers!(sys, neighbors, step_n; parallel=parallel)

        if step_n != n_steps
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                        parallel=parallel)
        end
    end
    return sys
end


text:
"""
    @setup_procs

Lorem ipsum dolor sit amet.
"""
macro setup_procs()
  load_input("input.jl", true)

  [
    f([randn(3,3) for i=1:2], *)
    for f in (prefix_serial!, prefix_parallel!)
  ]

  quote
    r = [
      @spawnat p randn(n,n)
      for p in procs()[1:nprocs()]
    ]

    s = fetch(r) #These are in-memory matrices

    r, s
  end
end


text:
include("information-plane.jl")

epochs = 10000
number_nets = 50
step_len = 1
#=
x_train, y_train, x_test, y_test = get_datasetDemo("../datasets/Demo", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(datasetDemo_nfeatures, datasetDemo_labels, [2], number_nets, data, epochs, step_len, "IP-dtDemo-MLP.png")

x_train, y_train, x_test, y_test = get_datasetDemoOverlap("../datasets/Demo-overlap", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(datasetOverlap_nfeatures, datasetOverlap_labels, [2], number_nets, data, epochs, step_len, "IP-dtOverlap-MLP.png")

x_train, y_train, x_test, y_test = get_dataset1("../datasets/Drug-Classification", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(dataset1_nfeatures, dataset1_labels, [5], number_nets, data, epochs, step_len, "IP-dt1-MLP.png")

x_train, y_train, x_test, y_test = get_dataset2("../datasets/Income-Classification", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(dataset2_nfeatures, dataset2_labels, [26,13,6], number_nets, data, epochs, step_len, "IP-dt2-MLP.png")
=#
x_train, y_train, x_test, y_test = get_dataset3("../datasets/Wine-Classification", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(dataset3_nfeatures, dataset3_labels, [22,10,15,11], number_nets, data, epochs, step_len, "IP-dt3-MLP.png")

x_train, y_train, x_test, y_test = get_dataset4("../datasets/Drug-Consumism-Classification", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(dataset4_nfeatures, dataset4_labels, [24,6,12,7], number_nets, data, epochs, step_len, "IP-dt4-MLP.png")

x_train, y_train, x_test, y_test = get_dataset5("../datasets/Heart-Disease-Classification", normalise=true)
data = zip(x_train, y_train)
plot_IP_average(dataset5_nfeatures, dataset5_labels,  [26,13,7,5], number_nets, data, epochs, step_len, "IP-dt5-MLP.png")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
    decasteljau(control, t)

Evaluate de Casteljau's algorithm.
"""
function decasteljau(control::AbstractArray{T,2}, t::T) where {T<:AbstractFloat}
    n = size(control, 2) - 1
    β = zeros(size(control, 1), n+1, n+1)
    @inbounds β[:,:,1] = control
    for j = 2:n+1, i = 1:n-j+2
        @inbounds β[:,i,j] = β[:,i,j-1]*(1-t) + β[:,i+1,j-1]*t
    end
    β
end

"""
    subdivide(control, t)

Subdivide a Bezier curve using de Casteljau's algorithm.
"""
function subdivide(control::AbstractArray{T,2}, t::T) where {T<:AbstractFloat}
    β = decasteljau(control, t)
    β₁ = β[:,1,:]
    β₂ = zeros(control)
    n = size(control, 2) - 1
    for i = 1:n+1
        @inbounds β₂[:,i] = β[:,i,n-i+2]
    end
    β₁, β₂
end


text:
# Export Segments API
export Segments, handle

# Export Segment API
export Segment, deref, segment_name, segment_offset, segment_file_size,
       segment_memory_size, segment_address

export SegmentRef, segment_number

# Import Base stuff for extension
import Base: getindex, length, eltype, findfirst, findall, iterate, keys, lastindex

"""
    Segments

An abstraction over the concept of a collection of `Segment` types within an
object file.  One can think of the `Segments` object containing the table of
segment headers within the object file, whereas the `Segment`/`SegmentRef`
objects contain the actual segment data itself.  The list of available API
operations is given below, with methods that subclasses must implement marked
in emphasis:

### Creation
  - *Segments()*

### Iteration
  - *getindex()*
  - *lastindex()*
  - length()
  - iterate()
  - keys()
  - eltype()

### Search
  - findall()
  - findfirst()

### Misc.
  - *handle()*
"""
abstract type Segments{H<:ObjectHandle} end

# Fairly simple iteration interface specification
@mustimplement lastindex(segs::Segments)
length(segs::Segments) = lastindex(segs)
keys(segs::Segments) = 1:length(segs)
iterate(segs::Segments, idx=1) = idx > length(segs) ? nothing : (segs[idx], idx+1)
eltype(::Type{S}) where {S <: Segments} = SegmentRef

function getindex(segs::Segments, idx)
    # Punt off to `getindex_ref`
    oh = handle(segs)
    return getindex_ref(
        segs,
        segment_header_offset(oh),
        segment_header_size(oh),
        segment_header_type(oh),
        SegmentRef,
        idx
    )
end

"""
    Segment

An abstraction over the concept of a `Segment` within an object file.  A
`Segment` is a portion of an object file that is given instruction on its
layout in virtual memory; this is in contrast to a `Section`, which delineates
different portions of an object file on disk.  ELF files have the strictest
separation here, with a single executable file containing multiple `Segment`
and `Section` objects, with `Section`s being assigned to one or more `Segment`s
for virtual memory placement.  Mach-O files typically have two `Segment`s, one
called `__TEXT`, one called `__DATA`.  COFF files do not have `Segment`.

Just like with `Section` objects, many operations upon segments require global
operations (access to the string table, knowledge of position within the file,
etc...) which causes some operations to be defined only upon the `SegmentRef`
datatype.  As a user, the `SegmentRef` type should be the primary method of
interacting with segments, as a developer adding new object file formats,
some methods must support `Segment`s, others must support only `SegmentRef`s.
Note that any method that works on a `Segment` must also work with a
`SegmentRef`, see the `@derefmethod` macro for a convenient helper macro to
generate `SegmentRef` -> `Section` wrapper methods. The list of available API
operations is given below, with methods that subclasses must implement marked
in emphasis:

### Creation:
  - *read()*

### Utility:
  - deref()

### Format-specific properties:
  - *segment_name()*
  - *segment_offset()*
  - *segment_file_size()*
  - *segment_memory_size()*
  - *segment_address()*
"""
abstract type Segment{H <: ObjectHandle} end

deref(seg::Segment) = seg
@mustimplement read(oh::ObjectHandle, ::Type{Segment})
@mustimplement segment_name(seg::Segment)
@mustimplement segment_offset(seg::Segment)
@mustimplement segment_file_size(seg::Segment)
@mustimplement segment_memory_size(seg::Segment)
@mustimplement segment_address(seg::Segment)

"""
    SegmentRef

Provides a reference to a `Segment`, along with a reference to the
`ObjectHandle` this `Segment` comes from.  This should be the primary method by
which users interact with segments inside object files.  The list of available
API operations is given below, with methods that subclasses must implement
marked in emphasis.  Note that this overlaps heavily with the `Segment` object
API, this is by design as many of the methods are simply passthroughs to the
underlying `Segment` API calls for ease of use.

### Creation:
  - *SegmentRef()*

### Utility
  - *deref()*
  - *Segments()*
  - handle()

### Format-specific properties:
  - *segment_name()*
  - *segment_number()*
  - segment_offset()
  - segment_file_size()
  - segment_memory_size()
  - segment_address()
"""
abstract type SegmentRef{H<:ObjectHandle} end

@mustimplement SegmentRef(segs::Segments, seg::Segment, idx)

"""
    deref(seg::SegmentRef)

Dereference the given `SegmentRef` object to a `Segment`.
"""
@mustimplement deref(seg::SegmentRef)

"""
    Segments(seg::SegmentRef)

Return the `Segments` collection this `Segment` belongs to.
"""
@mustimplement Segments(seg::SegmentRef)

"""
    handle(seg::SegmentRef)

Return the `ObjectHandle` this `SegmentRef` belongs to.  This method is
`SegmentRef`-only.
"""
handle(seg::SegmentRef) = handle(Segments(seg))

"""
    segment_name(seg::SegmentRef)

The name of the given `Segment`, returned as a string.  This method often
performs some kind of lookup within the string table of the object to get the
full name of the segment.
"""
@mustimplement segment_name(seg::SegmentRef)

"""
    segment_number(seg::SegmentRef)

Return the index of the referred segment.
"""
@mustimplement segment_number(seg::SegmentRef)

@derefmethod segment_type(seg::SegmentRef)
@derefmethod segment_offset(seg::SegmentRef)
@derefmethod segment_file_size(seg::SegmentRef)
@derefmethod segment_memory_size(seg::SegmentRef)
@derefmethod segment_address(seg::SegmentRef)

## Printing
show(io::IO, segs::Segments{H}) where {H <: ObjectHandle} = show_collection(io, segs, H)


text:
using Pkg;
Pkg.activate("");
using CUTEst, Dates, NLPModels, SolverBenchmark

using ARCTR, JSOSolvers, NLPModelsIpopt

nmax = 10
problems = readlines("list_problems_$nmax.dat")
cutest_problems = (CUTEstModel(p) for p in problems)

max_time = 3600.0
max_ev = typemax(Int)
max_iter = typemax(Int)
atol = 1e-5
rtol = 1e-6

solvers = Dict(
    :ARCqKOp =>
        nlp -> ARCqKOp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
        ),
    :ST_TROp =>
        nlp -> ST_TROp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
        ),
)
stats = bmark_solvers(solvers, cutest_problems)

using JLD2
@save "$(today())_$(prod(String.(keys(solvers)) .* "_"))cutest_$(string(length(problems)))_$(nmax).jld2" stats


text:
# Grabbed from Sfxs.cs
# These are the constant fields in the class and they seem to be correct

module Songs
macro songnames(s::Symbol, expr::Expr)
    return Expr(:block, esc(:($s = $expr)), [esc(:($(Meta.parse(k)) = $v)) for (k, v) in eval(expr)]...)
end

@songnames(songs, Dict{String, String}(
    "music_levelselect" => "event:/music/menu/level_select",
    "music_credits" => "event:/music/menu/credits",
    "music_complete_area" => "event:/music/menu/complete_area",
    "music_complete_summit" => "event:/music/menu/complete_summit",
    "music_complete_bside" => "event:/music/menu/complete_bside",
    "music_prologue_intro_vignette" => "event:/game/00_prologue/intro_vignette",
    "music_prologue_beginning" => "event:/music/lvl0/intro",
    "music_prologue_collapse" => "event:/music/lvl0/bridge",
    "music_prologue_title_ping" => "event:/music/lvl0/title_ping",
    "music_city" => "event:/music/lvl1/main",
    "music_city_theo" => "event:/music/lvl1/theo",
    "music_oldsite_beginning" => "event:/music/lvl2/beginning",
    "music_oldsite_mirror" => "event:/music/lvl2/mirror",
    "music_oldsite_dreamblock_sting_pt1" => "event:/music/lvl2/dreamblock_sting_pt1",
    "music_oldsite_dreamblock_sting_pt2" => "event:/music/lvl2/dreamblock_sting_pt2",
    "music_oldsite_evil_maddy" => "event:/music/lvl2/evil_madeline",
    "music_oldsite_chase" => "event:/music/lvl2/chase",
    "music_oldsite_payphone_loop" => "event:/music/lvl2/phone_loop",
    "music_oldsite_payphone_end" => "event:/music/lvl2/phone_end",
    "music_oldsite_awake" => "event:/music/lvl2/awake",
    "music_resort_intro" => "event:/music/lvl3/intro",
    "music_resort_explore" => "event:/music/lvl3/explore",
    "music_resort_clean" => "event:/music/lvl3/clean",
    "music_resort_clean_extended" => "event:/music/lvl3/clean_extended",
    "music_resort_oshiro_theme" => "event:/music/lvl3/oshiro_theme",
    "music_resort_oshiro_chase" => "event:/music/lvl3/oshiro_chase",
    "music_cliffside_main" => "event:/music/lvl4/main",
    "music_cliffside_heavywinds" => "event:/music/lvl4/heavy_winds",
    "music_cliffside_panicattack" => "event:/music/lvl4/minigame",
    "music_temple_normal" => "event:/music/lvl5/normal",
    "music_temple_middle" => "event:/music/lvl5/middle_temple",
    "music_temple_mirror" => "event:/music/lvl5/mirror",
    "music_temple_mirrorcutscene" => "event:/music/lvl5/mirror_cutscene",
    "music_reflection_maddietheo" => "event:/music/lvl6/madeline_and_theo",
    "music_reflection_starjump" => "event:/music/lvl6/starjump",
    "music_reflection_fall" => "event:/music/lvl6/the_fall",
    "music_reflection_fight" => "event:/music/lvl6/badeline_fight",
    "music_reflection_fight_glitch" => "event:/music/lvl6/badeline_glitch",
    "music_reflection_fight_finish" => "event:/music/lvl6/badeline_acoustic",
    "music_reflection_main" => "event:/music/lvl6/main",
    "music_reflection_secretroom" => "event:/music/lvl6/secret_room",
    "music_summit_main" => "event:/music/lvl7/main",
    "music_summit_finalascent" => "event:/music/lvl7/final_ascent",
    "music_epilogue_main" => "event:/music/lvl8/main",
    "music_core_main" => "event:/music/lvl9/main",
    "music_pico8_title" => "event:/classic/pico8_mus_00",
    "music_pico8_area1" => "event:/classic/pico8_mus_01",
    "music_pico8_area2" => "event:/classic/pico8_mus_02",
    "music_pico8_area3" => "event:/classic/pico8_mus_03",
    "music_pico8_wind" => "event:/classic/sfx61",
    "music_pico8_end" => "event:/classic/sfx62",
    "music_pico8_boot" => "event:/classic/pico8_boot",
    "music_rmx_01_forsaken_city" => "event:/music/remix/01_forsaken_city",
    "music_rmx_02_old_site" => "event:/music/remix/02_old_site",
    "music_rmx_03_resort" => "event:/music/remix/03_resort",
    "music_rmx_04_cliffside" => "event:/music/remix/04_cliffside",
    "music_rmx_05_mirror_temple" => "event:/music/remix/05_mirror_temple",
    "music_rmx_06_reflection" => "event:/music/remix/06_reflection",
    "music_rmx_07_summit" => "event:/music/remix/07_summit",
    "music_rmx_09_core" => "event:/music/remix/09_core",
    "cas_01_forsaken_city" => "event:/music/cassette/01_forsaken_city",
    "cas_02_old_site" => "event:/music/cassette/02_old_site",
    "cas_03_resort" => "event:/music/cassette/03_resort",
    "cas_04_cliffside" => "event:/music/cassette/04_cliffside",
    "cas_05_mirror_temple" => "event:/music/cassette/05_mirror_temple",
    "cas_06_reflection" => "event:/music/cassette/06_reflection",
    "cas_07_summit" => "event:/music/cassette/07_summit",
    "cas_08_core" => "event:/music/cassette/09_core",
    "music_farewell_part01" => "event:/new_content/music/lvl10/part01",
    "music_farewell_part02" => "event:/new_content/music/lvl10/part02",
    "music_farewell_part03" => "event:/new_content/music/lvl10/part03",
    "music_farewell_intermission_heartgroove" => "event:/new_content/music/lvl10/intermission_heartgroove",
    "music_farewell_intermission_powerpoint" => "event:/new_content/music/lvl10/intermission_powerpoint",
    "music_farewell_reconciliation" => "event:/new_content/music/lvl10/reconciliation",
    "music_farewell_cassette" => "event:/new_content/music/lvl10/cassette_rooms",
    "music_farewell_final_run" => "event:/new_content/music/lvl10/final_run",
    "music_farewell_end_cinematic" => "event:/new_content/music/lvl10/cinematic/end",
    "music_farewell_end_cinematic_intro" => "event:/new_content/music/lvl10/cinematic/end_intro",
    "music_farewell_firstbirdcrash_cinematic" => "event:/new_content/music/lvl10/cinematic/bird_crash_first",
    "music_farewell_secondbirdcrash_cinematic" => "event:/new_content/music/lvl10/cinematic/bird_crash_second",
    "music_farewell_granny" => "event:/new_content/music/lvl10/granny_farewell",
    "music_farewell_golden_room" => "event:/new_content/music/lvl10/golden_room",
))
end

module CassetteSongs
macro songnames(s::Symbol, expr::Expr)
    return Expr(:block, esc(:($s = $expr)), [esc(:($(Meta.parse(k)) = $v)) for (k, v) in eval(expr)]...)
end

@songnames(songs, Dict{String, String}(
    "cas_01_forsaken_city" => "event:/music/cassette/01_forsaken_city",
    "cas_02_old_site" => "event:/music/cassette/02_old_site",
    "cas_03_resort" => "event:/music/cassette/03_resort",
    "cas_04_cliffside" => "event:/music/cassette/04_cliffside",
    "cas_05_mirror_temple" => "event:/music/cassette/05_mirror_temple",
    "cas_06_reflection" => "event:/music/cassette/06_reflection",
    "cas_07_summit" => "event:/music/cassette/07_summit",
    "cas_08_core" => "event:/music/cassette/09_core",
    "music_farewell_cassette" => "event:/new_content/music/lvl10/cassette_rooms"
))
end

module EnvironmentSounds
macro soundnames(s::Symbol, expr::Expr)
    return Expr(:block, esc(:($s = $expr)), [esc(:($(Meta.parse(k)) = $v)) for (k, v) in eval(expr)]...)
end

@soundnames(sounds, Dict{String, String}(
    "env_loc_02_lamp" => "event:/env/local/02_old_site/phone_lamp",
    "env_loc_03_brokenwindow_large_loop" => "event:/env/local/03_resort/broken_window_large",
    "env_loc_03_brokenwindow_small_loop" => "event:/env/local/03_resort/broken_window_small",
    "env_loc_03_pico8machine_loop" => "event:/env/local/03_resort/pico8_machine",
    "env_loc_07_flag_flap" => "event:/env/local/07_summit/flag_flap",
    "env_loc_09_conveyer_idle" => "event:/env/local/09_core/conveyor_idle",
    "env_loc_09_fireball_idle" => "event:/env/local/09_core/fireballs_idle",
    "env_loc_09_lavagate_idle" => "event:/env/local/09_core/lavagate_idle",
    "env_loc_campfire_loop" => "event:/env/local/campfire_loop",
    "env_loc_campfire_start" => "event:/env/local/campfire_start",
    "env_loc_waterfall_big_in" => "event:/env/local/waterfall_big_in",
    "env_loc_waterfall_big_main" => "event:/env/local/waterfall_big_main",
    "env_loc_waterfall_small_in_deep" => "event:/env/local/waterfall_small_in_deep",
    "env_loc_waterfall_small_in_shallow" => "event:/env/local/waterfall_small_in_shallow",
    "env_loc_waterfall_small_main" => "event:/env/local/waterfall_small_main",
    "env_loc_10_cafe_computer" => "event:/new_content/env/local/cafe_computer",
    "env_loc_10_cafe_sign" => "event:/new_content/env/local/cafe_sign",
    "env_loc_10_tutorial_static_left" => "event:/new_content/env/local/tutorial_static_left",
    "env_loc_10_tutorial_static_right" => "event:/new_content/env/local/tutorial_static_right",
    "env_loc_10_kevinpc" => "event:/new_content/env/local/kevinpc",
))
end

module AmbientSounds
macro soundnames(s::Symbol, expr::Expr)
    return Expr(:block, esc(:($s = $expr)), [esc(:($(Meta.parse(k)) = $v)) for (k, v) in eval(expr)]...)
end

@soundnames(sounds, Dict{String, String}(
    "env_amb_00_main" => "event:/env/amb/00_prologue",
    "env_amb_01_main" => "event:/env/amb/01_main",
    "env_amb_02_awake" => "event:/env/amb/02_awake",
    "env_amb_02_dream" => "event:/env/amb/02_dream",
    "env_amb_03_exterior" => "event:/env/amb/03_exterior",
    "env_amb_03_interior" => "event:/env/amb/03_interior",
    "env_amb_03_pico8_closeup" => "event:/env/amb/03_pico8_closeup",
    "env_amb_04_main" => "event:/env/amb/04_main",
    "env_amb_05_interior_dark" => "event:/env/amb/05_interior_dark",
    "env_amb_05_interior_main" => "event:/env/amb/05_interior_main",
    "env_amb_05_mirror_sequence" => "event:/env/amb/05_mirror_sequence",
    "env_amb_06_lake" => "event:/env/amb/06_lake",
    "env_amb_06_main" => "event:/env/amb/06_main",
    "env_amb_06_prehug" => "event:/env/amb/06_prehug",
    "env_amb_09_main" => "event:/env/amb/09_main",
    "env_amb_worldmap" => "event:/env/amb/worldmap",
    "env_amb_10_rain" => "event:/new_content/env/10_rain",
    "env_amb_10_electricity" => "event:/new_content/env/10_electricity",
    "env_amb_10_endscene" => "event:/new_content/env/10_endscene",
    "env_amb_10_rushingvoid" => "event:/new_content/env/10_rushingvoid",
    "env_amb_10_space_underwater" => "event:/new_content/env/10_space_underwater",
    "env_amb_10_voidspiral" => "event:/new_content/env/10_voidspiral",
    "env_amb_10_grannyclouds" => "event:/new_content/env/10_grannyclouds",
))
end

const tileset_sound_ids = Dict{String, Int}(
    "Default" => -1,
    "Null" => 0,
    "Asphalt" => 1,
    "Car" => 2,
    "Dirt" => 3,
    "Snow" => 4,
    "Wood" => 5,
    "Bridge" => 6,
    "Girder" => 7,
    "Brick" => 8,
    "Zip Mover" => 9,
    "Space Jam (Inactive)" => 11,
    "Space Jam (Active)" => 12,
    "Resort Wood" => 13,
    "Resort Roof" => 14,
    "Resort Platform" => 15,
    "Resort Basement" => 16,
    "Resort Laundry" => 17,
    "Resort Boxes" => 18,
    "Resort Books" => 19,
    "Resort Forcefield" => 20,
    "Resort Clutterswitch" => 21,
    "Resort Elevator" => 22,
    "Cliffside Snow" => 23,
    "Cliffside Grass" => 25,
    "Cliffside Whiteblock" => 27,
    "Gondola" => 28,
    "Glass" => 32,
    "Grass" => 33,
    "Cassette Block" => 35,
    "Core Ice" => 36,
    "Core Rock" => 37,
    "Glitch" => 40,
    "Internet Café" => 42,
    "Cloud" => 43,
    "Moon" => 44
)

const wind_patterns = String[
    "None",
    "Left",
    "Right",
    "LeftStrong",
    "RightStrong",
    "LeftOnOff",
    "RightOnOff",
    "LeftOnOffFast",
    "RightOnOffFast",
    "Alternating",
    "LeftGemsOnly",
    "RightCrazy",
    "Down",
    "Up",
    "Space"
]

const temple_gate_modes = String[
    "NearestSwitch",
    "CloseBehindPlayer",
    "CloseBehindPlayerAlways",
    "HoldingTheo",
    "TouchSwitches",
    "CloseBehindPlayerAndTheo"
]

const bird_npc_modes = String[
    "ClimbingTutorial",
    "DashingTutorial",
    "DreamJumpTutorial",
    "SuperWallJumpTutorial",
    "HyperJumpTutorial",
    "FlyAway",
    "Sleeping",
    "MoveToNodes",
    "WaitForLightningOff",
    "None"
]

const everest_bird_tutorial_tutorials = String[
    "TUTORIAL_CLIMB",
    "TUTORIAL_HOLD",
    "TUTORIAL_DASH",
    "TUTORIAL_DREAMJUMP",
    "TUTORIAL_CARRY",
    "hyperjump/tutorial00",
    "hyperjump/tutorial01"
]

const bonfire_modes = String[
    "Unlit",
    "Lit",
    "Smoking"
]

const clutter_block_colors = String[
    "Red",
    "Green",
    "Yellow"
]

# Core and Rainbow color added by Everest
const crystal_colors = String[
    "Blue",
    "Red",
    "Purple",
    "Core",
    "Rainbow"
]

const dash_switch_sides = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const tentacle_sides = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const tentacle_fear_distance = String[
    "close",
    "medium",
    "far"
]

const move_block_directions = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const spike_directions = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const tentacle_effect_directions = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const planet_effect_sizes = String[
    "big",
    "small"
]

const trigger_spike_directions = String[
    "Up",
    "Down",
    "Left",
    "Right"
]

const spring_orientations = String[
    "Floor",
    "WallLeft",
    "WallRight"
]

const track_spinner_speeds = String[
    "Slow",
    "Normal",
    "Fast"
]

const fake_wall_modes = String[
    "Wall",
    "Block"
]

const condition_block_conditions = String[
    "Key",
    "Button",
    "Strawberry"
]

const spike_types = String[
    "default",
    "outline",
    "cliffside",
    "tentacles",
    "reflection"
]

const crumble_block_textures = String[
    "default",
    "cliffside"
]

const wood_platform_textures = String[
    "default",
    "cliffside"
]

const kevin_axes = String[
    "both",
    "horizontal",
    "vertical"
]

const npc_npcs = String[
    "granny_00_house",
    "theo_01_campfire",
    "theo_02_campfire",
    "theo_03_escaping",
    "theo_03_vents",
    "oshiro_03_lobby",
    "oshiro_03_hallway",
    "oshiro_03_hallway2",
    "oshiro_03_bigroom",
    "oshiro_03_breakdown",
    "oshiro_03_suite",
    "oshiro_03_rooftop",
    "granny_04_cliffside",
    "theo_04_cliffside",
    "theo_05_entrance",
    "theo_05_inmirror",
    "evil_05",
    "theo_06_plateau",
    "granny_06_intro",
    "badeline_06_crying",
    "granny_06_ending",
    "theo_06_ending",
    "granny_07x",
    "theo_08_inside",
    "granny_08_inside",
    "granny_09_outside",
    "granny_09_inside",
    "gravestone_10",
    "granny_10_never"
]

const seeker_statue_hatches = String[
    "Distance",
    "PlayerRightOfX"
]

const badeline_boss_shooting_patterns = Int[
    0, 1, 2, 3, 4,
    5, 6, 7, 8, 9,
    10, 11, 12, 13, 14,
    15
]

const tile_entity_legal_tiles = Char[
    '1', '3', '4', '5', '6', '7',
    '8', '9', 'a', 'b', 'c', 'd',
    'e', 'f', 'g', 'h', 'i', 'j',
    'k', 'l'
]

const slider_surfaces = String[
    "Ceiling",
    "LeftWall",
    "RightWall",
    "Floor"
]

const trigger_position_modes = String[
    "HorizontalCenter",
    "VerticalCenter",
    "TopToBottom",
    "BottomToTop",
    "LeftToRight",
    "RightToLeft",
    "NoEffect"
]

const event_trigger_events = String[
    "end_city",
    "end_oldsite_dream",
    "end_oldsite_awake",
    "ch5_see_theo",
    "ch5_found_theo",
    "ch5_mirror_reflection",
    "cancel_ch5_see_theo",
    "ch6_boss_intro",
    "ch6_reflect",
    "ch7_summit",
    "ch8_door",
    "ch9_goto_the_future",
    "ch9_goto_the_past",
    "ch9_moon_intro",
    "ch9_hub_intro",
    "ch9_hub_transition_out",
    "ch9_badeline_helps",
    "ch9_farewell",
    "ch9_ending",
    "ch9_end_golden",
    "ch9_final_room",
    "ch9_ding_ding_ding",
    "ch9_golden_snapshot"
]

const mini_textbox_trigger_modes = String[
    "OnPlayerEnter",
    "OnLevelStart",
    "OnTheoEnter"
]

const everest_flag_trigger_modes = String[
    "OnPlayerEnter",
    "OnPlayerLeave",
    "OnLevelStart"
]

const everest_crystal_shatter_trigger_modes = String[
    "All",
    "Contained"
]

const everest_fake_heart_colors = String[
    "Normal",
    "BSide",
    "CSide",
    "Random"
]

const music_fade_trigger_directions = String[
    "leftToRight",
    "topToBottom"
]

const black_hole_trigger_strengths = String[
    "Mild",
    "Medium",
    "High",
    "Wild"
]

const moon_glitch_background_trigger_durations = String[
    "Short",
    "Medium",
    "Long"
]

const spawn_facing_trigger_facings = String[
    "Right",
    "Left"
]

const zip_mover_themes = String[
    "Normal",
    "Moon"
]

const swap_block_themes = String[
    "Normal",
    "Moon"
]

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
names = ["FormDefine", "Form", "FormCategory", "FormInstance", "FileItem"]
for it in names
    print("Query<$it> query$it()\n")
end

max_stars_repo_name: itang/launch-speed-test

max_issues_repo_name: itang/launch-speed-test

max_forks_repo_name: itang/launch-speed-test


text:
#=

  Gaussian mixture model.

  cf ~/blog/gaussian_mixture_model.blog
     ~/psi/gaussian_mixture_model.psi
     ~/webppl/gaussian_mixture_model.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function gaussian_mixture_model()
    x = [0.2,1.0,0.5,0.6]
    n = length(x)

    p ~ Beta(0.5, 1)
    a ~ Uniform(-1, 1)
    b ~ Uniform(-1, 1)

    z = tzeros(n)
    x = tzeros(n)
    for i in 1:n
        z[i] ~ Bernoulli(p)
        if z[i] == 1
            x[i] ~ Normal(a, 1.0)
        else
            x[i] ~ Normal(b, 1.0)
        end
    end

    a_is_greater_than_b ~ Dirac(a > b)

end

model = gaussian_mixture_model()

num_chns = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), MCMCThreads(), 100_000, num_chns)
# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), 100_000)

# chns = sample(model, PG(15), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
chns = sample(model, SMC(), 40_000)
# chns = sample(model, IS(), 10_000)

#
# chns = sample(model, Gibbs(HMC(0.1,5,:a,:b),PG(15,:p,:a_is_greater_than_b)), 10_000)
# chns = sample(model, Gibbs(NUTS(1000,0.65,:a,:b),PG(15,:p,:a_is_greater_than_b)), 10_000)

display(chns)
# display(plot(chns))

# show_var_dist_pct(chns,:a,20)
# show_var_dist_pct(chns,:b,20)
# show_var_dist_pct(chns,:p,20)

show_var_dist_pct(chns,:a_is_greater_than_b)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module WebSockets

using HttpCommon
using HttpServer
using Codecs
using Nettle
export WebSocket,
       WebSocketHandler,
       write,
       read,
       close,
       send_ping,
       send_pong

import Base: read, write, isopen, close

# A WebSocket is just a wrapper over a TcpSocket
# All it does is wrap outgoing data in the protocol
# and unwrapp incoming data.
type WebSocket
  id::Int
  socket::Base.TcpSocket
  is_closed::Bool
  sent_close::Bool

  function WebSocket(id::Int,socket::Base.TcpSocket)
    new(id,socket, !isopen(socket), false)
  end
end

#
# WebSocket Packets
#

#      0                   1                   2                   3
#      0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
#     +-+-+-+-+-------+-+-------------+-------------------------------+
#     |F|R|R|R| opcode|M| Payload len |    Extended payload length    |
#     |I|S|S|S|  (4)  |A|     (7)     |             (16/64)           |
#     |N|V|V|V|       |S|             |   (if payload len==126/127)   |
#     | |1|2|3|       |K|             |                               |
#     +-+-+-+-+-------+-+-------------+ - - - - - - - - - - - - - - - +
#     |     Extended payload length continued, if payload len == 127  |
#     + - - - - - - - - - - - - - - - +-------------------------------+
#     |                               |Masking-key, if MASK set to 1  |
#     +-------------------------------+-------------------------------+
#     | Masking-key (continued)       |          Payload Data         |
#     +-------------------------------- - - - - - - - - - - - - - - - +
#     :                     Payload Data continued ...                :
#     + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +
#     |                     Payload Data continued ...                |
#     +---------------------------------------------------------------+
#

# Opcode values
#  *  %x0 denotes a continuation frame
#  *  %x1 denotes a text frame
#  *  %x2 denotes a binary frame
#  *  %x3-7 are reserved for further non-control frames
#
#  *  %x8 denotes a connection close
#  *  %x9 denotes a ping
#  *  %xA denotes a pong
#  *  %xB-F are reserved for further control frames

# Internal function for wrapping one
# piece of data into a WS header
# sending it out over the TcpSocket.
function send_fragment(ws::WebSocket, islast::Bool, data::Array{Uint8}, opcode=0b0001)
  l = length(data)
  b1::Uint8 = (islast ? 0b1000_0000 : 0b0000_0000) | opcode
  if l <= 125
    write(ws.socket,b1)
    write(ws.socket,uint8(l))
    write(ws.socket,data)
  elseif l <= typemax(Uint16)
    write(ws.socket,b1)
    write(ws.socket,uint8(126))
    write(ws.socket,hton(uint16(l)))
    write(ws.socket,data)
  elseif l <= typemax(Uint64)
    write(ws.socket,b1)
    write(ws.socket,uint8(127))
    write(ws.socket,hton(uint64(l)))
    write(ws.socket,data)
  else
    error("Attempted to send too much data for one websocket fragment\n")
  end
end
send_fragment(ws::WebSocket, islast::Bool, data::ByteString, opcode=0b0001) =
  send_fragment(ws, islast, data.data, opcode)

# Exported function for sending data into a websocket
# data should allow length(data) and write(TcpSocket,data)
# all protocol details are taken care of.
function write(ws::WebSocket,data::ByteString)
  if ws.is_closed
    @show ws
    error("attempt to write to closed WebSocket\n")
  end

  #assume data fits in one fragment
  send_fragment(ws,true,data)
end

function write(ws::WebSocket, data::Array{Uint8})
  if ws.is_closed
    @show ws
    error("attempt to write to closed WebSocket\n")
  end

  #indicate this is binary data
  send_fragment(ws,true,data, 0b0010)
end

function send_ping(ws::WebSocket, data = "")
  send_fragment(ws, true, data, 0x9)
end

function send_pong(ws::WebSocket, data = "")
  send_fragment(ws, true, data, 0xA)
end

function close(ws::WebSocket)
    send_fragment(ws, true, "", 0x8)
    ws.is_closed = true
    while true
      wsf = read_frame(ws)
      is_control_frame(wsf) || continue
      wsf.opcode == 0x8 || continue
      break
    end
    close(ws.socket)
end

isopen(ws::WebSocket) = !ws.is_closed

# represents one received message fragment
# (headers + data)
type WebSocketFragment
  is_last::Bool
  rsv1::Bool
  rsv2::Bool
  rsv3::Bool
  opcode::Uint8 #really, Uint4
  is_masked::Bool
  payload_len::Uint64
  maskkey::Vector{Uint8} #Union{Array{Uint8,4}, Nothing}
  data::Vector{Uint8} #ByteString
end

# constructor to do some conversions from bits to Bool.
function WebSocketFragment(
   fin::Uint8
  ,rsv1::Uint8
  ,rsv2::Uint8
  ,rsv3::Uint8
  ,opcode::Uint8
  ,masked::Uint8
  ,payload_len::Uint64
  ,maskkey::Vector{Uint8}
  ,data::Vector{Uint8})

  WebSocketFragment(
      bool(fin)
    , bool(rsv1)
    , bool(rsv2)
    , bool(rsv3)
    , opcode
    , bool(masked)
    , payload_len
    , maskkey
    , data)
end

# A message frame/fragment can be
# either a control frame or a data frame
# this function determines which it is
# according to the opcode.
# if that bit is set (1), then this is a control frame.
is_control_frame(msg::WebSocketFragment) = (msg.opcode & 0b0000_1000) > 0

function handle_control_frame(ws::WebSocket,wsf::WebSocketFragment)
  if wsf.opcode == 0x8 #  %x8 denotes a connection close
    send_fragment(ws, true, "", 0x8)
    ws.is_closed = true
    wait(ws.socket.closenotify)
  elseif wsf.opcode == 0x9 #  %x9 denotes a ping
    send_pong(ws,wsf.data)
  elseif wsf.opcode == 0xA #  %xA denotes a pong
    # Nothing to do here
  else #  %xB-F are reserved for further control frames
    error("Unknown opcode $(wsf.opcode)")
  end
end

function read_frame(ws::WebSocket)
  a = read(ws.socket,Uint8)
  fin    = a & 0b1000_0000 >>> 7 #if fin, then is final fragment
  rsv1   = a & 0b0100_0000 #if not 0, fail.
  rsv2   = a & 0b0010_0000 #if not 0, fail.
  rsv3   = a & 0b0001_0000 #if not 0, fail.
  opcode = a & 0b0000_1111 #if not known code, fail.
  #TODO: add validation somewhere to ensure rsv,opcode,mask,etc are valid.

  b = read(ws.socket,Uint8)
  mask = b & 0b1000_0000 >>> 7 #if not 1, fail.

  payload_len::Uint64 = b & 0b0111_1111
  if payload_len == 126
    payload_len = ntoh(read(ws.socket,Uint16)) #2 bytes
  elseif payload_len == 127
    payload_len = ntoh(read(ws.socket,Uint64)) #8 bytes
  end

  maskkey = Array(Uint8,4)
  for i in 1:4
   maskkey[i] = read(ws.socket,Uint8)
  end

  data = Array(Uint8, payload_len)
  for i in 1:payload_len
    d = read(ws.socket, Uint8)
    d = d $ maskkey[mod(i - 1, 4) + 1]
    data[i] = d
  end

  return WebSocketFragment(fin,rsv1,rsv2,rsv3,opcode,mask,payload_len,maskkey,data)
end

import Base.read
# Read data from a WebSocket.
# This will block until a full message has been received.
# The headers will be stripped and only the data will be returned.
function read(ws::WebSocket)
  if ws.is_closed
    error("Attempt to read from closed WebSocket")
  end

  frame = read_frame(ws)

  #handle control (non-data) messages
  if is_control_frame(frame)
    @show handle_control_frame(ws,frame)
    return read(ws)
  end

  #handle data that uses multiple fragments
  if !frame.is_last
    return concatenate(frame.data,read(ws))
  end

  return frame.data
end

#
# WebSocket Handshake
#

# get key out of request header
get_websocket_key(request::Request) = begin
  return request.headers["Sec-WebSocket-Key"]
end

# the protocol requires that a special key
# be processed and sent back with the handshake response
# to prove that received the HTTP request
# and that we *really* know what webSockets means.
function generate_websocket_key(key)
  h = HashState(SHA1)
  update!(h, key*"258EAFA5-E914-47DA-95CA-C5AB0DC85B11")
  bytestring(encode(Base64, digest!(h)))
end

# perform the handshake assuming it's a websocket request
websocket_handshake(request,client) = begin
  key = get_websocket_key(request)
  resp_key = generate_websocket_key(key)

  #TODO: use a proper HTTP response type
  response = "HTTP/1.1 101 Switching Protocols\r\nUpgrade: websocket\r\nConnection: Upgrade\r\nSec-WebSocket-Accept: "
  Base.write(client.sock,"$response$resp_key\r\n\r\n")
end

# Implement the WebSocketInterface
# so that this implementation can be used
# in Http's server implementation.
immutable WebSocketHandler <: HttpServer.WebSocketInterface
    handle::Function
end

import HttpServer: handle, is_websocket_handshake
function handle(handler::WebSocketHandler, req::Request, client::HttpServer.Client)
    websocket_handshake(req, client)
    sock = WebSocket(client.id, client.sock)
    handler.handle(req, sock)
    isopen(sock) && close(sock)
end
function is_websocket_handshake(handler::WebSocketHandler, req::Request)
    get(req.headers, "Upgrade", false) == "websocket"
end

end # module WebSockets


text:
@testset "deprecated.jl" begin
    @testset "calibrationerror" begin
        ece = ECE(UniformBinning(10))
        skce1 = UnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())
        skce2 = BlockUnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())
        skce3 = BiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())
        ucme = UCME(
            SqExponentialKernel() ⊗ WhiteKernel(),
            [[1, 0, 0, 0, 0], [0, 0, 0, 0, 1]],
            [2, 1],
        )

        predictions = map(1:10) do _
            x = rand(5)
            x ./= sum(x)
            return x
        end
        targets = rand(1:5, 10)

        for estimator in (ece, skce1, skce2, skce3, ucme)
            estimate = estimator(predictions, targets)
            @test @test_deprecated(calibrationerror(estimator, predictions, targets)) ==
                  estimate
            @test @test_deprecated(calibrationerror(estimator, (predictions, targets))) ==
                  estimate
            @test @test_deprecated(
                calibrationerror(estimator, reduce(hcat, predictions), targets)
            ) == estimate
            @test @test_deprecated(
                calibrationerror(estimator, map(tuple, predictions, targets))
            ) == estimate
        end
    end

    @testset "BiasedSKCE" begin
        skce = BiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())

        skce2 = @test_deprecated BiasedSKCE(SqExponentialKernel(), WhiteKernel())
        @test typeof(skce2) === typeof(skce)
        @test skce2.kernel == skce.kernel
    end

    @testset "UnbiasedSKCE" begin
        skce = UnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())

        skce2 = @test_deprecated UnbiasedSKCE(SqExponentialKernel(), WhiteKernel())
        @test typeof(skce2) === typeof(skce)
        @test skce2.kernel == skce.kernel
    end

    @testset "BlockUnbiasedSKCE" begin
        skce = BlockUnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())

        skce2 = @test_deprecated BlockUnbiasedSKCE(SqExponentialKernel(), WhiteKernel())
        @test typeof(skce2) === typeof(skce)
        @test skce2.kernel == skce.kernel
        @test skce2.blocksize == skce.blocksize
        skce3 = @test_deprecated BlockUnbiasedSKCE(SqExponentialKernel(), WhiteKernel(), 2)
        @test typeof(skce3) === typeof(skce)
        @test skce3.kernel == skce.kernel
        @test skce3.blocksize == skce.blocksize
    end

    @testset "UCME" begin
        # three test locations
        ucme = UCME(
            SqExponentialKernel() ⊗ WhiteKernel(),
            [[1.0, 0], [0.5, 0.5], [0.0, 1]],
            [1, 1, 2],
        )

        # Deprecations
        ucme2 = @test_deprecated UCME(
            SqExponentialKernel(),
            WhiteKernel(),
            [[1.0, 0], [0.5, 0.5], [0.0, 1]],
            [1, 1, 2],
        )
        @test typeof(ucme2) === typeof(ucme)
        @test ucme2.kernel == ucme.kernel
        @test ucme2.testpredictions == ucme.testpredictions
        @test ucme2.testtargets == ucme.testtargets

        ucme3 = @test_deprecated UCME(
            SqExponentialKernel() ⊗ WhiteKernel(), [1.0 0.5 0.0; 0 0.5 1], [1, 1, 2]
        )
        @test typeof(ucme3) === typeof(ucme)
        @test ucme3.kernel == ucme.kernel
        @test ucme3.testpredictions == ucme.testpredictions
        @test ucme3.testtargets == ucme.testtargets

        ucme4 = @test_deprecated UCME(
            SqExponentialKernel() ⊗ WhiteKernel(),
            ([[1.0, 0], [0.5, 0.5], [0.0, 1]], [1, 1, 2]),
        )
        @test typeof(ucme4) === typeof(ucme)
        @test ucme4.kernel == ucme.kernel
        @test ucme4.testpredictions == ucme.testpredictions
        @test ucme4.testtargets == ucme.testtargets

        ucme5 = @test_deprecated UCME(
            SqExponentialKernel() ⊗ WhiteKernel(),
            [([1.0, 0], 1), ([0.5, 0.5], 1), ([0.0, 1], 2)],
        )
        @test typeof(ucme5) === typeof(ucme)
        @test ucme5.kernel == ucme.kernel
        @test ucme5.testpredictions == ucme.testpredictions
        @test ucme5.testtargets == ucme.testtargets
    end

    @testset "predictions_targets" begin
        predictions = map(1:10) do _
            x = rand(5)
            x ./= sum(x)
            return x
        end
        targets = rand(1:5, 10)
        data = (predictions, targets)

        @test @test_deprecated(
            CalibrationErrors.predictions_targets(predictions, targets)
        ) == data
        @test @test_deprecated(
            CalibrationErrors.predictions_targets((predictions, targets))
        ) == data
        @test @test_deprecated(
            CalibrationErrors.predictions_targets(reduce(hcat, predictions), targets)
        ) == data
        @test @test_deprecated(
            CalibrationErrors.predictions_targets(map(tuple, predictions, targets))
        ) == data
    end

    @testset "unbiasedskce" begin
        kernel = SqExponentialKernel() ⊗ WhiteKernel()
        predictions = map(1:10) do _
            x = rand(5)
            x ./= sum(x)
            return x
        end
        targets = rand(1:5, 10)
        @test @test_deprecated(
            CalibrationErrors.unbiasedskce(kernel, predictions, targets)
        ) == UnbiasedSKCE(kernel)(predictions, targets)
    end

    @testset "TVExponentialKernel" begin
        kernel = @test_deprecated(TVExponentialKernel())
        @test kernel isa ExponentialKernel{TotalVariation}
    end
end


text:
"""
Returns true if there was a collision between original and sprite
"""
collision(original, sprite) = (original & sprite) != 0

"""
Return the single bits of line as an iterator
"""
decode(line::Integer) = Iterators.reverse(digits(line, base = 2, pad = sizeof(line)*8))


text:
module AxisKeys

include("struct.jl")
export KeyedArray, axiskeys

include("lookup.jl")

include("names.jl")
export NamedDimsArray, dimnames

include("wrap.jl")
export wrapdims

include("selectors.jl")
export Near, Index, Interval, Not, Key

include("functions.jl")
export sortkeys

include("push.jl")

include("broadcast.jl")

include("show.jl")

include("tables.jl") # Tables.jl

include("stack.jl") # LazyStack.jl

include("fft.jl") # AbstractFFTs.jl

end


text:
##############
# SumOfTerms #
##############

struct SumOfLocalTerms{N, T1<:NTuple{N,ScalarFunction},
                            T2<:NTuple{N,LocalOperator}} <: LocalOperator
    coefficients::T1
    operators::T2
end

operators(op::SumOfLocalTerms) = op.operators
coefficients(op::SumOfLocalTerms) = op.coefficients

Base.:*(α::ScalarFunction, op::LocalOperator) = SumOfLocalTerms((α,), (op,))
Base.:*(op::LocalOperator, α::ScalarFunction) = *(α, op)

Base.:*(α::ScalarFunction, op::SumOfLocalTerms) =
    SumOfLocalTerms(α .* op.coefficients, op.operators)
Base.:*(op::SumOfLocalTerms, α::ScalarFunction) = *(α, op)

Base.:-(op::LocalOperator) = SumOfLocalTerms((-1,), (op,))
Base.:-(op::SumOfLocalTerms) = SumOfLocalTerms(.-(op.coefficients), op.operators)

Base.:+(op1::LocalOperator, op2::LocalOperator) = SumOfLocalTerms((1,1), (op1, op2))
Base.:-(op1::LocalOperator, op2::LocalOperator) = SumOfLocalTerms((1,-1), (op1, op2))
Base.:+(op1::SumOfLocalTerms, op2::LocalOperator) =
    SumOfLocalTerms((op1.coefficients..., 1), (op1.operators..., op2))
Base.:-(op1::SumOfLocalTerms, op2::LocalOperator) =
    SumOfLocalTerms((op1.coefficients..., -1), (op1.operators..., op2))
Base.:+(op1::LocalOperator, op2::SumOfLocalTerms) =
    SumOfLocalTerms((1, op2.coefficients...), (op1, op2.operators...))
Base.:-(op1::LocalOperator, op2::SumOfLocalTerms) = +(op1, -(op2))
Base.:+(op1::SumOfLocalTerms, op2::SumOfLocalTerms) =
    SumOfLocalTerms((op1.coefficients..., op2.coefficients...),
                    (op1.operators..., op2.operators...))
Base.:-(op1::SumOfLocalTerms, op2::SumOfLocalTerms) = +(op1, -(op2))

Base.:*(op1::SumOfLocalTerms, op2::LocalOperator) =
    SumOfLocalTerms(op1.coefficients, op1.operators .* (op2,))
Base.:*(op1::LocalOperator, op2::SumOfLocalTerms) =
    SumOfLocalTerms(op2.coefficients, (op1,) .* op2.operators)

Base.:*(op1::SumOfLocalTerms{1}, op2::SumOfLocalTerms) =
    op1.coefficients[1] * (op1.operators[1] * op2)

function Base.:*(op1::SumOfLocalTerms, op2::SumOfLocalTerms)
    op1.coefficients[1] * (op1.operators[1] * op2) +
    SumOfLocalTerms(Base.tail(op1.coefficients), Base.tail(op1.operators)) * op2
end

Base.adjoint(op::SumOfLocalTerms) =
    SumOfLocalTerms(adjoint.(op.coefficients), adjoint.(op.operators))


text:
@testset "LP Atoms: $solver" for solver in solvers
    @testset "abs atom" begin
        x = Variable()
        p = minimize(abs(x), x<=-1)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 1 atol=TOL
        @test evaluate(abs(x)) ≈ 1 atol=TOL

        x = Variable(2,2)
        p = minimize(sum(abs(x)), x[2,2]>=1, x[1,1]>=1, x>=0)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 2 atol=TOL
        @test evaluate(sum(abs(x))) ≈ 2 atol=TOL
    end

    @testset "maximum atom" begin
        x = Variable(10)
        a = shuffle(collect(0.1:0.1:1.0))
        p = minimize(maximum(x), x >= a)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ maximum(a) atol=TOL
        @test evaluate(maximum(x)) ≈ maximum(a) atol=TOL
    end

    @testset "minimum atom" begin
        x = Variable(1)
        a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
        p = maximize(minimum(x), x <= a)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ minimum(a) atol=TOL
        @test evaluate(minimum(x)) ≈ minimum(a) atol=TOL

        x = Variable(4, 4)
        y = Variable(4, 6)
        z = Variable(1)
        c = ones(4, 1)
        d = fill(2.0, (6, 1))
        constraints = [[x y] <= 2, z <= 0, z <= x, 2z >= -1]
        objective = sum(x + z) + minimum(y) + c' * y * d
        p = maximize(objective, constraints)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 130 atol=TOL
        @test (evaluate(objective))[1] ≈ 130 atol=TOL
    end

    @testset "max atom" begin
        x = Variable(10, 10)
        y = Variable(10, 10)
        a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
        b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
        p = minimize(maximum(max(x, y)), [x >= a, y >= b])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        max_a = maximum(a)
        max_b = maximum(b)
        @test p.optval ≈ max(max_a, max_b) atol=10TOL
        @test evaluate(maximum(max(x, y))) ≈ max(max_a, max_b) atol=10TOL
    end

    @testset "min atom" begin
        x = Variable(10, 10)
        y = Variable(10, 10)
        a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
        b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10))
        p = maximize(minimum(min(x, y)), [x <= a, y <= b])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        min_a = minimum(a)
        min_b = minimum(b)
        @test p.optval ≈ min(min_a, min_b) atol=10TOL
        @test evaluate(minimum(min(x, y))) ≈ min(min_a, min_b) atol=10TOL
    end

    @testset "pos atom" begin
        x = Variable(3)
        a = [-2; 1; 2]
        p = minimize(sum(pos(x)), [x >= a, x <= 2])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 3 atol=TOL
        @test evaluate(sum(pos(x))) ≈ 3 atol=TOL
    end

    @testset "neg atom" begin
        x = Variable(3)
        p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 1 atol=TOL
        @test evaluate(sum(neg(x))) ≈ 6 atol=TOL
    end

    @testset "sumlargest atom" begin
        x = Variable(2)
        p = minimize(sumlargest(x, 2), x >= [1; 1])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 2 atol=TOL
        @test evaluate(sumlargest(x, 2)) ≈ 2 atol=TOL

        x = Variable(4, 4)
        p = minimize(sumlargest(x, 3), x >= eye(4), x[1, 1] >= 1.5, x[2, 3] >= 2.1)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 4.6 atol=TOL
        @test evaluate(sumlargest(x, 2)) ≈ 3.6 atol=TOL
    end

    @testset "sumsmallest atom" begin
        x = Variable(4, 4)
        p = minimize(sumlargest(x, 2), sumsmallest(x, 4) >= 1)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 0.5 atol=TOL
        @test evaluate(sumsmallest(x, 4)) ≈ 1 atol=TOL

        x = Variable(3, 2)
        p = maximize(sumsmallest(x, 3), x >= 2, x <= 5, sumlargest(x, 3) <= 12)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 12 atol=TOL
        @test evaluate(sumsmallest(x, 3)) ≈ 12 atol=TOL
    end

    @testset "dotsort atom" begin
        x = Variable(4, 1)
        p = minimize(dotsort(x, [1, 2, 3, 4]), sum(x) >= 7, x >= 0, x <= 2, x[4] <= 1)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 19 atol=TOL
        @test vec(x.value) ≈ [2; 2; 2; 1] atol=TOL
        @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=TOL

        x = Variable(2, 2)
        p = minimize(dotsort(x, [1 2; 3 4]), sum(x) >= 7, x >= 0, x <= 2, x[2, 2] <= 1)
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 19 atol=TOL
        @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=TOL
    end

    @testset "hinge loss atom" begin
        # TODO: @davidlizeng. We should finish this someday.
    end

    @testset "norm inf atom" begin
        x = Variable(3)
        p = minimize(norm_inf(x), [-2 <= x, x <= 1])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 0 atol=TOL
        @test evaluate(norm_inf(x)) ≈ 0 atol=TOL
    end

    @testset "norm 1 atom" begin
        x = Variable(3)
        p = minimize(norm_1(x), [-2 <= x, x <= 1])
        @test vexity(p) == ConvexVexity()
        solve!(p, solver)
        @test p.optval ≈ 0 atol=TOL
        @test evaluate(norm_1(x)) ≈ 0 atol=TOL
    end
end


text:
using Catalyst
using PathWeightSampling: ComplexSystem
using StaticArrays
using Test

sn = @reaction_network begin
    κ, ∅ --> 2L
    λ, L --> ∅
end κ λ

rn = @reaction_network begin
    ρ, L + R --> L + LR
    μ, LR --> R
    ξ, R + CheY --> R + CheYp
    ν, CheYp --> CheY
end ρ μ ξ ν

xn = @reaction_network begin
    δ, CheYp --> CheYp + X
    χ, X --> ∅
end δ χ

u0 = SA[10, 30, 0, 50, 0, 0]
dtimes = 0:0.5:10.0
ps = [5.0, 1.0]
pr = [1.0, 4.0, 1.0, 2.0]
px = [1.0, 1.0]

system = ComplexSystem(sn, rn, xn, u0, ps, pr, px, dtimes)

algorithms = [DirectMCEstimate(128), SMCEstimate(128)]

for algorithm in algorithms
    result = mutual_information(system, algorithm, num_samples = 5)

    for v in result[!, :MutualInformation]
        @test length(v) == length(system.dtimes)
    end
end


text:
using MagneticReadHead
using Test
@testset "only" begin
    @test MagneticReadHead.only([10]) == 10
    @test MagneticReadHead.only((15 for i in 1:1)) == 15

    @test_throws AssertionError MagneticReadHead.only([10,20])
    @test_throws AssertionError MagneticReadHead.only((i for i in 1:10))
end