text:
using StateMachineIterators
using Base.Test
import Iterators: chain
import IJulia

# Test non-iterator interface:

let
    machine = StateMachine(:going_up, :final,
       [(:going_up, x -> x >= 10, :going_down),
        (:going_down, x -> x <= 1, :final)])

    behaviors = Dict(:going_up => x -> x + 1,
                     :going_down => x -> x - 1,
                     :final => x -> x)

    values = []
    x = 1
    state = initial(machine)
    while state != final(machine)
        push!(values, x)
        state = nextstate(machine, state, x)
        x = behaviors[state](x)
    end

    @test all(values .== [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1])

end

# Test iterator interface:

let
    machine = StateMachine(:going_up, :final,
       [(:going_up, x -> x >= 5, :going_down),
        (:going_down, x -> x <= 1, :final)])

    behaviors = Dict(:going_up => x -> x + 1,
                     :going_down => x -> x - 1,
                     :final => x -> x)

    values = []
    x = 1
    for state in StateMachineIterator(machine, () -> x)
        push!(values, x)
        x = behaviors[state](x)
    end

    @test all(values .== [1, 2, 3, 4, 5, 4, 3, 2, 1])
end

let
    x_machine = StateMachine(:going_up, :final,
       [(:going_up, x -> x >= 10, :going_down),
        (:going_down, x -> x <= 1, :final)])

    y_machine = StateMachine(:going_up, :final,
       [(:going_up, y -> y >= 3, :going_down),
        (:going_down, y -> y <= 1, :going_up)])

    behaviors = Dict(:going_up => a -> a + 1,
                     :going_down => a -> a - 1,
                     :final => a -> a)

    values = []
    x = 1
    y = 1
    for (x_state, y_state) in zip(StateMachineIterator(x_machine, () -> x),
                                  StateMachineIterator(y_machine, () -> y))
        push!(values, (x, y))
        x = behaviors[x_state](x)
        y = behaviors[y_state](y)
    end

    @test all(values .== [
        (1, 1),
        (2, 2),
        (3, 3),
        (4, 2),
        (5, 1),
        (6, 2),
        (7, 3),
        (8, 2),
        (9, 1),
        (10, 2),
        (9, 3),
        (8, 2),
        (7, 1),
        (6, 2),
        (5, 3),
        (4, 2),
        (3, 1),
        (2, 2),
        (1, 3)])
end

let
    machine = StateMachine(:going_up, :final,
       [(:going_up, x -> x >= 5, :going_down),
        (:going_down, x -> x <= 1, :final)])

    behaviors = Dict(:going_up => x -> x + 1,
                     :going_down => x -> x - 1,
                     :final => x -> x)

    values = []

    x = 1
    for state in chain(repeated(StateMachineIterator(machine, () -> x), 2)...)
        push!(values, x)
        x = behaviors[state](x)
    end
    @test all(values .== [1, 2, 3, 4, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5, 4, 3, 2, 1])
end

jupyter = IJulia.jupyter
notebook = "../examples/state_machines.ipynb"
kernel = "julia-$(VERSION.major).$(VERSION.minor)"
run(`$jupyter nbconvert --to notebook --execute $notebook --output $notebook --ExecutePreprocessor.kernel_name="$(kernel)"`)


text:
module Muon

using Random
using SparseArrays
import LinearAlgebra: Adjoint

using HDF5
using DataFrames
using CategoricalArrays
using StructArrays
import CompressHashDisplace: FrozenDict
using FileIO

export readh5mu, readh5ad, writeh5mu, writeh5ad, isbacked, update_obs!, update_var!, update!
export AnnData, MuData

import Pkg
# this executes only during precompilation
let
    pkg = Pkg.Types.read_package(joinpath(@__DIR__, "..", "Project.toml"))
    global VERSION = pkg.version
    global NAME = pkg.name * ".jl"
end
MUDATAVERSION = v"0.1.0"
ANNDATAVERSION = v"0.1.0"

include("index.jl")
include("sparsedataset.jl")
include("transposeddataset.jl")
include("hdf5_io.jl")
include("alignedmapping.jl")
include("anndata.jl")
include("mudata.jl")
include("util.jl")

end # module

ext: jl


text:
mutable struct DecisionChannel{A <: AbstractArray} <: AbstractChannel{A}
    decisions::Dict{Int,A}
    cond_take::Condition
    DecisionChannel(decisions::Dict{Int,A}) where A <: AbstractArray = new{A}(decisions, Condition())
end

function put!(channel::DecisionChannel, t, x)
    channel.decisions[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::DecisionChannel, t)
    x = fetch(channel, t)
    delete!(channel.decisions, t)
    return x
end

isready(channel::DecisionChannel) = length(channel.decisions) > 1
isready(channel::DecisionChannel, t) = haskey(channel.decisions, t)

function fetch(channel::DecisionChannel, t)
    wait(channel, t)
    return channel.decisions[t]
end

function wait(channel::DecisionChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

Decisions{A} = RemoteChannel{DecisionChannel{A}}


text:
module MetacommunityDynamics
using DynamicGrids
using DynamicGrids:
    CellRule, SetCellRule, applyrule, Rule, Ruleset, Neighborhood, Moore, Chain
using Dispersal: OutwardsDispersal, DispersalKernel
using Distributions
using Crayons
using EcologicalNetworks

include(joinpath("types.jl"))
export Measurement

include(joinpath("layers", "layer.jl"))
export Layer
include(joinpath("layers", "occupancy.jl"))
export Occupancy
include(joinpath("layers", "biomass.jl"))
export Biomass

include(joinpath("generators", "poissonprocess.jl"))
export generate, PoissonProcess, grid, pointstogrid
include(joinpath("layers", "traits.jl"))
export StaticTraitLayer

include(joinpath("generators", "environment.jl"))
export StaticEnvironmentalLayer

include(joinpath("mechanisms", "_combinemechanisms.jl"))
include(joinpath("mechanisms", "selection.jl"))
export FitnessFunction,
    GaussianFitness,
    ExponentialFitness,
    DensityDependentFitness,
    AbioticSelection,
    BioticSelection

include(joinpath("mechanisms", "colonization.jl"))
export Colonization,
    RandomColonization, SpatiallyExplicitLevinsColonization, LevinsColonization

include(joinpath("mechanisms", "extinction.jl"))
export RandomExtinction, AbioticExtinction

include(joinpath("mechanisms", "eating.jl"))
export Eating,
    FunctionalResponse,
    HollingTypeI,
    HollingTypeII,
    HollingTypeIII,
    LinearFunctionalResponse,
    LotkaVolterra,
    MichaelisMenten,
    CrowleyMartin,
    FoodWebEating

include(joinpath("mechanisms", "mortality.jl"))
export Mortality, LinearMortality

include(joinpath("mechanisms", "dispersal.jl"))
export AdjacentBernoulliDispersal

include(joinpath("mechanisms", "growth.jl"))
export LogisticGrowth

include(joinpath("species", "speciespool.jl"))
export ContinuousSpeciesPool, DiscreteSpeciesPool
export SingleSpecies,
    DiscreteUnipartiteSpeciesPool,
    DiscreteKpartiteSpeciesPool,
    ContinuousUnipartiteSpeciesPool,
    ContinuousKpartiteSpeciesPool
export speciespool, specieslayers, nspecies, species, metaweb
include(joinpath("species", "printing.jl"))

end # module


text:
# Including Plots output into a question can
# be achieved in various ways, as illustrated
# here.

using JuliaBlackBoard
using Plots
gr()

OPEN = JuliaBlackBoard.OPEN(@__FILE__)

OPEN() do io

    ## With Markdown (the default), we include a plot
    ## through the image markup: `![caption]({{{:plot}}})` where
    ## `{{{:plot}}}` is an encoded plot given by passing the
    ## value:
    ## * `q(plot=Plot(p))` with `p` a `Plots` object, or
    ## * `q(plot=File(fname))` where `fname` is the name of a file where the plot is
    ##    written. This can be useful, say, when a different program is used to generate
    ##    the plot (e.g., `RCall`).
    p = plot(sin, 0, 2pi)

    q = mt"""
# include via `Plot` and the markdown syntax for a plot:

![caption]({{{:plot}}})
"""

    question(io, ESS, q(plot=Plot(p)))

    p = plot(cos, 0, 2pi)
    pngfile = tempname() * ".png"
    savefig(p, pngfile)  # save to a file

    q = mt"""
# include a plot which is stored in a file

![caption]({{{:plot}}})
"""
    question(io, ESS, q(plot=File(pngfile)))

    ## With LaTeX there are some subtleties:
    ##
    ## The easiest way to generate pdf output is with the `ltx"""` string macro,
    ## but this won't allow easy substitution into an `\includegraphics`, as the
    ## spaces aren't allowed and `\includegraphics{{{{:plot}}}}` isn't parsed correctly for
    ## substitution.
    ##
    ## So, we have two ways. You can avoid any string macro and use a raw string, then
    ## regular interpolation (via the $) is possible.
    ## Alternatively, we show how to use a different tag in the Mustache template.
    ## For both, you need to explicitly call `LaTeX()`, as below.
    ##
    ## First without a string Maco
    q = """
This is \\LaTeX with inline math (\$\\sin(x)^2\$) and here we have a plot
(note the necessary escaping.)

\\includegraphics[width=4in, height=3in]{$pngfile}
"""
    question(io, ESS, LaTeX(q))

    ## Next we use `MT"""` which parses the tags `<<`, `>>`. We have an extra `{}` to
    ## ensure no HTML substitution is performed.
    q = MT"""
This is cumbersome, as we need different tags for Mustache substitution.

\includegraphics[width=4in, height=3in]{<<{:plot}>>}

This shows how one might do a second round of substitution, for example a variable {{:a}}.

"""

    qq = q(plot=pngfile)
    qqq = JuliaBlackBoard.Mustache.render(qq, a="(a)")
    question(io, ESS, LaTeX(qqq))

    ## Finally, using the `pgfplotsx` backend, plots can be written as `tikz` files for
    ## inclusion. For this the contents of the file are included. (Not the plot object.)
    ## The `Tikz` function makes this straightforward. This substitution allows the use
    ## of the `ltx"""` string macro and the usual substitution:

    pgfplotsx() # use PGFPlotsX backend for Plots
    p = plot(asin, -1/2, 1/2);
    tikzfile = tempname() * ".tikz"
    savefig(p, tikzfile)

    q = ltx"""
A Tikz plot requires the contents of a file to be included, as illustrated here.\\

{{{:plot}}}

"""

    question(io, ESS, q(plot=Tikz(tikzfile)))

end


text:
"""
	version(vstr::String=j4pr_version)

Prints the current J4PR version and ASCII art. The version is
contained in the global constant `j4pr_version`.
"""
function version(vstr::String=j4pr_version)
	commit = ""
	date = ""
	try
		readmehtod = x->nothing
		if v"0.6" <= VERSION < v"0.7-"
			readmethod = x->readstring(x) 			# Julia 0.6
		else
			readmethod = x->read(x,String) 			# Julia 0.7+
		end
		commit = open(`git show --oneline -s`) do x
				readmethod(x)[1:7]
		end

		date = open(`git show -s --format="%ci"`) do x
			readmethod(x)[1:10]
		end
	catch
		commit = "67a5c10+"
		date = "2017-10-27"
	end

	vers ="
 _  _
(_\\/_)                  |  This is a small library and package wrapper written at 0x0α Research.
(_/\\_)                  |  Type \"?j4pr\" for general documentation.
   _ _   _  _____ _ _   |  Look inside src/j4pr.jl for a list of available algorithms.
  | | | | |/____ / ` |  |
  | | |_| | |  | | /-/  |  Version $(vstr) \"The Monolith\" commit: $(commit) ($(date))
 _/ |\\__  | |  | | |    |
|__/    |_|_|  |_|_|    |  License: MIT, view ./LICENSE.md for details.
\n\n"

	print(vers)
end


text:
function accept(obj::AmbiguityAcceptance, arg0::Vector{IntegerLeastSquareSolution})
    return jcall(obj, "accept", IntegerLeastSquareSolution, (Vector{IntegerLeastSquareSolution},), arg0)
end

function number_of_candidates(obj::AmbiguityAcceptance)
    return jcall(obj, "numberOfCandidates", jint, ())
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: December 2020
# Copyright: Georgia Institute of Technology, 2020

export filter_normalization, sol_traj

function filter_normalization(theta::Flux.Params)

    d = deepcopy(theta)

    for (p, q) in zip(theta, d)

        if ndims(p) == 4

            q .= randn(Float32, size(p))

            for i in size(p, 4)
                q[:, :, :, i] .= q[:, :, :, i] / (Flux.norm(q[:, :, :, i]) + 1f-8)
                q[:, :, :, i] .= q[:, :, :, i] * Flux.norm(p[:, :, :, i])
            end

        end
    end

    return d
end

function +(theta::Flux.Params, delta::Flux.Params)

    d = deepcopy(theta)

    for (p, q) in zip(d, delta)

        if ndims(p) > 1
            p .= p + q
        end

    end

    return d
end

function *(alpha::Float32, theta::Flux.Params)

    d = deepcopy(theta)

    for p in d

        if ndims(p) > 1
            p .= alpha*p
        end

    end

    return d
end

function sol_traj(delta, eta, ws, w0)

    delta_ = cat([delta[i] for i=1:2:length(delta)]..., dims=3)[:, :, :, 1]
    eta_ = cat([eta[i] for i=1:2:length(eta)]..., dims=3)[:, :, :, 1]

    A = cat(vec(delta_), vec(eta_), dims=2)
    proj(v) = (A'*A)^(-1) * (A'*vec(v - w0))

    proj_a = []
    proj_b = []

    for i=1:length(ws)
        pa, pb = proj(ws[i])
        push!(proj_a, pa)
        push!(proj_b, pb)
    end

    return proj_a, proj_b
end


text:
using GLM, DataFrames
glm(y ~ 1, DataFrame(y = float(bitrand(10))), Binomial())

n = 2_500_000; srand(1234321)
df2 = DataFrame(x1 = rand(Normal(), n),
                x2 = rand(Exponential(), n),
                ss = pool(rand(DiscreteUniform(50), n)));
mf = ModelFrame(x1 ~ x1 + x2 + ss, df2)
mm = ModelMatrix(mf)
beta = unshift!(rand(Normal(),52), 0.5); # "true" parameter values

## Create linear predictor and mean response
eta = mm.m * beta;
mu = [linkinv(LogitLink(), x) for x in eta]
y = Float64[rand() < μ for μ in mu];        # simulate observed responses

gm6 = glm(mm.m, y, Binomial())
@time glm(mm.m, y, Binomial());
#@profile glm(mm.m, y, Binomial());
#using ProfileView
#ProfileView.view()


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#Function to do indirect multiple shooting in CRTBP
#
#This version:
#   - Fixed endpoints.
#   - Impulsive maneuvers not allowed.
#   - Uses sigmoid control smoothing
#   - Variable step integration.
#   - Uses dual numbers for automatic differentiation of the Jacobian (only the
#     non-zero elements are computed). ~Same speed as finite differences,
#     but accuracy is perfect.
#   - Exploits sparsity of Jacobian for least squares solution. (speeds up
#     matrix inversion by ~10-100x).
#   - Mesh refinement not implemented yet.
#   - When close to converging, uses second-order-correction step to get a
#     "free" iteration
#
#Inputs:
#   - XC_all    = [12 x n_nodes] Initial guess array of states and costates.
#                 Each column is the state and costate at time t[i].
#                 [Nondimensional units of CRTBP]
#   - t_TU      = Guess vector of times [Nondimensional time units of CRTBP]
#   - MU        = μ for CRTBP
#   - DU        = Distance unit for CRTBP [km]
#   - TU        = Time unit for CRTBP [sec]
#   - n_nodes   = Number of nodes
#   - mass      = Initial spacecraft mass [kg]
#   - thrustLimit = Thrust limit [Newtons]
#   - plot_yn   = Boolean, true to plot solution at every iteration.
#   - flag_adjointsOnly = Boolean, true to only change the adjoints (not the
#                 states). This is helpful to run for a few iterations when we
#                 have a good guess of the states, but a poor guess of the
#                 adjoints.
#   - maxIter   = Max number of iterations.
#   - p         = Chooses the control law.
#                 With p=2, minimum energy objective.
#                 With p=1, minimum fuel objective.
#   - rho       = Smoothing term on sigmoid-smoothed control. Usually start with
#                 rho=1 or rho=0.5. Then use continuation to get down to
#                 rho=1e-3 or rho=1e-4.
#
#Outputs: (XC_all, defect, status_flag)
#   - XC_all    = [12 x n_nodes] Solution array of states and costates.
#                 Each column is the state and costate at time t[i].
#                 [Nondimensional units of CRTBP]
#   - defect    = Array of defect constraint violations.
#   - status_flag =
#                 0 for successful convergence.
#                 1 for failure.

function multiShoot_CRTBP_indirect(XC_all, t_TU, MU, DU, TU, n_nodes,
    mass0, thrustLimit, plot_yn, flag_adjointsOnly, maxIter, p, rho)

    ## DEFINE SUB-FUNCTIONS

    function defectCalc(XC_all, t_TU, nstate, n_nodes, odefun, params)
        #Calculate defect constraint violations

        #pre-allocate defect:
        defect1 = zeros(nstate*2, n_nodes-1)

        #pre-allocate integration errors:
        errors = zeros(n_nodes-1)
        for i = 1:n_nodes-1
            #propagate node i forward from time (i) to time (i+1)

            ########### propagate forward:
            x0 = copy(XC_all[:,i])
            tspan = (t_TU[i], t_TU[i+1])

            prob = ODEProblem(odefun, x0, tspan, params)
            sol_for = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)

            ########## defect calc
            defect1[:,i] = sol_for[:,end] - XC_all[:,i+1]

            #save largest error:
            errors[i] = 0.
        end

        #outputs:
        (defect1, errors)
    end

    function jacobianCalc(XC_all, t_TU, nstate, n_nodes, odefun, params)
        #Compute the Jacobian sparsely, with automatic differentation

        #Notes:
        #- The i'th defect vector is affected by XC_all[:,i]
        #- Size/shape of Jacobian is [# defects,  # variables]
        #- uses Automatic Differentation (takes a while to compile, but runtime
        #  performance is good)

        #simple function used by automatic differentiation:
        function f(XC0)
            #need to use 'eltype' because 'ForwardDiff' will call this function with
            #dual numbers and with regular numbers, and it needs to work with either
            #one.
            prob = ODEProblem(odefun, eltype(XC0).(XC0), eltype(XC0).(tspan), eltype(XC0).(params))

            #propagate forward:
            sol_for = DifferentialEquations.solve(prob, Vern8(), reltol=1e-13, abstol=1e-13)[:,end]
        end

        #Build up non-zero elements of sparse Jacobian:
        Jac_temp = zeros(nstate*2*(n_nodes-1), 2*nstate*2)
        tspan = (t_TU[1], t_TU[2]) #elevate scope
        for i = 1:n_nodes-1
            x0 = copy(XC_all[:,i])
            tspan = (t_TU[i], t_TU[i+1])

            #automatic differentiation with dual numbers:
            Jac_AD_ForwardDiff_prop = ForwardDiff.jacobian(f,x0)

            Jac_temp[(i-1)*nstate*2 + 1 : i*nstate*2, :] = hcat(Jac_AD_ForwardDiff_prop, -eye(nstate*2))
        end

        #re-create full Jacobian (band diagonal):
        Jac_full = zeros( nstate*2*(n_nodes-1), n_nodes*nstate*2 )
        for ind = 1:size(Jac_full,1) #loop through rows
            #Which node's defect are we on?
            idx = Int32(ceil(ind/(nstate*2)));

            #offset (in # of columns) for the state-related ones
            idx2 = (idx-1)*(2*nstate);

            #set the appropriate columns of Jac_full from Jac_temp:
            Jac_full[ind, (idx2+1):(idx2+2*2*nstate)] = Jac_temp[ind,:]
        end

        #Remove the first and last states:
        Jac_full[:, 1:nstate] = 0.
        Jac_full[:, (end-2*nstate+1):(end-nstate)] = 0.

        #outputs:
        Jac_full
    end

    function optimizeTraj_OLS(XC_all, t_TU, defect, Jac_full, nstate, n_nodes, odefun, params)
        #Solve the linearized problem with Ordinary Least Squares (makes use of
        #sparse Jacobian)

        #use Ordinary Least Squares to find a solution that satisfies the constraints
        #Remove the endpoints from the problem. Also remove the initial mass from
        #optimization (just leave it fixed at mass0).

        #concatenate all the defect constraints into one vector:
        if nstate == 6
            defect_vec = defect[:]
        elseif nstate == 7
            defect_vec = defect[:][1:end-1] #leave off final mass defect
        end

        #If true, we only adjust the adjoints. This is helpful in the case where we
        #do not have a good guess for the adjoints, but we do have a good guess for
        #the states. If false, adjust the states & adjoints.
        (r,c) = size(Jac_full)
        temp = trues(c)
        if flag_adjointsOnly
            #Remove the columns of Jac_full corresponding to the all states
            idx = 1:nstate
            for ind = 1:(n_nodes-1)
                temp[idx] = false
                idx += 2*nstate
            end

            Jac_full = Jac_full[:, temp]
        end

        #output:
        Jac_sparse = sparse(Jac_full) #create sparse Jacobian (~2-3 ms)
        xc_update = -Jac_sparse \ defect_vec #sparse matrix pseudo-inverse
        xc_update2 = zeros(size(temp))
        xc_update2[temp] = xc_update
        xc_update = reshape(xc_update2, 2*nstate, n_nodes)

        #Second Order Correction step:
        #We only do the SOC when the update is small (necessary for assumptions
        #to be valid)
        if norm(xc_update, Inf) < 1e-1
            #this trick approximates using the second derivatives, but in reality
            #just uses the same first derivatives again. In some cases, it
            #essentially gives us another iteration for ~zero computational cost.
            XC_all_soc = XC_all + xc_update

            #calculate the defects at the first order step:
            (defect, errors) = defectCalc(XC_all_soc, t_TU, nstate, n_nodes, odefun, params)

            #concatenate all the defect constraints into one vector:
            if nstate == 6
                defect_vec = defect[:]
            elseif nstate == 7
                defect_vec = defect[:][1:end-1] #leave off final mass defect
            end

            #solve for SOC step:
            xc_update_soc = -Jac_sparse \ defect_vec #sparse matrix pseudo-inverse

            #apply SOC step:
            xc_update_soc2 = zeros(size(temp))
            xc_update_soc2[temp] = xc_update_soc
            xc_update_soc = reshape(xc_update_soc2, 2*nstate, n_nodes)
            xc_update += xc_update_soc
        end

        #Output:
        xc_update
    end

    function lineSearch(XC_all, xc_update, t_TU, nstate, n_nodes, odefun, params)
        #Very simple line search function. Seeks to minimize the sum of squared
        #defect constraint violations, while forcing some step to be taken.

        alpha = 1.0

        alpha_all = linspace(0.1, 1, 20)

        er = zeros(size(alpha_all))
        cost = zeros(size(alpha_all))

        for ind = 1:length(alpha_all)
            alpha = alpha_all[ind]

            XC_all_trial = XC_all + xc_update*alpha

            #check defects
            (defect, errors) = defectCalc(XC_all_trial, t_TU, nstate, n_nodes, odefun, params)

            er[ind] = sum(defect[:].^2)
        end

        #Output:
        alpha = alpha_all[ er .== minimum(er) ]
        alpha = alpha[1]
    end

################################################################################

    #number of state elements (should be 6 or 7)
    nstate = Int32(size(XC_all, 1)/2)

    #derivatives function to use
    odefun = CRTBP_stateCostate_deriv!
    time_direction = 1.0 #forward
    params = (MU, DU, TU, thrustLimit, mass0, time_direction, p, rho)

    #default to success:
    status_flag = 0

    t0_TU = copy(t_TU[1])
    tf_TU = copy(t_TU[end])
    tau = (t_TU - t0_TU) / (tf_TU - t0_TU) * 2 - 1; #transform t -> tau

    state_0 = XC_all[1:nstate, 1]
    state_f = XC_all[1:nstate, end]

    ############ First nominal run:
    (defect, errors) = defectCalc(XC_all, t_TU, nstate, n_nodes, odefun, params)

    #Iterate until converged
    iterCount = 0;
    er = 1.0 #initialize error
    while er > 1e-10
        iterCount += 1
        if iterCount > maxIter
            warn("Reached max iteration count at ", iterCount," iterations")
            status_flag = 1
            break
        end

        ############ Compute Jacobian for dynamics constraints
        Jac_full = jacobianCalc(XC_all, t_TU, nstate, n_nodes, odefun, params)

        ############ OPTIMIZE LINEAR SUBPROBLEM
        xc_update = optimizeTraj_OLS(XC_all, t_TU, defect, Jac_full, nstate, n_nodes, odefun, params)

        ############ LINE SEARCH
        alpha = 1.0
        #Call line search:
        if iterCount > 3
            alpha = lineSearch(XC_all, xc_update, t_TU, nstate, n_nodes, odefun, params)
        end

        XC_all = XC_all + xc_update * alpha

        ############ PLOT UPDATE
        if plot_yn
            (XC_dense, t_dense) = densify(XC_all, t_TU, params, 100)

            #compute control, or just use zeros (plotting is faster when not
            #drawing control arrows)
            u_all = zeros(3, size(XC_dense,2))
            # for ind = 1:size(XC_dense,2)
            #     u_all[:,ind] = controlLaw_cart(XC_dense[nstate+4:nstate+6,ind], thrustLimit, p, rho, mass0)
            # end

            plotlyjs(legend=false)
            plotTrajPlotly_indirect(XC_dense, u_all, X0_states, Xf_states, 0.2)
            gui()
        end

        #Make sure the end states didn't change:
        XC_all[1:6, 1] = state_0
        XC_all[1:6, end] = state_f

        ############ CHECK UPDATE
        (defect, errors) = defectCalc(XC_all, t_TU, nstate, n_nodes, odefun, params)

        #print progress report
        er = norm(defect[:], Inf)
        @printf "Iter %d. Max defect = %.2e. tf = %.2f days. α = %.3f.\n" iterCount er tf_TU*TU/day alpha
        if er > 1e3
            warn("Not likely to converge. Aborting.")
            iterCount += 100
        end
    end

    if isnan(XC_all[1])
        status_flag = 2
    end

    #outputs:
    (XC_all, defect, status_flag)
end

function plotTrajPlotly_indirect(XC_all, u_all, X0_states, Xf_states, u_scale)
    #Plot in Julia with PlotlyJS package

    Plots.plot(X0_states[1,:], X0_states[2,:], X0_states[3,:], w=3,
        xlabel="X (DU)", ylabel="Y (DU)", zlabel="Z (DU)")
    Plots.plot!(Xf_states[1,:], Xf_states[2,:], Xf_states[3,:], w=3)

    Plots.plot!(XC_all[1,:],XC_all[2,:],XC_all[3,:], w=3, color=:black)

    arrows_start = XC_all[1:3,:]
    arrows_end = XC_all[1:3,:] + u_all*u_scale

    #remove zero thrust vectors:
    arrows_start = arrows_start[:, norm_many(u_all) .!== 0.]
    arrows_end   = arrows_end[:, norm_many(u_all) .!== 0.]

    arrows_x = vcat(arrows_start[1,:]', arrows_end[1,:]')
    arrows_y = vcat(arrows_start[2,:]', arrows_end[2,:]')
    arrows_z = vcat(arrows_start[3,:]', arrows_end[3,:]')

    Plots.plot!(arrows_x, arrows_y, arrows_z, color=:red)

    #invisible points to get axis scaling right
    plot!([1-MU,], [0,], [.1,], marker=1, markeralpha=0.)
    plot!([1-MU,], [0,], [-.1,], marker=1, markeralpha=0.)

    #create data points for a sphere surface:
    (x,y,z) = sphere(32)

    #scale and move the sphere:
    r_moon = 1737. #km
    x = (x * r_moon/DU) + (1-MU)
    y = y * r_moon/DU
    z = z * r_moon/DU

    #plot the sphere as a surface:
    clibrary(:cmocean)
    Plots.surface!(x,y,z, c=:gray)
end

function controlLaw_cart(lambda_v, thrustLimit, p, rho, mass)
    #Calculate control (only used for plotting)

    #Control law for indirect multiple shooting (compute the control as a function
    #of the costates).
    #
    #Valid for any Cartesian representation (inertial or rotating frame).
    #
    #Inputs:
    #     lambda_v = [3 x 1]
    #     thrustLimit
    #     p = exponent to determine which control law to use
    #             (p=2 for min energy, p=1 for min fuel)
    #Outputs:
    #     controls = [3 x 1]

    #magnitude of primer vector:
    lambda_v_mag = norm(lambda_v);

    #direction of thrust, based on primer vector:
    uhat = - lambda_v ./ lambda_v_mag

    accelLimit = thrustLimit / mass / 1e3 * TU^2 / DU #(kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration

    #Thrust magnitude from control law:
    if p == 0
      #Force thrust to always be on, in the optimal direction. This limits how big
      #of a search we need to make for the neural net stuff.
      umag = accelLimit #DU/TU^2

    elseif p == 1 #minimum mass problem
        umag = 1/2 * (1 + tanh( (lambda_v_mag - 1) / (2*rho) )) * accelLimit #DU/TU^2

    elseif p > 1 #p is in the range (1, 2]
        umag = (1/p * norm(lambda_v))^(1 / (p-1)) #DU/TU^2

        if umag > accelLimit
          umag = accelLimit #DU/TU^2
        end
    else
      error("Invalid value of p!")
    end

    if isnan(umag)
        umag = 0.
    end

    #units of umag are acceleration (DU/TU^2)

    #Control vector, scaled properly:
    control = -umag * lambda_v ./ norm(lambda_v) * mass * DU * 1e3 / TU^2 #(DU/TU^2) -> (km/s^2) -> (kg*m/s^2) force
end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
# Include everything
include("md.jl")

# File to load and save to
fname = "../../swp/s_auto.mdv"

# Load trial s_auto
s = loadmdtrial(fname)

# Run for another 10000 steps
@time for i in 1:20
    println("starting another round")
    verlet!(s, 500)
    save(fname, s)
    println("saved round $i.")
end


text:
abstract type AbstractOptimizer end

@params mutable struct Optimizer{Tdev} <: AbstractOptimizer
    model::AbstractModel
    alg
    workspace
    dev::Tdev
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::Optimizer) = println("TopOpt optimizer")

multol!(o::Optimizer, m::Real) = o.options.tol *= m
function setbounds!(o::Optimizer, x, w)
    o.model.box_min .= max.(0, x .- w)
    o.model.box_max .= min.(1, x .+ w)
end

function Optimizer(
    obj,
    constr,
    vars,
    opt = Nonconvex.MMA87(),
    device::Tdev = CPU();
    options = Nonconvex.MMAOptions(),
    convcriteria = Nonconvex.KKTCriteria(),
) where {Tdev <: AbstractDevice}

    T = eltype(vars)
    nvars = length(vars)
    x0 = copy(vars)
    model = Nonconvex.Model(obj)
    addvar!(model, zeros(T, nvars), ones(T, nvars))
    add_ineq_constraint!(model, Nonconvex.FunctionWrapper(constr, length(constr(x0))))
    workspace = Nonconvex.Workspace(model, opt, x0; options = options, convcriteria = convcriteria)
    return Optimizer(model, opt, workspace, device)
end

Utilities.getpenalty(o::Optimizer) = getpenalty(o.solver)
function Utilities.setpenalty!(o::Optimizer, p)
    setpenalty!(o.solver, p)
    return o
end

function (o::Optimizer)(x0::AbstractVector)
    @unpack workspace = o
    @unpack options, model = workspace
    reset!(workspace, x0)
    setoptions!(workspace, options)
    mma_results = Nonconvex.optimize!(workspace)
    return mma_results
end
function setoptions!(workspace::Nonconvex.Workspace, options)
    workspace.options = options
    return workspace
end

function reset!(w::Nonconvex.Workspace, x0 = nothing)
    if x0 !== nothing
        w.x0 .= x0
    end
    return w
end
function reset!(w::Nonconvex.MMAWorkspace, x0 = nothing)
    @unpack solution = w
    outer_iter, iter, fcalls = 0, 0, 0, 0
    @pack! w = fcalls, iter, outer_iter
    if x0 !== nothing
        w.x0 .= x0
        w.tempx .= solution.prevx
        solution.prevx .= solution.x
        solution.x .= x0
        Nonconvex.updateapprox!(w.dualmodel, x0)
    end
    Nonconvex.assess_convergence!(w)
    return w
end

function reset!(w::Nonconvex.PercivalWorkspace, x0 = nothing)
    w.counter[] = 0
    if x0 !== nothing
        w.x0 .= x0
    end
    return w
end

# For adaptive SIMP
function (o::Optimizer)(workspace::Nonconvex.Workspace)
    mma_results = Nonconvex.optimize!(workspace)
    return mma_results
end

@params struct MMAOptionsGen
    maxiter
    outer_maxiter
    tol
    s_init
    s_incr
    s_decr
    store_trace
    show_trace
    auto_scale
    dual_options
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::MMAOptionsGen) = println("TopOpt MMA options generator")
function (g::MMAOptionsGen)(i)
    Nonconvex.MMAOptions(
        maxiter = g.maxiter(i),
        outer_maxiter = g.outer_maxiter(i),
        tol = g.tol(i),
        s_init = g.s_init(i),
        s_incr = g.s_incr(i),
        s_decr = g.s_decr(i),
        store_trace = g.store_trace(i),
        show_trace = g.show_trace(i),
        auto_scale = g.auto_scale(i),
        dual_options = g.dual_options(i)
    )
end

function (g::MMAOptionsGen)(options, i)
    Nonconvex.MMAOptions(
        maxiter = optionalcall(g, :maxiter, options, i),
        outer_maxiter = optionalcall(g, :outer_maxiter, options, i),
        tol = optionalcall(g, :tol, options, i),
        s_init = optionalcall(g, :s_init, options, i),
        s_incr = optionalcall(g, :s_incr, options, i),
        s_decr = optionalcall(g, :s_decr, options, i),
        store_trace = optionalcall(g, :store_trace, options, i),
        show_trace = optionalcall(g, :show_trace, options, i),
        auto_scale = optionalcall(g, :auto_scale, options, i),
        dual_options = optionalcall(g, :dual_options, options, i)
    )
end
function optionalcall(g, s, options, i)
    if getproperty(g, s) isa Nothing
        return getproperty(options, s)
    else
        return getproperty(g, s)(i)
    end
end


text:
struct ROS3PConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  C21::T
  C31::T
  C32::T
  b1::T
  b2::T
  b3::T
  btilde1::T
  btilde2::T
  btilde3::T
  gamma::T2
  c2::T2
  c3::T2
  d1::T
  d2::T
  d3::T
end

function ROS3PConstantCache(T::Type,T2::Type)
  gamma = T(1/2 + sqrt(3)/6)
  igamma = inv(gamma)
  a21 = T(igamma)
  a31 = T(igamma)
  a32 = T(0)
  C21 = T(-igamma^2)
  tmp = -igamma*(2 - (1/2)*igamma)
  C31 = -igamma*(1-tmp)
  C32 = tmp
  tmp = igamma*(2/3 - (1/6)*igamma)
  b1 = igamma*(1+tmp)
  b2 = tmp
  b3 = (1/3)*igamma
  # btilde1 = T(2.113248654051871)
  # btilde2 = T(1.000000000000000)
  # btilde3 = T(0.4226497308103742)
  btilde1 = b1 - T(2.113248654051871)
  btilde2 = b2 - T(1.000000000000000)
  btilde3 = b3 - T(0.4226497308103742)
  c2 = T(1)
  c3 = T(1)
  d1 = T(0.7886751345948129)
  d2 = T(-0.2113248654051871)
  d3 = T(-1.077350269189626)
  ROS3PConstantCache(a21,a31,a32,C21,C31,C32,b1,b2,b3,btilde1,btilde2,btilde3,gamma,c2,c3,d1,d2,d3)
end

struct Rodas3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  b1::T
  b2::T
  b3::T
  b4::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  gamma::T2
  c2::T2
  c3::T2
  d1::T
  d2::T
  d3::T
  d4::T
end

function Rodas3ConstantCache(T::Type,T2::Type)
  gamma = T(1//2)
  a21 = T(0)
  a31 = T(2)
  a32 = T(0)
  a41 = T(2)
  a42 = T(0)
  a43 = T(1)
  C21 = T(4)
  C31 = T(1)
  C32 = T(-1)
  C41 = T(1)
  C42 = T(-1)
  C43 = T(-8//3)
  b1 = T(2)
  b2 = T(0)
  b3 = T(1)
  b4 = T(1)
  btilde1 = T(0.0)
  btilde2 = T(0.0)
  btilde3 = T(0.0)
  btilde4 = T(1.0)
  c2 = T(0.0)
  c3 = T(1.0)
  c4 = T(1.0)
  d1 = T(1//2)
  d2 = T(3//2)
  d3 = T(0)
  d4 = T(0)
  Rodas3ConstantCache(a21,a31,a32,a41,a42,a43,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

struct Ros4ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  b1::T
  b2::T
  b3::T
  b4::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  gamma::T2
  c2::T2
  c3::T2
  d1::T
  d2::T
  d3::T
  d4::T
end

function RosShamp4ConstantCache(T::Type,T2::Type)
  a21=T(2)
  a31=T(48//25)
  a32=T(6//25)
  C21=T(-8)
  C31=T(372//25)
  C32=T(12//5)
  C41=T(-112//125)
  C42=T(-54//125)
  C43=T(-2//5)
  b1=T(19//9)
  b2=T(1//2)
  b3=T(25//108)
  b4=T(125//108)
  btilde1=T(17//54)
  btilde2=T(7//36)
  btilde3=T(0)
  btilde4=T(125//108)
  gamma=T2(1//2)
  c2= T2(1)
  c3= T2(3//5)
  d1=T( 1//2)
  d2=T(-3//2)
  d3=T( 2.42)
  d4=T( 0.116)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

function Veldd4ConstantCache(T::Type,T2::Type)
  a21= T(2.000000000000000)
  a31= T(4.812234362695436)
  a32= T(4.578146956747842)
  C21=-T(5.333333333333331)
  C31= T(6.100529678848254)
  C32= T(1.804736797378427)
  C41=-T(2.540515456634749)
  C42=-T(9.443746328915205)
  C43=-T(1.988471753215993)
  b1= T(4.289339254654537)
  b2= T(5.036098482851414)
  b3= T(0.6085736420673917)
  b4= T(1.355958941201148)
  btilde1= T(2.175672787531755)
  btilde2= T(2.950911222575741)
  btilde3=-T(.7859744544887430)
  btilde4=-T(1.355958941201148)
  gamma= T2(0.2257081148225682)
  c2= T2(0.4514162296451364)
  c3= T2(0.8755928946018455)
  d1= T(0.2257081148225682)
  d2=-T(0.04599403502680582)
  d3= T(0.5177590504944076)
  d4=-T(0.03805623938054428)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

function Velds4ConstantCache(T::Type,T2::Type)
  a21= T(2)
  a31= T(7//4)
  a32= T(1//4)
  C21=-T(8)
  C31=-T(8)
  C32=-T(1)
  C41= T(1//2)
  C42=-T(1//2)
  C43= T(2)
  b1= T(4//3)
  b2= T(2//3)
  b3=-T(4//3)
  b4= T(4//3)
  btilde1=-T(1//3)
  btilde2=-T(1//3)
  btilde3=-T(0)
  btilde4=-T(4//3)
  gamma= T2(1//2)
  c2= T2(1)
  c3= T2(1//2)
  d1= T(1//2)
  d2=-T(3//2)
  d3=-T(3//4)
  d4= T(1//4)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

function GRK4TConstantCache(T::Type,T2::Type)
  a21= T(2)
  a31= T(4.524708207373116)
  a32= T(4.163528788597648)
  C21=-T(5.071675338776316)
  C31= T(6.020152728650786)
  C32= T(0.1597506846727117)
  C41=-T(1.856343618686113)
  C42=-T(8.505380858179826)
  C43=-T(2.084075136023187)
  b1= T(3.957503746640777)
  b2= T(4.624892388363313)
  b3= T(0.6174772638750108)
  b4= T(1.282612945269037)
  btilde1= T(2.302155402932996)
  btilde2= T(3.073634485392623)
  btilde3=-T(0.8732808018045032)
  btilde4=-T(1.282612945269037)
  gamma= T2(0.231)
  c2= T2(0.462)
  c3= T2(0.8802083333333334)
  d1= T(0.2310000000000000)
  d2=-T(0.03962966775244303)
  d3= T(0.5507789395789127)
  d4=-T(0.05535098457052764)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

function GRK4AConstantCache(T::Type,T2::Type)
  a21= T(1.108860759493671)
  a31= T(2.377085261983360)
  a32= T(0.1850114988899692)
  C21=-T(4.920188402397641)
  C31= T(1.055588686048583)
  C32= T(3.351817267668938)
  C41= T(3.846869007049313)
  C42= T(3.427109241268180)
  C43=-T(2.162408848753263)
  b1= T(1.845683240405840)
  b2= T(0.1369796894360503)
  b3= T(0.7129097783291559)
  b4= T(0.6329113924050632)
  btilde1= T(0.04831870177201765)
  btilde2=-T(0.6471108651049505)
  btilde3= T(0.2186876660500240)
  btilde4=-T(0.6329113924050632)
  gamma= T2(0.3950000000000000)
  c2= T2(0.4380000000000000)
  c3= T2(0.8700000000000000)
  d1= T(0.3950000000000000)
  d2=-T(0.3726723954840920)
  d3= T(0.06629196544571492)
  d4= T(0.4340946962568634)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

function Ros4LStabConstantCache(T::Type,T2::Type)
  a21= T(2.000000000000000)
  a31= T(1.867943637803922)
  a32= T(0.2344449711399156)
  C21=-T(7.137615036412310)
  C31= T(2.580708087951457)
  C32= T(0.6515950076447975)
  C41=-T(2.137148994382534)
  C42=-T(0.3214669691237626)
  C43=-T(0.6949742501781779)
  b1= T(2.255570073418735)
  b2= T(0.2870493262186792)
  b3= T(0.4353179431840180)
  b4= T(1.093502252409163)
  btilde1=-T(0.2815431932141155)
  btilde2=-T(0.07276199124938920)
  btilde3=-T(0.1082196201495311)
  btilde4=-T(1.093502252409163)
  gamma= T2(0.5728200000000000)
  c2= T2(1.145640000000000)
  c3= T2(0.6552168638155900)
  d1= T(0.5728200000000000)
  d2=-T(1.769193891319233)
  d3= T(0.7592633437920482)
  d4=-T(0.1049021087100450)
  Ros4ConstantCache(a21,a31,a32,C21,C31,C32,C41,C42,C43,b1,b2,b3,b4,btilde1,btilde2,btilde3,btilde4,gamma,c2,c3,d1,d2,d3,d4)
end

struct RodasConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  C51::T
  C52::T
  C53::T
  C54::T
  C61::T
  C62::T
  C63::T
  C64::T
  C65::T
  gamma::T
  c2::T2
  c3::T2
  c4::T2
  d1::T
  d2::T
  d3::T
  d4::T
  h21::T
  h22::T
  h23::T
  h24::T
  h25::T
  h31::T
  h32::T
  h33::T
  h34::T
  h35::T
end

function Rodas4ConstantCache(T::Type,T2::Type)
  gamma=T(1//4)
  #BET2P=0.0317D0
  #BET3P=0.0635D0
  #BET4P=0.3438D0
  a21= T(1.544000000000000)
  a31= T(0.9466785280815826)
  a32= T(0.2557011698983284)
  a41= T(3.314825187068521)
  a42= T(2.896124015972201)
  a43= T(0.9986419139977817)
  a51= T(1.221224509226641)
  a52= T(6.019134481288629)
  a53= T(12.53708332932087)
  a54=-T(0.6878860361058950)
  C21=-T(5.668800000000000)
  C31=-T(2.430093356833875)
  C32=-T(0.2063599157091915)
  C41=-T(0.1073529058151375)
  C42=-T(9.594562251023355)
  C43=-T(20.47028614809616)
  C51= T(7.496443313967647)
  C52=-T(10.24680431464352)
  C53=-T(33.99990352819905)
  C54= T(11.70890893206160)
  C61= T(8.083246795921522)
  C62=-T(7.981132988064893)
  C63=-T(31.52159432874371)
  C64= T(16.31930543123136)
  C65=-T(6.058818238834054)

  c2=T2(0.386)
  c3=T2(0.21)
  c4=T2(0.63)

  d1= T(0.2500000000000000)
  d2=-T(0.1043000000000000)
  d3= T(0.1035000000000000)
  d4=-T(0.03620000000000023)

  h21= T(10.12623508344586)
  h22=-T(7.487995877610167)
  h23=-T(34.80091861555747)
  h24=-T(7.992771707568823)
  h25= T(1.025137723295662)
  h31=-T(0.6762803392801253)
  h32= T(6.087714651680015)
  h33= T(16.43084320892478)
  h34= T(24.76722511418386)
  h35=-T(6.594389125716872)

  RodasConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

function Rodas42ConstantCache(T::Type,T2::Type)
  gamma= T(1//4)
  #BET2P=0.0317D0
  #BET3P=0.0047369D0
  #BET4P=0.3438D0
  a21= T(1.402888400000000)
  a31= T(0.6581212688557198)
  a32=-T(1.320936088384301)
  a41= T(7.131197445744498)
  a42= T(16.02964143958207)
  a43=-T(5.561572550509766)
  a51= T(22.73885722420363)
  a52= T(67.38147284535289)
  a53=-T(31.21877493038560)
  a54= T(0.7285641833203814)
  C21=-T(5.104353600000000)
  C31=-T(2.899967805418783)
  C32= T(4.040399359702244)
  C41=-T(32.64449927841361)
  C42=-T(99.35311008728094)
  C43= T(49.99119122405989)
  C51=-T(76.46023087151691)
  C52=-T(278.5942120829058)
  C53= T(153.9294840910643)
  C54= T(10.97101866258358)
  C61=-T(76.29701586804983)
  C62=-T(294.2795630511232)
  C63= T(162.0029695867566)
  C64= T(23.65166903095270)
  C65=-T(7.652977706771382)
  c2=T2(0.3507221)
  c3=T2(0.2557041)
  c4=T2(0.6817790)
  d1= T(0.2500000000000000)
  d2=-T(0.06902209999999998)
  d3=-T(0.0009671999999999459)
  d4=-T(0.08797900000000025)

  h21=-T(38.71940424117216)
  h22=-T(135.8025833007622)
  h23= T(64.51068857505875)
  h24=-T(4.192663174613162)
  h25=-T(2.531932050335060)
  h31=-T(14.99268484949843)
  h32=-T(76.30242396627033)
  h33= T(58.65928432851416)
  h34= T(16.61359034616402)
  h35=-T(0.6758691794084156)

  RodasConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

function Rodas4PConstantCache(T::Type,T2::Type)
  gamma = T(1//4)
	#BET2P=0.D0
	#BET3P=c3*c3*(c3/6.d0-GAMMA/2.d0)/(GAMMA*GAMMA)
  #BET4P=0.3438D0
  a21= T(3)
  a31= T(1.831036793486759)
  a32= T(0.4955183967433795)
  a41= T(2.304376582692669)
  a42=-T(0.05249275245743001)
  a43=-T(1.176798761832782)
  a51=-T(7.170454962423024)
  a52=-T(4.741636671481785)
  a53=-T(16.31002631330971)
  a54=-T(1.062004044111401)
  C21=-T(12)
  C31=-T(8.791795173947035)
  C32=-T(2.207865586973518)
  C41= T(10.81793056857153)
  C42= T(6.780270611428266)
  C43= T(19.53485944642410)
  C51= T(34.19095006749676)
  C52= T(15.49671153725963)
  C53= T(54.74760875964130)
  C54= T(14.16005392148534)
  C61= T(34.62605830930532)
  C62= T(15.30084976114473)
  C63= T(56.99955578662667)
  C64= T(18.40807009793095)
  C65=-T(5.714285714285717)
  c2=T2(3*gamma)
  c3=T2(0.21)
  c4=T2(0.63)
  d1=T( 0.2500000000000000)
  d2=T(-0.5000000000000000)
  d3=T(-0.0235040000000000)
  d4=T(-0.0362000000000000)

  h21= T(25.09876703708589)
  h22= T(11.62013104361867)
  h23= T(28.49148307714626)
  h24=-T(5.664021568594133)
  h25= T(0)
  h31= T(1.638054557396973)
  h32=-T(0.7373619806678748)
  h33= T(8.477918219238990)
  h34= T(15.99253148779520)
  h35=-T(1.882352941176471)

  RodasConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63,C64,C65,
                    gamma,c2,c3,c4,d1,d2,d3,d4,
                    h21,h22,h23,h24,h25,h31,h32,h33,h34,h35)
end

struct Rodas5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  C21::T
  C31::T
  C32::T
  C41::T
  C42::T
  C43::T
  C51::T
  C52::T
  C53::T
  C54::T
  C61::T
  C62::T
  C63::T
  C64::T
  C65::T
  C71::T
  C72::T
  C73::T
  C74::T
  C75::T
  C76::T
  C81::T
  C82::T
  C83::T
  C84::T
  C85::T
  C86::T
  C87::T
  gamma::T2
  d1::T
  d2::T
  d3::T
  d4::T
  d5::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
end

function Rodas5ConstantCache(T::Type,T2::Type)
  gamma = T2(.19)
  a21 = T(2.0)
  a31 = T(3.040894194418781 )
  a32 = T(1.041747909077569 )
  a41 = T(2.576417536461461 )
  a42 = T(1.622083060776640 )
  a43 = T(-.9089668560264532)
  a51 = T(2.760842080225597 )
  a52 = T(1.446624659844071 )
  a53 = T(-.3036980084553738)
  a54 = T(.2877498600325443 )
  a61 = T(-14.09640773051259)
  a62 = T(6.925207756232704 )
  a63 = T(-41.47510893210728)
  a64 = T(2.343771018586405 )
  a65 = T(24.13215229196062 )
  C21 = T(-10.31323885133993)
  C31 = T(-21.04823117650003)
  C32 = T(-7.234992135176716)
  C41 = T(32.22751541853323 )
  C42 = T(-4.943732386540191)
  C43 = T(19.44922031041879 )
  C51 = T(-20.69865579590063)
  C52 = T(-8.816374604402768)
  C53 = T(1.260436877740897 )
  C54 = T(-.7495647613787146)
  C61 = T(-46.22004352711257)
  C62 = T(-17.49534862857472)
  C63 = T(-289.6389582892057)
  C64 = T(93.60855400400906 )
  C65 = T(318.3822534212147 )
  C71 = T(34.20013733472935 )
  C72 = T(-14.15535402717690)
  C73 = T(57.82335640988400 )
  C74 = T(25.83362985412365 )
  C75 = T(1.408950972071624 )
  C76 = T(-6.551835421242162)
  C81 = T(42.57076742291101 )
  C82 = T(-13.80770672017997)
  C83 = T(93.98938432427124 )
  C84 = T(18.77919633714503 )
  C85 = T(-31.58359187223370)
  C86 = T(-6.685968952921985)
  C87 = T(-5.810979938412932)
  c2 = T2(0.38              )
  c3 = T2(0.3878509998321533)
  c4 = T2(0.4839718937873840)
  c5 = T2(0.4570477008819580)
  d1 = T(gamma                 )
  d2 = T(-0.1823079225333714636)
  d3 = T(-0.319231832186874912 )
  d4 = T( 0.3449828624725343   )
  d5 = T(-0.377417564392089818 )

  #=
  a71 = -14.09640773051259
  a72 = 6.925207756232704
  a73 = -41.47510893210728
  a74 = 2.343771018586405
  a75 = 24.13215229196062
  a76 = T(1)
  a81 = -14.09640773051259
  a82 = 6.925207756232704
  a83 = -41.47510893210728
  a84 = 2.343771018586405
  a85 = 24.13215229196062
  a86 = T(1)
  a87 = T(1)
  b1 = -14.09640773051259
  b2 = 6.925207756232704
  b3 = -41.47510893210728
  b4 = 2.343771018586405
  b5 = 24.13215229196062
  b6 = T(1)
  b7 = T(1)
  b8 = T(1)
  =#

  Rodas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                      a61,a62,a63,a64,a65,
                      C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,
                      C61,C62,C63,C64,C65,C71,C72,C73,C74,C75,C76,
                      C81,C82,C83,C84,C85,C86,C87,
                      gamma,d1,d2,d3,d4,d5,c2,c3,c4,c5)
end

#=
# alpha_ij
A = [0 0 0 0 0 0 0 0
     big"0.38" 0 0 0 0 0 0 0
     big"0.1899188971074152"    big"0.1979321027247381"  0 0 0 0 0 0
     big"0.1110729281178426"    big"0.5456026683145674"  big"-0.1727037026450261" 0 0 0 0 0
     big"0.2329444418850307"    big"0.025099380960713898" big"0.1443314046300300"  big"0.054672473406183418" 0 0 0 0
     big"-0.036201017843430883" big"4.208448872731939"   big"-7.549674427720996"  big"-0.2076823626400282" big"4.585108935472517" 0 0 0
     big"7.585261698003052"     big"-15.57426208319938"  big"-8.814406895608121"  big"1.534698996826085"   big"16.07870828397837" big"0.19" 0 0
     big"0.4646018839086969"    big"0"                   big"-1.720907508837576"  big"0.2910480220957973"  big"1.821778861539924" big"-0.046521258706842056" big"0.19" 0]
# bi
B = [big"0.464601884",0,big"-1.72090751",big"0.29104802",big"1.82177886",big"-0.02674488",big"-0.01977638",big"0.19"]

# Beta_ij

Beta = [big"0.19" 0 0 0 0 0 0 0
        big"0.0076920774666285364"  big"0.19" 0 0 0 0 0 0
        big"-0.058129718999580252"  big"-0.063251113355141360"  big"0.19" 0 0 0 0 0
        big"0.7075715596134048"     big"-0.5980299539145789"   big"0.5294131505610923"    big"0.19" 0 0 0 0
        big"-0.034975026573934865"  big"-0.1928476085817357"    big"0.089839586125126941"  big"0.027613185520411822" big"0.19" 0 0 0
        big"7.585261698003052"      big"-15.57426208319938"     big"-8.814406895608121"    big"1.534698996826085"    big"16.07870828397837"  big"0.19" 0 0
        big"0.4646018839086969"     0                           big"-1.720907508837576"    big"0.2910480220957973"   big"1.821778861539924"  big"-0.046521258706842056" big"0.19" 0
        big"0.4646018839086969"     0                           big"-1.720907508837576"    big"0.2910480220957973"   big"1.821778861539924"  big"-0.026744882930135193" big"-0.019776375776706864" big"0.19"]

Gamma = Beta - A
a = A*Gamma
m = B'*inv(Gamma) # = b_i
C = inv(Diagonal(diag(Gamma))) - inv(Gamma)
c = sum(A,2)
D = Beta*inv(Gamma)
d = sum(Gamma,2)

# Dense output
D_i = tanspose(D)*k_i
y1(Θ) = y₀*(1-Θ) + Θ*(y₁ + (Θ-1)*(D₂ + D₄ + (Θ + 1)*(D₃ + ΘD₄)))

# Determining coefficients
gamma = 0.19
c3 = 0.3878509998321533 == alpha3
c4 = 0.4839718937873840 == alpha4
c5 = 0.4570477008819580 == alpha5
beta3 = 6.8619167645278386e-2
beta4 = 0.8289547562599182
beta5 = 7.9630136489868164e-2
alpha64 = -0.2076823627400282
=#


text:
function loadcavelist(filename::AbstractString)
    links = Dict{SubString, Vector{SubString}}()
    function addlink(l, r)
        if haskey(links, l)
            push!(links[l], r)
        else
            links[l] = [r]
        end
    end
    for line in readlines(filename)
        l, r = split(line, '-')
        addlink(l, r)
        addlink(r, l)
    end
    links
end

isbigcave(cave::AbstractString) = first(cave) < 'a'

function part1(filename::AbstractString)
    length(takepath(loadcavelist(filename), "start"))
end

function takepath(links::Dict, here::AbstractString, history = [])
    push!(history, here)
    if here == "end"
        [history]
    else
        options = filter(links[here]) do cave
            isbigcave(cave) || !(cave in history)
        end
        vcat((takepath(links, options, copy(history)) for options in options)...)
    end
end

if basename(pwd()) == "aoc"
    cd("2021/12")
end

@assert part1("example.txt") == 10
@assert part1("example2.txt") == 19
@assert part1("example3.txt") == 226
# clipboard(@show part1("input.txt"))

function part2(filename::AbstractString)
    length(takepath2(loadcavelist(filename), "start"))
end

function takepath2(links::Dict, here::AbstractString, history = [])
    push!(history, here)
    if here == "end"
        [history]
    else
        options = filter(links[here]) do cave
            cave != "start" && (
                isbigcave(cave) ||
                !(cave in history) ||
                allunique(filter(!isbigcave, history))
            )
        end
        vcat((takepath2(links, option, copy(history)) for option in options)...)
    end
end

@assert part2("example.txt") == 36
@assert part2("example2.txt") == 103
@assert part2("example3.txt") == 3509
clipboard(@show part2("input.txt"))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
export MPS,
       position!,
       prime!,
       primelinks!,
       simlinks!,
       inner,
       randomMPS,
       maxDim,
       linkindex,
       siteindex,
       siteinds

mutable struct MPS
  N_::Int
  A_::Vector{ITensor}
  llim_::Int
  rlim_::Int

  MPS() = new(0,Vector{ITensor}(),0,0)

  function MPS(N::Int,
               A::Vector{ITensor},
               llim::Int=0,
               rlim::Int=N+1)
    new(N,A,llim,rlim)
  end

  function MPS(sites::SiteSet)
    N = length(sites)
    v = Vector{ITensor}(undef, N)
    l = [Index(1, "Link,l=$ii") for ii=1:N-1]
    @inbounds for ii in eachindex(sites)
      s = sites[ii]
      if ii == 1
        v[ii] = ITensor(l[ii], s)
      elseif ii == N
        v[ii] = ITensor(l[ii-1], s)
      else
        v[ii] = ITensor(l[ii-1], l[ii], s)
      end
    end
    new(N,v)
  end

  function MPS(is::InitState)
    N = length(is)
    As = Vector{ITensor}(undef,N)
    links  = Vector{Index}(undef,N)
    @inbounds for n in eachindex(is)
      s = sites(is)[n]
      links[n] = Index(1,"Link,l=$n")
      if n == 1
        A = ITensor(s,links[n])
        A[state(is,n),links[n](1)] = 1.0
      elseif n == N
        A = ITensor(links[n-1],s)
        A[links[n-1](1),state(is,n)] = 1.0
      else
        A = ITensor(links[n-1],links[n],s)
        A[links[n-1](1),state(is,n),links[n](1)] = 1.0
      end
      As[n] = A
    end
    new(N,As,0,2)
  end
end
MPS(N::Int, d::Int, opcode::String) = MPS(InitState(Sites(N,d), opcode))
MPS(N::Int) = MPS(N,Vector{ITensor}(undef,N),0,N+1)
MPS(s::SiteSet, opcode::String) = MPS(InitState(s, opcode))

function randomMPS(sites)
  M = MPS(sites)
  @inbounds for i ∈ eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  M.llim_ = 1
  M.rlim_ = length(M)
  return M
end

length(m::MPS) = m.N_
leftLim(m::MPS) = m.llim_
rightLim(m::MPS) = m.rlim_

getindex(m::MPS, n::Integer) = getindex(m.A_,n)
setindex!(m::MPS,T::ITensor,n::Integer) = setindex!(m.A_,T,n)

copy(m::MPS) = MPS(m.N_,copy(m.A_),m.llim_,m.rlim_)

eachindex(m::MPS) = 1:length(m)

"""
    dag(m::MPS)

Hermitian conjugation of a matrix product state `m`.
"""
function dag(m::MPS)
  N = length(m)
  mdag = MPS(N)
  @inbounds for i ∈ eachindex(m)
    mdag[i] = dag(m[i])
  end
  return mdag
end

function show(io::IO, M::MPS)
  print(io,"MPS")
  (length(M) > 0) && print(io,"\n")
  @inbounds for (i, m) ∈ enumerate(inds.(M.A_))
    println(io,"$i  $m")
  end
end

function linkindex(M::MPS,j::Integer)
  N = length(M)
  j ≥ length(M) && error("No link index to the right of site $j (length of MPS is $N)")
  li = commonindex(M[j],M[j+1])
  if isdefault(li)
    error("linkindex: no MPS link index at link $j")
  end
  return li
end

function siteindex(M::MPS,j::Integer)
  N = length(M)
  if j == 1
    si = uniqueindex(M[j],M[j+1])
  elseif j == N
    si = uniqueindex(M[j],M[j-1])
  else
    si = uniqueindex(M[j],(M[j-1],M[j+1]))
  end
  return si
end

function siteinds(M::MPS)
  N = length(M)
  is = IndexSet(N)
  @inbounds for j ∈ eachindex(M)
    is[j] = siteindex(M,j)
  end
  return is
end

function replacesites!(M::MPS,sites)
  N = length(M)
  @inbounds for j in eachindex(M)
    sj = siteindex(M,j)
    replaceindex!(M[j],sj,sites[j])
  end
  return M
end

function position!(M::MPS,
                   j::Integer)
  N = length(M)
  while leftLim(M) < (j-1)
    ll = leftLim(M)+1
    s = findindex(M[ll],"Site")
    if ll == 1
      (Q,R) = qr(M[ll],s)
    else
      li = linkindex(M,ll-1)
      (Q,R) = qr(M[ll],s,li)
    end
    M[ll] = Q
    M[ll+1] *= R
    M.llim_ += 1
  end

  while rightLim(M) > (j+1)
    rl = rightLim(M)-1
    s = findindex(M[rl],"Site")
    if rl == N
      (Q,R) = qr(M[rl],s)
    else
      ri = linkindex(M,rl)
      (Q,R) = qr(M[rl],s,ri)
    end
    M[rl] = Q
    M[rl-1] *= R
    M.rlim_ -= 1
  end
  M.llim_ = j-1
  M.rlim_ = j+1
end

"""
inner(ψ::MPS, ϕ::MPS)

Compute <ψ|ϕ>
"""
function inner(M1::MPS, M2::MPS)::Number
  N = length(M1)
  if length(M2) != N
      throw(DimensionMismatch("inner: mismatched lengths $N and $(length(M2))"))
  end
  M1dag = dag(M1)
  simlinks!(M1dag)
  O = M1dag[1]*M2[1]
  @inbounds for j ∈ eachindex(M1)[2:end]
    O *= M1dag[j]*M2[j]
  end
  return O[]
end

function replaceBond!(M::MPS,
                      b::Int,
                      phi::ITensor;
                      kwargs...)
  FU,FV = factorize(phi,inds(M[b]); which_factorization="automatic", kwargs...)
  M[b]   = FU
  M[b+1] = FV
end


text:
using LinearAlgebra

mutable struct Adam{T, TF} <: AbstractOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    α::T
    f::TF
    β1::T
    β2::T
    ϵ::T
    t::Int
    m::Array{T}
    v::Array{T}
    Δ::Array{T}

    function Adam(fun::AbstractFuncApprox, α::T; β1::T = 0.9, β2::T = 0.999, ϵ::T =1e-8) where {T<:Real}
        tp = eltype(get_params(fun))
        m = zeros(tp, size(get_params(fun)))
        v = zeros(tp, size(get_params(fun)))
        Δ = zeros(tp, size(get_params(fun)))
        new{tp, typeof(fun)}(tp(α), fun, β1, β2, ϵ, 0, m, v, Δ)
    end
end

function update!(opt::Adam{T}, g::Array{T}) where {T<:Real}
    opt.t += 1
    t, β1, β2, ϵ = opt.t, opt.β1, opt.β2, opt.ϵ
    @. opt.m *= β1
    @. opt.m += (1.0 - β1) * g
    @. opt.v *= β2
    @. opt.v += (1.0 - β2) * g^2
    α = opt.α * √(1.0 - β2^t) / (1.0 - β1^t)
    @. opt.Δ = α * opt.m / (√opt.v + ϵ)
    # add_to_params!(opt.f, @. α * opt.m / (√opt.v + ϵ))
    add_to_params!(opt.f, opt.Δ)
end

function new_episode!(opt::Adam, rng=nothing)
    nothing
end

function clone(opt::Adam)
    opt2 = Adam(opt.f, opt.α, β1=opt.β1, β2=opt.β2, ϵ=opt.ϵ)
    opt2.t = opt.t
    opt2.m .= deepcopy(opt.m)
    opt2.v .= deepcopy(opt.v)
    opt2.Δ .= deepcopy(opt.Δ)
    return opt2
end

function clone(opt::Adam, f::AbstractFuncApprox)
    opt2 = Adam(f, opt.α, β1=opt.β1, β2=opt.β2, ϵ=opt.ϵ)
    opt2.t = opt.t
    opt2.m .= deepcopy(opt.m)
    opt2.v .= deepcopy(opt.v)
    opt2.Δ .= deepcopy(opt.Δ)
    return opt2
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((4, true),),
    ((4, false), :SOCExpPSD),
    ]
insts["fast"] = [
    ((6, true),),
    ((6, false), :SOCExpPSD),
    ((8, true),),
    ]
insts["various"] = vcat(insts["fast"], [
    ((8, false), :SOCExpPSD),
    ((10, true),),
    ((10, false), :SOCExpPSD),
    ])
return (RegionOfAttrJuMP, insts)


text:
@testset "IntSetVar.jl" begin
    @testset "isbound()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntSetVar(2, 6, "x", trailer)

        @test !SeaPearl.isbound(x)

        y = SeaPearl.IntSetVar(2, 3, "x", trailer)

        SeaPearl.exclude!(y.domain, 3)
        SeaPearl.require!(y.domain, 2)

        @test SeaPearl.isbound(y)
    end

    @testset "assignedValue()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntSetVar(5, 6, "x", trailer)

        SeaPearl.require!(x.domain, 5)
        SeaPearl.require!(x.domain, 6)

        @test SeaPearl.assignedValue(x) == Set{Int}([5, 6])

        y = SeaPearl.IntSetVar(5, 8, "y", trailer)
        @test_throws AssertionError SeaPearl.assignedValue(y)
    end

    @testset "IntSetVar()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntSetVar(5, 8, "x", trailer)

        @test SeaPearl.possible_not_required_values(x.domain) == Set{Int}([5, 6, 7, 8])
    end
end


text:
using Interact, Jedi, CairoMakie, CSV, DataFrames, Dates, ColorSchemes, Statistics

Jedi.styles.default_makie!()

##
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-4]...)

# Find date
workdir = split(dir, '/')[end]
DATE = parse(Int64, split(workdir, '_')[1])
RUN_NO = parse(Int64, split(workdir, '_')[2][end])
##

# Read output
data_df = CSV.read("$dir/output/$(DATE)_r$(RUN_NO)_growth_plate.csv", DataFrame)

df = DataFrames.DataFrame()

fig = Figure()
ax1 = Axis(fig[1, 1])
ax1.xlabel = "time [min]"
ax1.ylabel = "normalized OD600"

ColorSchemes.seaborn_colorblind

for (strain, color) in zip(data_df.strain |> unique, ColorSchemes.seaborn_colorblind)
    sub_df = data_df[data_df.strain .== strain, :]
    gdf = groupby(sub_df, :time_min)
    mean_df = combine(gdf, :OD600_norm => (x -> (OD_mean=mean(x), OD_std=std(x))) => AsTable)
    insertcols!(mean_df, "strain"=>strain)
    append!(df, mean_df)
    lines!(ax1, mean_df.time_min, mean_df.OD_mean, label=strain)
    scatter!(ax1, mean_df.time_min, mean_df.OD_mean)
    errorbars!(ax1, mean_df.time_min, mean_df.OD_mean, mean_df.OD_std)
end

axislegend(ax1, position=:lt)#, merge = merge, unique = unique)
fig
save("$dir/output/$(DATE)_r$(RUN_NO)_all_curves.pdf", fig)
##
# Read summary output
sum_df = CSV.read("$dir/output/gp_analysis_summary.csv", DataFrame)
fig = Figure(resolution=(300,300))
ax = Axis(fig[1, 1])
scatter!(ax, sum_df.lambda_max, sum_df.exp_growth_rate)
ax.xlabel = "Gaussian Process Growth Rate [1/min]"
ax.ylabel = "Exponential Growth Rate [1/min]"
MIN = minimum(hcat(sum_df.lambda_max, sum_df.exp_growth_rate))
MAX = maximum(hcat(sum_df.lambda_max, sum_df.exp_growth_rate))
lines!(ax, [MIN, MAX], [MIN, MAX])
fig

##
fig = Figure(resolution=(600,300))
ax_min = Axis(fig[1, 1])
ax_supp = Axis(fig[1, 2])

ax_min.xticks = ([1,2], ["Gaussian Process", "Exponential Growth"])
ax_supp.xticks = ([1,2], ["Gaussian Process", "Exponential Growth"])

ax_min.title = "Minimal Medium"
ax_supp.title = "Minimal Medium with supplements"

blue = ColorSchemes.seaborn_colorblind[1]
orange = ColorSchemes.seaborn_colorblind[2]
for i in 1:3

    strain_min = "MG1655_$i"
    strain_supp = "MG1655_supp_$i"
    y1_min = sum_df[sum_df.strain .== strain_min, "lambda_max"]
    y2_min = sum_df[sum_df.strain .== strain_min, "exp_growth_rate"]
    y1_supp = sum_df[sum_df.strain .== strain_supp, "lambda_max"]
    y2_supp = sum_df[sum_df.strain .== strain_supp, "exp_growth_rate"]

    scatter!(ax_min, ones(length(y1_min)) .+randn(length(y1_min))./10 , y1_min, color=blue)
    scatter!(ax_min, ones(length(y2_min)) .* 2 .+randn(length(y1_min))./10, y2_min, color=orange)
    scatter!(ax_supp, ones(length(y1_supp)) .+randn(length(y1_min))./10, y1_supp, color=blue)
    scatter!(ax_supp, ones(length(y2_supp)) .* 2 .+randn(length(y1_min))./10, y2_supp, color=orange)

end
ylims!(ax_min, (0.006, 0.009))
ylims!(ax_supp, (0.006, 0.009))
fig
save("$dir/output/growth_rate_comparison.pdf", fig)


text:
@testset "Constructors" begin

  @test_throws ArgumentError NFFTPlan(zeros(1,4), (2,2))
  @test_throws ArgumentError NFFTPlan(zeros(2,4), (2,3))

  p = NFFTPlan(zeros(2,4), (2,2))
  pCopy = copy(p)

  for n in fieldnames(typeof(p))
    println(n)
    if n ∉ [:tmpVec, :forwardFFT, :backwardFFT, :blocks, :nodesInBlock]
      @test getfield(p,n) == getfield(pCopy,n)
    end
  end

  pAdjCopy = copy(adjoint(p)) # just ensure that copy does not error

  @show p
  @show adjoint(p)

  ## test range error

  x = [-0.6  0.9; 0.5  -0.5]
  @test_throws ArgumentError NFFTPlan(x, (2,2))
  x = [-0.3  0.3; 0.3  NaN]
  @test_throws Exception NFFTPlan(x, (2,2))
  # The previous test throws an ArgumentError in the single-threaded case
  # and an TaskFailedException in the multi-threaded case. Needs some rethrow

  ## test nodes!(p, tr)
  Nx = 32
  trj1 = rand(2, 1000) .- 0.5
  trj2 = rand(2, 1000) .- 0.5

  p1 = NFFT.NFFTPlan(trj1, (Nx, Nx))
  p2 = NFFT.NFFTPlan(trj2, (Nx, Nx))
  NFFT.nodes!(p2, trj1)

  for n in fieldnames(typeof(p1))
    println(n)
    if n ∉ [:tmpVec, :forwardFFT, :backwardFFT, :blocks, :nodesInBlock, :params]
      @test getfield(p1,n) == getfield(p2,n)
    end
  end

  for n in fieldnames(typeof(p1.params))
    println(n)
    @test getfield(p1.params,n) == getfield(p2.params,n)
  end

  for n in fieldnames(typeof(p1.forwardFFT))
      println(n)
      if n ∉ [:pinv, :plan]
          @test getfield(p1.forwardFFT,n) == getfield(p2.forwardFFT,n)
          @test getfield(p1.backwardFFT,n) == getfield(p2.backwardFFT,n)
      end
  end

end


text:
export step

################################################################################
#                                     step
################################################################################

"""
### function step(sys::StateSpace, t::Vector = []; solver=Tsit5(), plot=true)

Compute the system response to the step input.

##### Args

* sys: Dynamic system in state space form.
* t: Time vector. The solver will add the values in this vector to those in
  which the solution will be evaluated. If it is empty, then the simulation will
  occur up to 10x the slowest pole.

##### Returns

* y: The output of the system.
* t: Array with the time instants in which the solver evaluated the system.
* x: State vector at each time instant.

##### Remarks

If the system has more than one input, then the step simulation will be applied
to each input. In this case, the output vectors `y`, `t`, and `x` will be arrays
of the solution data. Hence, in this case, the data from the first input can be
accessed by mean of:

    y[1]
    t[1]
    x[1]

Moreover, a new plot window will be created for each simulated input.

"""

function step(sys::StateSpace, t::Vector = []; solver=Tsit5(), plot=true)
    # Check if the time vector was specified.
    if (isempty(t))
        # If not, simulate the system up to 10x the slowest pole.
        poles = pole(sys)

        # Get the slowest frequency of the system.
        slowest_freq = minimum(abs.(real(poles)))

        t_final = (slowest_freq != 0) ? 10/slowest_freq : 40

        # Do not simulate the entire interval if t_final is too large.
        (t_final > 1000) && (t_final = 1000)

        t = [0,t_final]
    end

    # Simulate the system considering each input.
    num_inputs = size(sys.B,2)

    if num_inputs == 1
        lsim(sys, (t)->(t > 0) ? 1 : 0, t, [];
             solver=solver, plot=plot, title="Step Response")
    else
        # Outputs.
        yo = []
        to = []
        uo = []

        for i = 1:num_inputs
            yi, ti, ui = lsim(sys, (t)->begin
                                  u = zeros(num_inputs)
                                  (t > 0) && (u[i] = 1)
                                  u
                              end, t, [];
                              solver=solver, plot=plot, title="Step Response - From Input #$i")
            push!(yo,yi)
            push!(to,ti)
            push!(uo,ui)
        end
        yo, to, uo
    end
end

"""
### function step(G::TransferFunction, t::Vector = []; solver=Tsit5(), plot=true)

Compute the system response to the step input.

##### Args

* G: Transfer function of the dynamic system.
* t: Time vector. The solver will add the values in this vector to those in
  which the solution will be evaluated. If it is empty, then the simulation will
  occur up to 10x the slowest pole.

##### Returns

* y: The output of the system.
* t: Array with the time instants in which the solver evaluated the system.
* x: State vector at each time instant.

"""

function step(G::TransferFunction, t::Array = []; solver=Tsit5(), plot=true)
    # Convert the transfer function to state space and call step.
    step(tf2ss(G), t; plot=plot, solver=solver)
end


text:
# Test multihypo computation assembly

using Test

using IncrementalInference

@testset "test IncrementalInference.assembleHypothesesElements! with only nullhypothesis..." begin

# n2_1 == (certainidx, allelements, activehypo, mhidx)
n2_1_gt1 = 1:2
n2_1_gt2_ = (3,3,0)
n2_1_gt3 = [(0,Int[1;]); (1,1:2); (2,Int[])]
n2_1_gt4_ = 20
n2_1 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 1, 2, ones(Bool, 2), 0.5 )

@test sum([n2_1_gt1;] - [n2_1[1];]) == 0
@test length(n2_1[2][1]) > n2_1_gt2_[1]
@test length(n2_1[2][2]) > n2_1_gt2_[2]
@test length(n2_1[2][3]) == n2_1_gt2_[3]
@test length(n2_1[2][1]) + length(n2_1[2][2]) == n2_1_gt4_
@test n2_1_gt3[1][1] == n2_1[3][1][1]
@test n2_1_gt3[2][1] == n2_1[3][2][1]
@test n2_1_gt3[3][1] == n2_1[3][3][1]
@test sum(n2_1_gt3[1][2] .- n2_1[3][1][2]) == 0
@test sum([n2_1_gt3[2][2];] .- [n2_1[3][2][2];]) == 0
@test sum(n2_1_gt3[3][2] .- n2_1[3][3][2]) == 0
@test sum(n2_1[4] .== 0) > n2_1_gt2_[1]
@test sum(n2_1[4] .== 1) > n2_1_gt2_[1]
@test sum( [1:n2_1_gt4_;][n2_1[4] .== 0] .== n2_1[2][1] ) == length(n2_1[2][1])
@test sum( [1:n2_1_gt4_;][n2_1[4] .== 1] .== n2_1[2][2] ) == length(n2_1[2][2])
@test length(n2_1[4]) == n2_1_gt4_

# n2_1 == (certainidx, allelements, activehypo, mhidx)
n2_1_gt1 = 1:2
n2_1_gt2_ = (3,3,0)
n2_1_gt3 = [(0,Int[2;]); (1,1:2); (2,Int[])]
n2_1_gt4_ = 20
n2_1 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 2, 2, ones(Bool, 2), 0.5 )

@test sum([n2_1_gt1;] - [n2_1[1];]) == 0
@test length(n2_1[2][1]) > n2_1_gt2_[1]
@test length(n2_1[2][2]) > n2_1_gt2_[2]
@test length(n2_1[2][3]) == n2_1_gt2_[3]
@test length(n2_1[2][1]) + length(n2_1[2][2]) == n2_1_gt4_
@test n2_1_gt3[1][1] == n2_1[3][1][1]
@test n2_1_gt3[2][1] == n2_1[3][2][1]
@test n2_1_gt3[3][1] == n2_1[3][3][1]
@test sum(n2_1_gt3[1][2] .- n2_1[3][1][2]) == 0
@test sum([n2_1_gt3[2][2];] .- [n2_1[3][2][2];]) == 0
@test sum(n2_1_gt3[3][2] .- n2_1[3][3][2]) == 0
@test sum(n2_1[4] .== 0) > n2_1_gt2_[1]
@test sum(n2_1[4] .== 1) > n2_1_gt2_[1]
@test sum( [1:n2_1_gt4_;][n2_1[4] .== 0] .== n2_1[2][1] ) == length(n2_1[2][1])
@test sum( [1:n2_1_gt4_;][n2_1[4] .== 1] .== n2_1[2][2] ) == length(n2_1[2][2])
@test length(n2_1[4]) == n2_1_gt4_

end

@testset "test IncrementalInference.assembleHypothesesElements! without multihypothesis..." begin

# certainidx = 1 ## ??
# sfidx=1, mhidx=0:  ah = [1;]
# sfidx=1, mhidx=1:  ah = [1;2]

# s2_1 == (certainidx, allelements, activehypo, mhidx)
s2_1_gt1 = 1:2
s2_1_gt2 = (Int[],1:20,Int[])
s2_1_gt3 = [(0,Int[1;]); (1,1:2); (2,Int[])]
s2_1_gt4 = ones(20)
s2_1 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 1, 2 )
@test sum([s2_1_gt1;] .- [s2_1[1];]) == 0
@test sum( s2_1_gt2[1] .- s2_1[2][1]) == 0
@test sum([s2_1_gt2[2];] .- [s2_1[2][2];]) == 0
@test sum( s2_1_gt2[3] .- s2_1[2][3]) == 0
@test s2_1_gt3[1][1] == s2_1[3][1][1]
@test sum(s2_1_gt3[1][2] .- s2_1[3][1][2]) == 0
@test s2_1_gt3[2][1] == s2_1[3][2][1]
@test sum([s2_1_gt3[2][2];] .- [s2_1[3][2][2];]) == 0
@test s2_1_gt3[3][1] == s2_1[3][3][1]
@test sum(s2_1_gt3[3][2] .- s2_1[3][3][2]) == 0
@test sum(s2_1_gt4 .- s2_1[4])  == 0

s2_2_gt1 = 1:2
s2_2_gt2 = (Int[],1:20,Int[])
s2_2_gt3 = [(0,Int[2;]); (1,1:2); (2,Int[])]
s2_2_gt4 = ones(20) # Int[]

s2_2 = IncrementalInference.assembleHypothesesElements!(nothing, 20, 2, 2 )

@test sum([s2_2_gt1;] .- [s2_2[1];]) == 0
@test sum([s2_2_gt2[1];] .- [s2_2[2][1];]) == 0
@test sum(s2_2_gt2[2] .- s2_2[2][2]) == 0
@test sum(s2_2_gt2[3] .- s2_2[2][3]) == 0
@test      s2_2_gt3[1][1]   ==  s2_2[3][1][1]
@test sum([s2_2_gt3[1][2];] .- [s2_2[3][1][2];]) == 0
@test s2_2_gt3[2][1] == s2_2[3][2][1]
@test sum([s2_2_gt3[2][2];] .- [s2_2[3][2][2];]) == 0
@test s2_2_gt3[3][1] == s2_2[3][3][1]
@test sum([s2_2_gt3[3][2];] .- [s2_2[3][3][2];]) == 0
@test sum(s2_2_gt4 .- s2_2[4])  == 0

end

@testset "assembleHypothesesElements! with bi-modality (certain variable)" begin

# certainidx = 1
# sfidx=1, mhidx=1:  ah = []
# sfidx=1, mhidx=2:  ah = [1;2]
# sfidx=1, mhidx=3:  ah = [1;3]
s3_1_gt1 = [1]
s3_1_gt2 = (0,3,3,40)
s3_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3])]
s3_1_gt4 = 40

s3_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 1, 3)

@test sum([s3_1_gt1;] - [s3_1[1];]) == 0
@test sum([s3_1_gt2[1];] .- [s3_1[2][1];]) == 0
@test length(s3_1[2][2]) > s3_1_gt2[2]
@test length(s3_1[2][3]) > s3_1_gt2[3]
@test length(s3_1[2][2]) + length(s3_1[2][3]) == s3_1_gt2[4]
@test s3_1_gt3[1][1] == s3_1[3][1][1]
@test s3_1_gt3[2][1] == s3_1[3][2][1]
@test s3_1_gt3[3][1] == s3_1[3][3][1]
@test sum(s3_1_gt3[1][2] .- s3_1[3][1][2]) == 0
@test sum(s3_1_gt3[2][2] .- s3_1[3][2][2]) == 0
@test sum(s3_1_gt3[3][2] .- s3_1[3][3][2]) == 0

@test sum(s3_1[4] .== 2) > s3_1_gt2[2]
@test sum(s3_1[4] .== 3) > s3_1_gt2[3]

@test sum( [1:40;][s3_1[4] .== 2] .== s3_1[2][2] ) == length(s3_1[2][2])
@test sum( [1:40;][s3_1[4] .== 3] .== s3_1[2][3] ) == length(s3_1[2][3])
@test length(s3_1[4]) == s3_1_gt4

end

@testset "assembleHypothesesElements! with bi-modality (fractional variable 1/2)" begin

# certainidx = 1
# sfidx=2, mhidx=1:  ah = [1;2]
# sfidx=2, mhidx=2:  ah = [2;3]
# sfidx=2, mhidx=3:  [2;3], 2 should take a value from 3
s3_2_gt1 = [1]
s3_2_gt2 = (0,3,3,40)
s3_2_gt3 = [(0, Int[2]); (1,Int[1;2]); (2,Int[1;2]); (3,Int[2;3])]
s3_2_gt4 = 40

s3_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 2, 3 )

@test sum(s3_2_gt1 - s3_2[1]) == 0
@test sum(s3_2_gt2[1] .- s3_2[2][2]) == 0
@test length(s3_2[2][1]) > 0.5*s3_2_gt2[2] # reuse test reference for bad-init nullhypo case
@test length(s3_2[2][2]) == 0
@test length(s3_2[2][3]) > s3_2_gt2[2]
@test length(s3_2[2][4]) > s3_2_gt2[3]
@test length(s3_2[2][1]) + length(s3_2[2][3]) + length(s3_2[2][4]) == s3_2_gt2[4]
@test s3_2_gt3[1][1] == s3_2[3][1][1]
@test s3_2_gt3[2][1] == s3_2[3][2][1]
@test s3_2_gt3[3][1] == s3_2[3][3][1]
@test sum(s3_2_gt3[1][2] .- s3_2[3][1][2]) == 0
@test sum(s3_2_gt3[2][2] .- s3_2[3][2][2]) == 0
@test sum(s3_2_gt3[3][2] .- s3_2[3][3][2]) == 0

@test sum(s3_2[4] .== 2) > s3_2_gt2[2]
@test sum(s3_2[4] .== 3) > s3_2_gt2[3]

@test sum( [1:40;][s3_2[4] .== 0] .== s3_2[2][1] ) == length(s3_2[2][1])
@test 0 == length(s3_2[2][2])
@test sum( [1:40;][s3_2[4] .== 2] .== s3_2[2][3] ) == length(s3_2[2][3])
@test sum( [1:40;][s3_2[4] .== 3] .== s3_2[2][4] ) == length(s3_2[2][4])
@test length(s3_2[4]) == s3_2_gt4

end

@testset "assembleHypothesesElements! with bi-modality (fractional variable 2/2)" begin

# certainidx = 1
# sfidx=3, mhidx=1:  ah = [1;3]
# sfidx=3, mhidx=2:  [2:3], 3 should take a value from 2
# sfidx=3, mhidx=3:  ah = [1;3]
s3_3_gt1 = [1]
s3_3_gt2 = (0,3,3,40)
s3_3_gt3 = [(0, Int[3]); (1,Int[1;3]); (2,Int[2;3]); (3,Int[1;3])]
s3_3_gt4 = 40

s3_3 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.5;0.5]), 40, 3, 3 )

@test sum(s3_3_gt1 - s3_3[1]) == 0
@test sum(s3_3_gt2[1] .- s3_3[2][2]) == 0
@test length(s3_3[2][1]) > 0.5*s3_3_gt2[2]
@test length(s3_3[2][2]) == 0
@test length(s3_3[2][3]) > s3_3_gt2[2]
@test length(s3_3[2][4]) > s3_3_gt2[3]
@test length(s3_3[2][1]) + length(s3_3[2][3]) + length(s3_3[2][4]) == s3_3_gt2[4]
@test s3_3_gt3[1][1] == s3_3[3][1][1]
@test s3_3_gt3[2][1] == s3_3[3][2][1]
@test s3_3_gt3[3][1] == s3_3[3][3][1]
@test sum(s3_3_gt3[1][2] .- s3_3[3][1][2]) == 0
@test sum(s3_3_gt3[2][2] .- s3_3[3][2][2]) == 0
@test sum(s3_3_gt3[3][2] .- s3_3[3][3][2]) == 0

@test sum(s3_3[4] .== 2) > s3_3_gt2[2]
@test sum(s3_3[4] .== 3) > s3_3_gt2[3]

@test sum( [1:40;][s3_3[4] .== 0] .== s3_3[2][1] ) == length(s3_3[2][1])
@test 0 == length(s3_3[2][2])
@test sum( [1:40;][s3_3[4] .== 2] .== s3_3[2][3] ) == length(s3_3[2][3])
@test sum( [1:40;][s3_3[4] .== 3] .== s3_3[2][4] ) == length(s3_3[2][4])
@test length(s3_3[4]) == s3_3_gt4

end

# @testset "test IncrementalInference.assembleHypothesesElements! with bi-modality backwards permutation..." begin

# certainidx = 1
# sfidx=1, mhidx=1:  ah = []
# sfidx=1, mhidx=2:  ah = [1;2]
# sfidx=1, mhidx=3:  ah = [1;3]
# s3_1_gt1 = [1]
# s3_1_gt2 = (0,3,3,20)
# s3_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3])]
# s3_1_gt4 = 20
#
# s3_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.5;0.5;0.0]), 20, 1, 3)

# @test sum(s3_1_gt1 - s3_1[1]) == 0
# @test sum(s3_1_gt2[1] .- s3_1[2][1]) == 0
# @test length(s3_1[2][2]) > s3_1_gt2[2]
# @test length(s3_1[2][3]) > s3_1_gt2[3]
# @test length(s3_1[2][2]) + length(s3_1[2][3]) == s3_1_gt2[4]
# @test s3_1_gt3[1][1] == s3_1[3][1][1]
# @test s3_1_gt3[2][1] == s3_1[3][2][1]
# @test s3_1_gt3[3][1] == s3_1[3][3][1]
# @test sum(s3_1_gt3[1][2] .- s3_1[3][1][2]) == 0
# @test sum(s3_1_gt3[2][2] .- s3_1[3][2][2]) == 0
# @test sum(s3_1_gt3[3][2] .- s3_1[3][3][2]) == 0
#
# @test sum(s3_1[4] .== 2) > s3_1_gt2[2]
# @test sum(s3_1[4] .== 3) > s3_1_gt2[3]
#
# @test sum( [1:20;][s3_1[4] .== 2] .== s3_1[2][2] ) == length(s3_1[2][2])
# @test sum( [1:20;][s3_1[4] .== 3] .== s3_1[2][3] ) == length(s3_1[2][3])
# @test length(s3_1[4]) == s3_1_gt4
#

# end

@testset "assembleHypothesesElements! with tri-modality (certain variable)" begin

N = 50
s4_1_gt1 = [1]
s4_1_gt2 = (0,3,3,3,N)
s4_1_gt3 = [(1,Int[]); (2,Int[1;2]); (3,Int[1;3]); (4,Int[1;4])]
s4_1_gt4 = N

s4_1 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 1, 4 )

@test sum(s4_1_gt1 - s4_1[1]) == 0
@test sum(s4_1_gt2[1] .- s4_1[2][1]) == 0
@test length(s4_1[2][2]) > s4_1_gt2[2]
@test length(s4_1[2][3]) > s4_1_gt2[3]
@test length(s4_1[2][4]) > s4_1_gt2[4]
@test length(s4_1[2][2]) + length(s4_1[2][3]) + length(s4_1[2][4]) == s4_1_gt2[5]

@test s4_1_gt3[1][1] == s4_1[3][1][1]
@test s4_1_gt3[2][1] == s4_1[3][2][1]
@test s4_1_gt3[3][1] == s4_1[3][3][1]
@test s4_1_gt3[4][1] == s4_1[3][4][1]
@test sum(s4_1_gt3[1][2] .- s4_1[3][1][2]) == 0
@test sum(s4_1_gt3[2][2] .- s4_1[3][2][2]) == 0
@test sum(s4_1_gt3[3][2] .- s4_1[3][3][2]) == 0
@test sum(s4_1_gt3[4][2] .- s4_1[3][4][2]) == 0

@test sum(s4_1[4] .== 2) > s4_1_gt2[2]
@test sum(s4_1[4] .== 3) > s4_1_gt2[3]
@test sum(s4_1[4] .== 4) > s4_1_gt2[4]

@test sum( [1:N;][s4_1[4] .== 2] .== s4_1[2][2] ) == length(s4_1[2][2])
@test sum( [1:N;][s4_1[4] .== 3] .== s4_1[2][3] ) == length(s4_1[2][3])
@test sum( [1:N;][s4_1[4] .== 4] .== s4_1[2][4] ) == length(s4_1[2][4])
@test length(s4_1[4]) == s4_1_gt4

end

@testset "assembleHypothesesElements! with tri-modality (fractional variable 1/3)" begin
## solve for fractional variable in trinary case

N = 70
s4_2_gt1 = [1]
s4_2_gt2 = (0,3,3,3,N)
s4_2_gt3 = [(0,Int[2]); (1,Int[1;2]); (2,Int[1;2]); (3,Int[2;3;4]); (4,Int[2;3;4])]
s4_2_gt4 = N

s4_2 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 2, 4 )

@test sum(s4_2_gt1 - s4_2[1]) == 0
@test length(s4_2[2][1]) > 0.5*s4_2_gt2[2]
@test sum(s4_2_gt2[2] .- s4_2[2][2]) == 0
@test length(s4_2[2][3]) > s4_2_gt2[2]
@test length(s4_2[2][4]) > s4_2_gt2[3]
@test length(s4_2[2][5]) > s4_2_gt2[4]
@test length(s4_2[2][1]) + length(s4_2[2][3]) + length(s4_2[2][4]) + length(s4_2[2][5]) == s4_2_gt2[5]

@test s4_2_gt3[1][1] == s4_2[3][1][1]
@test s4_2_gt3[2][1] == s4_2[3][2][1]
@test s4_2_gt3[3][1] == s4_2[3][3][1]
@test s4_2_gt3[4][1] == s4_2[3][4][1]
@test sum(s4_2_gt3[1][2] .- s4_2[3][1][2]) == 0
@test sum(s4_2_gt3[2][2] .- s4_2[3][2][2]) == 0
@test sum(s4_2_gt3[3][2] .- s4_2[3][3][2]) == 0
@test sum(s4_2_gt3[4][2] .- s4_2[3][4][2]) == 0

@test sum(s4_2[4] .== 0) > s4_2_gt2[2]
@test sum(s4_2[4] .== 2) > s4_2_gt2[2]
@test sum(s4_2[4] .== 3) > s4_2_gt2[3]
@test sum(s4_2[4] .== 4) > s4_2_gt2[4]

@test sum( [1:N;][s4_2[4] .== 0] .== s4_2[2][1] ) == length(s4_2[2][1])
@test 0 == length(s4_2[2][2])
@test sum( [1:N;][s4_2[4] .== 2] .== s4_2[2][3] ) == length(s4_2[2][3])
@test sum( [1:N;][s4_2[4] .== 3] .== s4_2[2][4] ) == length(s4_2[2][4])
@test sum( [1:N;][s4_2[4] .== 4] .== s4_2[2][5] ) == length(s4_2[2][5])
@test length(s4_2[4]) == s4_2_gt4

end

@testset "assembleHypothesesElements! with tri-modality (fractional variable 2/3)" begin

N = 70
s4_3_gt1 = [1]
s4_3_gt2 = (0,3,3,3,N)
s4_3_gt3 = [(0,Int[3]); (1,Int[1;3]); (2,Int[2;3;4]); (3,Int[1;3]); (4,Int[2;3;4])]
s4_3_gt4 = N

s4_3 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 3, 4 )

@test sum(s4_3_gt1 - s4_3[1]) == 0
@test length(s4_3[2][1]) > 0.5*s4_3_gt2[2]
@test sum(s4_3_gt2[2] .- s4_3[2][2]) == 0
@test length(s4_3[2][3]) > s4_3_gt2[2]
@test length(s4_3[2][4]) > s4_3_gt2[3]
@test length(s4_3[2][5]) > s4_3_gt2[4]
@test length(s4_3[2][1]) + length(s4_3[2][3]) + length(s4_3[2][4]) + length(s4_3[2][5]) == s4_3_gt2[5]

@test s4_3_gt3[1][1] == s4_3[3][1][1]
@test s4_3_gt3[2][1] == s4_3[3][2][1]
@test s4_3_gt3[3][1] == s4_3[3][3][1]
@test s4_3_gt3[4][1] == s4_3[3][4][1]
@test sum(s4_3_gt3[1][2] .- s4_3[3][1][2]) == 0
@test sum(s4_3_gt3[2][2] .- s4_3[3][2][2]) == 0
@test sum(s4_3_gt3[3][2] .- s4_3[3][3][2]) == 0
@test sum(s4_3_gt3[4][2] .- s4_3[3][4][2]) == 0

@test sum(s4_3[4] .== 0) > s4_3_gt2[2]
@test sum(s4_3[4] .== 2) > s4_3_gt2[2]
@test sum(s4_3[4] .== 3) > s4_3_gt2[3]
@test sum(s4_3[4] .== 4) > s4_3_gt2[4]

@test sum( [1:N;][s4_3[4] .== 0] .== s4_3[2][1] ) == length(s4_3[2][1])
@test 0 == length(s4_3[2][2])
@test sum( [1:N;][s4_3[4] .== 2] .== s4_3[2][3] ) == length(s4_3[2][3])
@test sum( [1:N;][s4_3[4] .== 3] .== s4_3[2][4] ) == length(s4_3[2][4])
@test sum( [1:N;][s4_3[4] .== 4] .== s4_3[2][5] ) == length(s4_3[2][5])
@test length(s4_3[4]) == s4_3_gt4

end

@testset "assembleHypothesesElements! with tri-modality (fractional variable 3/3)" begin

N = 70
s4_4_gt1 = [1]
s4_4_gt2 = (0,3,3,3,N)
s4_4_gt3 = [(0,Int[4]); (1,Int[1;4]); (2,Int[2;3;4]); (3,Int[2;3;4]); (4,Int[1;4])]
s4_4_gt4 = N

s4_4 = IncrementalInference.assembleHypothesesElements!(Categorical([0.0;0.33;0.33;0.34]), N, 4, 4 )

@test sum(s4_4_gt1 - s4_4[1]) == 0
@test length(s4_4[2][1]) > 0.5*s4_4_gt2[2]
@test sum(s4_4_gt2[2] .- s4_4[2][2]) == 0
@test length(s4_4[2][3]) > s4_4_gt2[2]
@test length(s4_4[2][4]) > s4_4_gt2[3]
@test length(s4_4[2][5]) > s4_4_gt2[4]
@test length(s4_4[2][1]) + length(s4_4[2][3]) + length(s4_4[2][4]) + length(s4_4[2][5]) == s4_4_gt2[5]

@test s4_4_gt3[1][1] == s4_4[3][1][1]
@test s4_4_gt3[2][1] == s4_4[3][2][1]
@test s4_4_gt3[3][1] == s4_4[3][3][1]
@test s4_4_gt3[4][1] == s4_4[3][4][1]
@test sum(s4_4_gt3[1][2] .- s4_4[3][1][2]) == 0
@test sum(s4_4_gt3[2][2] .- s4_4[3][2][2]) == 0
@test sum(s4_4_gt3[3][2] .- s4_4[3][3][2]) == 0
@test sum(s4_4_gt3[4][2] .- s4_4[3][4][2]) == 0

@test sum(s4_4[4] .== 0) > s4_4_gt2[2]
@test sum(s4_4[4] .== 2) > s4_4_gt2[2]
@test sum(s4_4[4] .== 3) > s4_4_gt2[3]
@test sum(s4_4[4] .== 4) > s4_4_gt2[4]

@test sum( [1:N;][s4_4[4] .== 0] .== s4_4[2][1] ) == length(s4_4[2][1])
@test 0 == length(s4_4[2][2])
@test sum( [1:N;][s4_4[4] .== 2] .== s4_4[2][3] ) == length(s4_4[2][3])
@test sum( [1:N;][s4_4[4] .== 3] .== s4_4[2][4] ) == length(s4_4[2][4])
@test sum( [1:N;][s4_4[4] .== 4] .== s4_4[2][5] ) == length(s4_4[2][5])
@test length(s4_4[4]) == s4_4_gt4

@warn "only partially testing tri-modality"

end


text:
using Subspaces
#using Convex, SCS
using LinearAlgebra
using Random
using Test

eye(n) = Matrix(1.0*I, (n,n))

@testset "Basic tests" begin
    a = random_subspace(ComplexF64, 1, 10)
    b = random_subspace(ComplexF64, 2, 10)
    c = random_subspace(ComplexF64, 3, 10)
    @test dim(a) == 1
    @test dim(b) == 2
    @test dim(c) == 3
    @test dim(a+b) == 3
    @test dim(a+c) == 4
    @test dim(b+c) == 5
    @test dim(b*c') == 6
    @test size(b*c') == (10, 10)
    @test dim(a+b+c) == 6
    @test dim((a+b) + (b+c)) == 6
    @test (a+b) + (b+c) == a+b+c
    @test (a+b) / a == Subspace([ projection(~a, x) for x in each_basis_element(b) ])
    @test (a+b) / b == Subspace([ projection(~b, x) for x in each_basis_element(a) ])
    @test (c+b) / b == Subspace([ projection(~b, x) for x in each_basis_element(c) ])
    @test (c+b) / b ⟂ b
    @test a ⊆ a
    @test a ⊇ a
    @test a ⊆ a+b
    @test a+b ⊇ a
    @test !(a ⊆ b)
    @test !(a ⊇ b)
    @test random_element(c) in c
    @test random_element(c) != random_element(c)
    x = random_element(c)
    @test projection(c, x) ≈ x
    @test frombasis(c, tobasis(c, x)) ≈ x
    @test dim(kron(b, c)) == 6
    @test size(kron(b, c)) == (100,)
    @test dim(hcat(b, c)) == 5
    @test size(hcat(b, c)) == (10, 2)
    @test dim(vcat(b, c)) == 5
    @test size(vcat(b, c)) == (20,)
    @test size(b') == (1, 10)
    @test size(b'') == (10, 1)
    @test dim(b + ones(10)) == 3
    @test dim(empty_subspace(Float64, (10,))) == 0
    @test dim(full_subspace(Float64, (10,))) == 10
    @test dim(full_subspace(Float64, (3, 4))) == 12
    @test size(full_subspace(Float64, (3, 4))) == (3, 4)
    @test dim(full_subspace(Float64, (4, 4)) / I) == 15
end

@testset "Empty spaces" begin
    a = empty_subspace(Float64, (3,4))
    b = empty_subspace(Float64, (4,5))
    f = full_subspace(Float64, (3,4))
    @test dim(a) == 0
    @test dim(~a) == 12
    @test dim(f) == 12
    @test dim(~f) == 0
    @test size(a') == (4,3)
    @test size(a+a) == (3,4)
    @test size(a*b) == (3,5)
    @test size(a&a) == (3,4)
    @test size(vcat(a)) == (3,4)
    @test size(vcat(a,a)) == (6,4)
    @test size(hcat(a)) == (3,4)
    @test size(hcat(a,a)) == (3,8)
    @test size(kron(a,b)) == (12,20)
    @test size(a/a) == (3,4)
    @test dim(a/a) == 0
    @test dim(f/a) == 12
    @test dim(f/f) == 0
    @test a ⊆ a
    @test a ⊆ f
    @test zeros((3,4)) in a
    @test !(ones((3,4)) in a)
    @test dim(random_subspace(ComplexF64, 0, (3,3))) == 0
end


text:
module Simple

using MLIR

test0 = () -> begin
    println("---- TEST 0 ----\n")

    # Constructors.
    ctx = MLIR.IR.Context()
    println(ctx)
    loc = MLIR.IR.Location(ctx)
    println(loc)
    mod = MLIR.IR.Module(loc)
    println(mod)
    op_state = MLIR.IR.OperationState("foo", loc)
    println(op_state)
    op = MLIR.IR.Operation(op_state)
    println(op)
    reg = MLIR.IR.Region()
    println(reg)
    t = MLIR.IR.Type(ctx, "index")
    println(t)
    blk = MLIR.IR.Block(t)
    arg = MLIR.IR.get_arg(blk, 0)
    println(arg)
    attr = MLIR.IR.Attribute(ctx, "\"add\"")
    println(attr)
    ident = MLIR.IR.Identifier(ctx, "type")
    println(ident)
    named_attr = MLIR.IR.NamedAttribute(ident, attr)
    println(named_attr)
end

test1 = () -> begin
    println("\n---- TEST 1 ----\n")

    # Create and destroy.
    ctx = MLIR.IR.create_context()
    MLIR.IR.num_loaded_dialects(ctx) |> y -> println("Num loaded dialects: $y")
    MLIR.IR.num_registered_dialects(ctx) |> y -> println("Num registered dialects: $y")
    MLIR.IR.destroy!(ctx)
end

test2 = () -> begin
    println("\n---- TEST 2 ----\n")

    # Create and register standard.
    ctx = MLIR.IR.create_context()
    MLIR.IR.register_standard_dialect!(ctx)
    MLIR.IR.load_standard_dialect!(ctx)
    MLIR.IR.num_loaded_dialects(ctx) |> y -> println("Num loaded dialects: $y")
    MLIR.IR.num_registered_dialects(ctx) |> y -> println("Num registered dialects: $y")
    MLIR.IR.destroy!(ctx)
end

test3 = () -> begin
    println("\n---- TEST 3 ----\n")

    # Create and dump an operation.
    ctx = MLIR.IR.create_context()
    loc = MLIR.IR.create_unknown_location(ctx)
    st = MLIR.IR.OperationState("std.add", loc)
    index_type = MLIR.IR.parse_type(ctx, "index")
    MLIR.IR.push!(st, index_type)
    op = MLIR.IR.Operation(st)
    MLIR.IR.dump(op)
end

test4 = () -> begin
    println("\n---- TEST 4 ----\n")

    # Create an operation and verify.
    ctx = MLIR.IR.create_context()
    loc = MLIR.IR.create_unknown_location(ctx)
    func_state = MLIR.IR.OperationState("func", loc)
    func_region = MLIR.IR.create_region()
    sym_name_ref = MLIR.IR.Identifier(ctx, "sym_name")
    func_name_attr = MLIR.IR.NamedAttribute(sym_name_ref, MLIR.IR.Attribute(ctx, "\"add\""))
    MLIR.IR.push!(func_state, func_region)
    MLIR.IR.push!(func_state, func_name_attr)
    type_ref = MLIR.IR.Identifier(ctx, "type")
    func_type_attr = MLIR.IR.Attribute(ctx, "(f32, f32) -> f32")
    named_func_type_attr = MLIR.IR.NamedAttribute(type_ref, func_type_attr)
    MLIR.IR.push!(func_state, named_func_type_attr)
    func = MLIR.IR.Operation(func_state)
    MLIR.IR.verify(func)
    MLIR.IR.dump(func)
end

test5 = () -> begin
    println("\n---- TEST 5 ----\n")

    # Create a more complex operation and verify.
    ctx = MLIR.IR.create_context()
    loc = MLIR.IR.create_unknown_location(ctx)
    module_op = MLIR.IR.Module(loc)
    module_body = MLIR.IR.get_body(module_op)
    memref_type = MLIR.IR.parse_type(ctx, "memref<?xf32>")
    func_body_arg_types = [memref_type, memref_type]
    func_region = MLIR.IR.create_region()
    func_body = MLIR.IR.create_block(func_body_arg_types)
    MLIR.IR.push!(func_region, func_body)
    func_type_attr = MLIR.IR.Attribute(ctx, "(memref<?xf32>, memref<?xf32>) -> ()")
    func_name_attr = MLIR.IR.Attribute(ctx, "\"add\"")
    type_ref = MLIR.IR.Identifier(ctx, "type")
    sym_name_ref = MLIR.IR.Identifier(ctx, "sym_name")
    func_attrs = [MLIR.IR.NamedAttribute(type_ref, func_type_attr), MLIR.IR.NamedAttribute(sym_name_ref, func_name_attr)]
    func_state = MLIR.IR.OperationState("func", loc)
    MLIR.IR.push!(func_state, func_attrs)
    MLIR.IR.push!(func_state, func_region)
    func = MLIR.IR.Operation(func_state)
    MLIR.IR.verify(func)
    MLIR.IR.dump(func)
end

test6 = () -> begin
    println("\n---- TEST 6 ----\n")

    # Do syntax.
    loc = MLIR.IR.Context() do ctx
        loc = MLIR.IR.Location(ctx)
        loc
    end
    println(loc)
end

test0()
test1()
test2()
test3()
test4()
test5()
test6()

end # module


text:
println("Add one to typemax:")
for t in S
    over = typemax(t) + one(t)
    @printf("%8s →  %-25s (%s)\n", t, over, typeof(over))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# included in module Wut

struct DateEncoder  <: AbstractEncoder

    season::Number
    dayOfWeek::Number
    weekend::Number
    holiday::Number
    timeOfDay::Number
    customDays::Number
    name::AbstractString
    forced::Bool

    w::Number            # cumulative width of the discontiguous 1 bits
    n::Number            # width of composite encoding

    sE::Nullable{ScalarEncoder} # SeasonEncoder
    sEo::Int64

    dE::Nullable{ScalarEncoder} # dayOfWeekEncoder
    dEo::Int64

    wE::Nullable{ScalarEncoder} # weekendEncoder
    wEo::Int64

    hE::Nullable{ScalarEncoder} # holidayEncoder
    hEo::Int64

    tE::Nullable{ScalarEncoder} # timeOfDayEncoder
    tEo::Int64

    cE::Nullable{ScalarEncoder} # customDaysEncoder
    cEo::Int64

    function DateEncoder(;season=0, dayOfWeek=0, weekend=0, holiday=0, timeOfDay=0, customDays=0, name = "", forced=true)
        w=0
        n=0
        z=Nullable{ScalarEncoder}()

        sE=Nullable{ScalarEncoder}()
        sE = season== 0 ? z : ScalarEncoder(w = season, minval=0, maxval=366,
                                         radius=91.5, periodic=true,
                                         name="season", forced=forced)
        n=isnull(sE) ? 0 : getWidth(sE)
        sEo=0

        dE=Nullable{ScalarEncoder}()
        dE = dayOfWeek== 0 ? z : ScalarEncoder(w = dayOfWeek, minval=0, maxval=7,
                                         radius=1, periodic=true,
                                         name="day of week", forced=forced)
        dEo=n
        n=isnull(dE) ? n : n + getWidth(dE)

        wE=Nullable{ScalarEncoder}()
        wE = weekend== 0 ? z : ScalarEncoder(w = weekend, minval=0, maxval=2,
                                         radius=1, periodic=true,
                                         name="weekend", forced=forced)
        wEo=n
        n=isnull(wE) ? n : n + getWidth(wE)

        cE=Nullable{ScalarEncoder}()
        cE = customDays== 0 ? z : ScalarEncoder(w = customDays, minval=0, maxval=7,
                                         radius=1, periodic=true,
                                         name="custom days", forced=forced)
        cEo=n
        n=isnull(cE) ? n : n + getWidth(cE)

        hE=Nullable{ScalarEncoder}()
        hE = holiday== 0 ? z : ScalarEncoder(w = holiday, minval=0, maxval=7,
                                         radius=1, periodic=true,
                                         name="holiday", forced=forced)
        hEo=n
        n=isnull(hE) ? n : n + getWidth(hE)

        tE=Nullable{ScalarEncoder}()
        tE = timeOfDay== 0 ? z : ScalarEncoder(w = timeOfDay, minval=0, maxval=7,
                                         radius=1, periodic=true,
                                         name="time of day", forced=forced)
        tEo=n
        n=isnull(tE) ? n : n + getWidth(tE)

            new(season, dayOfWeek, weekend, holiday, timeOfDay, customDays, name, forced,
                 w, n, sE, sEo, dE, dEo, wE, wEo, hE, hEo, tE, tEo, cE, cEo)

    end

end
export DateEncoder

function getWidth(e::DateEncoder)
    e.n
end
export getWidth

function getDescription(e::DateEncoder)
    @sprintf("%s%s%s%s%s%s",
        isnull(e.sE)?"":@sprintf(" %s %i","season",e.sEo),
        isnull(e.dE)?"":@sprintf(" %s %i","day of week",e.dEo),
        isnull(e.wE)?"":@sprintf(" %s %i","weekend",e.wEo),
        isnull(e.cE)?"":@sprintf(" %s %i","custom days",e.cEo),
        isnull(e.hE)?"":@sprintf(" %s %i","holiday",e.hEo),
        isnull(e.tE)?"":@sprintf(" %s %i","time of day",e.tEo))
end
export getDescription

function getName(e::DateEncoder)
    e.name
end
export getName

function getBucketIndices(e::DateEncoder)

end
export getBucketIndices

function encodeIntoArray(e::DateEncoder,d::DateTime,b::BitPat;learn=true, offset=0, length=0)
    fill!(b.b,false)

    timeOfDay=convert(Float32,Dates.value(Dates.Hour(d)))

    if ! isnull(e.sE) # seasson
        dy=Dates.dayofyear(d)
        b.b[e.sEo+1]=true
        encodeIntoArray(get(e.sE),dy-1,b,learn = true, offset = e.sEo)
    end
    if ! isnull(e.dE) # day of week
        dw=convert(Float32,Dates.dayofweek(d))+(timeOfDay/24)
        b.b[e.dEo+1]=true
        encodeIntoArray(get(e.dE),dw-1,b,learn = true, offset = e.dEo)
    end
    if ! isnull(e.wE) # weekend
        weekend = 0
        b.b[e.wEo+1]=true
    end
    if ! isnull(e.hE) # holiday
        isholiday = 0
        b.b[e.hEo+1]=true
    end
    if ! isnull(e.tE) # time of day
        b.b[e.tEo+1]=true
        encodeIntoArray(get(e.tE),timeOfDay,b,learn = true, offset = e.tEo)
    end
    if ! isnull(e.cE) # custom days
        customDay = 0
        b.b[e.cEo+1]=true
    end

    b
end
export encodeIntoArray

#function decode(e::DateEncoder, b::BitPat; parentFieldName="")
#
#end
#export decode

function getBucketValues(e::DateEncoder)

end
export getBucketValues

function getBucketInfo(e::DateEncoder)

end
export getBucketInfo

function topDownCompute(e::DateEncoder)

end
export topDownCompute

function closenessScores(e::DateEncoder)

end
export closenessScores

# end # module


text:
"Implements Zeckendorf's theorem and outputs a Zeckendorf sequence {cᵪ | c ∈ {0, 1}}"
function zeckendorf{T <: Integer}(n::T)
    @assert n >= 0 "n is required to be larger or equal than zero"

    x, z = 2, 2
    fib_cache = [1, 1]

    # find fib(z-1) <= n <= fib(z)
    while z <= n
        x += 1
        z = fib_cache[x-1] + fib_cache[x-2]
        push!(fib_cache, z)
    end

    c = falses(x-1)

    for i in (x-1):-1:1
        if fib_cache[i] <= n
            n -= fib_cache[i]
            c[i] = true
        end
    end

    return c
end

"Generates a binary string from a boolean array in reverse order"
binary_string{T <: AbstractArray{Bool}}(A::T) = join(map(Int, reverse(A)))

max_stars_repo_name: mooxmirror/julia-experiments

max_issues_repo_name: mooxmirror/julia-experiments

max_forks_repo_name: mooxmirror/julia-experiments