text:
import Statistics.cor

# Autocorrelation

#TODO: move to another file to maintain a more logical structure

function lagindex(x::Vector{Float64}, n::Int=1)::UnitRange{Int}
    @assert n < size(x,1)
    N = size(x, 1)
    idx = 1:N
    if n == 0
        return idx
    elseif n > 0
        return idx[n+1:N]
    elseif n < 0
        return idx[1:N+n]
    end
end

function lag(x::Vector{Float64}, n::Int=1; pad::Bool=true, padval::Float64=NaN)
    if n==0
        return x
    elseif n>0
        if pad
            return [ones(n)*padval; x[n+1:end]]
        else
            return x[n+1:end]
        end
    elseif n<0
        n *= -1
        if pad
            return [x[1:end-n]; ones(n)*padval]
        else
            return x[1:end-n]
        end
    end
end

function corlag(x::AbstractArray{T,1}, n::Int=1) where {T<:Number}
    if n == 0
        return 1.0
    end
    @assert n > 0
    @assert n < size(x,1) - 2
    idx = collect(lagindex(x, n))
    return cor(x[idx.-n], x[idx])
end

function acf(x::Vector{T}, maxlag::Int=15; lags::AbstractArray{Int,1}=0:maxlag)::Vector{Float64} where {T<:Number}
    @assert all(lags .< size(x,1)-2)
    return map((n) -> corlag(x, n), lags)
end

function acf(x::Matrix{T}, maxlag::Int=15; lags::AbstractVector{Int}=0:maxlag)::Matrix{Float64} where {T<:Number}
    k = size(x,2)
    out = zeros(Float64, (length(lags), k))
    @inbounds for j in 1:k
        out[:,j] = acf(x[:,j], maxlag, lags=lags)
    end
    return out
end

"""
    acf(x::TS,maxlag::Int=15;lags::AbstractArray{Int,1}=0:maxlag)
    acf(x::Vector{T},maxlag::Int=15;lags::AbstractArray{Int,1}=0:maxlag)::Vector{Float64}where{T<:Number}
    acf(x::Matrix{T},maxlag::Int=15;lags::AbstractVector{Int}=0:maxlag)::Matrix{Float64}where{T<:Number}

Compute the autocorrelation function of a time series.

The output will be a matrix with the same number of columns as the input time series `x`
and number of rows equal to the number of lags used as inputs to the autocorrelation function.

...
# Arguments
- `x::TS`: Time series object array containing columns on which to compute autocorrelation.
- `maxlag::Int=15`: Maximum lag of the autocorrelation series.
Optional args:
- `lags::AbstractArray{Int,1}=0:maxlag`: Explicitly specified list of lags to use (overrides use of `maxlag`).
...
"""
acf(x::TS, maxlag::Int=15; lags::AbstractArray{Int,1}=0:maxlag) = acf(x.values, maxlag; lags=lags)


text:
using DistroHelper
using Test

@testset "slug" begin
	@test DistroHelper.slug("./data/JSON/Package.toml", "./data/JSON/Versions.toml", "0.20.0") == "ebvl3"
	@test DistroHelper.slug("./data/JSON/Package.toml", "./data/JSON/Versions.toml") == "ebvl3"
end

ext: jl


text:
using Test
using Bifurcation

### Tests

function test_decimal_round()
  @test Bifurcation.decimal_round(0.001422, convert(Unsigned, 3)) == 0.001
  @test Bifurcation.decimal_round(0.000037, convert(Unsigned, 5)) == 0.00004
  @test Bifurcation.decimal_round(0.895, convert(Unsigned, 0)) == 1.0
  println("all decimal_round tests passed!")
end

function test_logistic_map()
  @test Bifurcation.logistic_map(.4, 2.6) == 0.624
  @test Bifurcation.logistic_map(.3, 1.0) == 0.21
  @test Bifurcation.decimal_round(
    Bifurcation.logistic_map(.2, 1.1),
    convert(Unsigned, 3)
    ) == 0.176
  @test Bifurcation.decimal_round(
    Bifurcation.logistic_map(.374, .7),
    convert(Unsigned, 4)
    ) == 0.1639
  println("all logistic_map tests passed!")
end

function test_logistic_map_series()
  s = Bifurcation.logistic_map_series(1, 2, 100, 100)
  @test size(s, 1) == 100

  # fancy broadcast syntax
  @test Bifurcation.decimal_round.(
    Bifurcation.logistic_map_series(0.77, 1.33, 100, 10),
    convert(Unsigned, 3)
   ) == fill(0.248, (10))
  println("all logistic_map_series tests passed!")
end

function run_tests()
  test_decimal_round();
  test_logistic_map();
  test_logistic_map_series();
end

run_tests()


text:
# Convenient methods for adding objectives to a MixedMOP
# Intended for use by the enduser.
# Requires concrete types `MixedMOP` and `VectorObjectiveFunction`.
# Internally refers to `_add_objective!`

# adding objectives to a mixed problem
@doc "Add a scalar objective to `mop::MixedMOP` modelled according to `model_config`."
function add_objective!(mop :: MixedMOP, func :: Function,
        model_config :: SurrogateConfig; can_batch = false, out_type = nothing )
    return _add_objective!( mop, VectorObjectiveFunction, func, model_config; n_out = 1, can_batch, out_type )
end

# FUNCTIONS TO ADD VECTOR OBJECTIVES TO A MixedMOP
@doc "Add a vector objective to `mop::MixedMOP` modelled according to `model_config`."
function add_vector_objective!(mop :: MixedMOP, func :: Function,
        model_config :: SurrogateConfig; n_out :: Int64, can_batch = false, out_type = nothing )
    return _add_objective!( mop, VectorObjectiveFunction, func, model_config; n_out, can_batch, out_type )
end

# OLD FUNCTIONS TO ADD OBJECTIVES TO AN MOP
# # these functions are less customizable and are mainly kept for convenience
@doc """
    add_objective!( mop :: MixedMOP, func :: T where{T <: Function}, type :: Symbol = :expensive, n_out :: Int64 = 1, can_batch :: Bool = false )

Add scalar-valued objective function `func` to `mop` structure.
`func` must take an `Vec` as its (first) argument, i.e. represent a function ``f: ℝ^n → ℝ``.
`type` must either be `:expensive` or `:cheap` to determine whether the function is replaced by a surrogate model or not.

If `type` is `:cheap` and `func` takes 1 argument only then its gradient is calculated by ForwardDiff.
A cheap function `func` with custom gradient function `grad` (representing ``∇f : ℝ^n → ℝ^n``) is added by

    add_objective!(mop, func, grad)

The optional argument `n_out` allows for the specification of vector-valued objective functions.
This is mainly meant to be used for *expensive* functions that are in some sense inter-dependent.

The flag `can_batch` defaults to false so that the objective function is simply looped over a bunch of arguments if required.
If `can_batch == true` then the objective function must be able to return an array of results when provided an array of input vectors
(whilst still returning a single result, not a singleton array containing the result, for a single input vector).

# Examples
```
# Define 2 scalar objective functions and a MOP ℝ^2 → ℝ^2

f1(x) =  x[1]^2 + x[2]

f2(x) = exp(sum(x))
∇f2(x) = exp(sum(x)) .* ones(2);

mop = MixedMOP()
add_objective!(mop, f1, :cheap)     # gradient will be calculated using ForwardDiff
add_objective!(mop, f2, ∇f2 )       # gradient is provided
```
"""
function add_objective!( mop :: MixedMOP, func :: Function, type :: Symbol = :expensive;
        n_out :: Int64 = 1, can_batch :: Bool = false, out_type = nothing)
    if type == :expensive
        objf_config = RbfConfig();
    elseif type == :cheap
        objf_config = ExactConfig(
            gradients = :autodiff
        )
    else
        error( "Type must be either `:expensive` or `:cheap`.")
    end

    if n_out > 1
    	add_vector_objective!( mop, func, objf_config; n_out, can_batch, out_type )
    else
        add_objective!(mop, func, objf_config; can_batch, out_type)
    end
end

@doc """
    add_objective!( mop :: MixedMOP, func :: T where{T <: Function}, grad :: T where{T <: Function})

Add scalar-valued objective function `func` and its vector-valued gradient `grad` to `mop` struture.
"""
function add_objective!( mop :: MixedMOP, func :: Function, grad :: Function; out_type = nothing)
    objf_config = ExactConfig(
        gradients = [grad,]
    )
    add_objective!(mop, func, objf_config; out_type)
end


text:
doc"""
    Gamma(α,θ)

The *Gamma distribution* with shape parameter `α` and scale `θ` has probability density
function

$f(x; \alpha, \beta) = \frac{x^{\alpha-1} e^{-x/\beta}}{\Gamma(\alpha) \beta^\alpha},
\quad x > 0$

```julia
Gamma()          # Gamma distribution with unit shape and unit scale, i.e. Gamma(1.0, 1.0)
Gamma(a)         # Gamma distribution with shape a and unit scale, i.e. Gamma(a, 1.0)
Gamma(a, b)      # Gamma distribution with shape a and scale b

params(d)        # Get the parameters, i.e. (a, b)
shape(d)         # Get the shape parameter, i.e. a
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Gamma distribution on Wikipedia](http://en.wikipedia.org/wiki/Gamma_distribution)

"""
immutable Gamma <: ContinuousUnivariateDistribution
    α::Float64
    θ::Float64

    function Gamma(α::Real, θ::Real)
        @check_args(Gamma, α > zero(α) && θ > zero(θ))
        new(α, θ)
    end
    function Gamma(α::Real)
        @check_args(Gamma, α > zero(α))
        new(α, 1.0)
    end
    Gamma() = new(1.0, 1.0)
end

@distr_support Gamma 0.0 Inf

#### Parameters

shape(d::Gamma) = d.α
scale(d::Gamma) = d.θ
rate(d::Gamma) = 1.0 / d.θ

params(d::Gamma) = (d.α, d.θ)

#### Statistics

mean(d::Gamma) = d.α * d.θ

var(d::Gamma) = d.α * d.θ^2

skewness(d::Gamma) = 2.0 / sqrt(d.α)

kurtosis(d::Gamma) = 6.0 / d.α

function mode(d::Gamma)
    (α, θ) = params(d)
    α >= 1.0 ? θ * (α - 1.0) : error("Gamma has no mode when shape < 1.0")
end

function entropy(d::Gamma)
    (α, θ) = params(d)
    α + lgamma(α) + (1.0 - α) * digamma(α) + log(θ)
end

mgf(d::Gamma, t::Real) = (1.0 - t * d.θ)^(-d.α)

cf(d::Gamma, t::Real) = (1.0 - im * t * d.θ)^(-d.α)

#### Evaluation & Sampling

@_delegate_statsfuns Gamma gamma α θ

gradlogpdf(d::Gamma, x::Float64) =
    insupport(Gamma, x) ? (d.α - 1.0) / x - 1.0 / d.θ : 0.0

rand(d::Gamma) = StatsFuns.Rmath.gammarand(d.α, d.θ)

#### Fit model

immutable GammaStats <: SufficientStats
    sx::Float64      # (weighted) sum of x
    slogx::Float64   # (weighted) sum of log(x)
    tw::Float64      # total sample weight

    GammaStats(sx::Real, slogx::Real, tw::Real) = new(sx, slogx, tw)
end

function suffstats{T<:Real}(::Type{Gamma}, x::AbstractArray{T})
    sx = 0.
    slogx = 0.
    for xi = x
        sx += xi
        slogx += log(xi)
    end
    GammaStats(sx, slogx, length(x))
end

function suffstats{T<:Real}(::Type{Gamma}, x::AbstractArray{T}, w::AbstractArray{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    sx = 0.
    slogx = 0.
    tw = 0.
    for i = 1:n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        sx += wi * xi
        slogx += wi * log(xi)
        tw += wi
    end
    GammaStats(sx, slogx, tw)
end

function gamma_mle_update(logmx::Float64, mlogx::Float64, a::Float64)
    ia = 1.0 / a
    z = ia + (mlogx - logmx + log(a) - digamma(a)) / (abs2(a) * (ia - trigamma(a)))
    1.0 / z
end

function fit_mle(::Type{Gamma}, ss::GammaStats;
    alpha0::Float64=NaN, maxiter::Int=1000, tol::Float64=1.0e-16)

    mx = ss.sx / ss.tw
    logmx = log(mx)
    mlogx = ss.slogx / ss.tw

    a::Float64 = isnan(alpha0) ? 0.5 / (logmx - mlogx) : alpha0
    converged = false

    t = 0
    while !converged && t < maxiter
        t += 1
        a_old = a
        a = gamma_mle_update(logmx, mlogx, a)
        converged = abs(a - a_old) <= tol
    end

    Gamma(a, mx / a)
end


text:
using ParallelUtilities

n = 1000
verbose=false
if verbose println("building $n x $n test matrix ...") end
@time A = Base.shmem_fill(zero(Float64),n,n)
if verbose println("done") end

# assertion fails for non-shared matrices, eg if A = zeros(n,n)

niters = 100

println("testing parallel repeat")
@time begin
    @sync begin
        @parallel_repeat niters for i=1:n
            A[:,i] += 1
        end
    end
end

@assert(A[5,5] == niters)
@assert(A[5,10] == niters)

println("testing parallel alternate")
@time begin
    @sync begin
        @parallel_alternate(niters,
            for i=1:n
                A[i,:] += i
            end,
            for i=1:n
                A[:,i] += i
            end
            )
    end
end

@assert(A[5,5] == 10*niters + niters)
@assert(A[5,10] == 15*niters + niters)


text:
import GLPK

using SparseArrays: spzeros

@testset "Minimum-cost flow" begin

    # bipartite oriented + source & sink
    # source 5, sink 6, v1 {1, 2} v2 {3, 4}
    g = Graphs.DiGraph(6)
    Graphs.add_edge!(g, 5, 1)
    Graphs.add_edge!(g, 5, 2)
    Graphs.add_edge!(g, 3, 6)
    Graphs.add_edge!(g, 4, 6)
    Graphs.add_edge!(g, 1, 3)
    Graphs.add_edge!(g, 1, 4)
    Graphs.add_edge!(g, 2, 3)
    Graphs.add_edge!(g, 2, 4)
    cost = zeros(6,6)
    cost[1,3] = 10.
    cost[1,4] = 5.
    cost[2,3] = 2.
    cost[2,4] = 2.
    # v2 -> sink have demand of one
    demand = spzeros(6,6)
    demand[3,6] = 1
    demand[4,6] = 1
    capacity = ones(6,6)

    o = GLPK.Optimizer
    node_demand = spzeros(Graphs.nv(g))
    flow = mincost_flow(g, node_demand, capacity, cost, o, edge_demand=demand, source_nodes=[5], sink_nodes=[6])
    @test flow[5,1] ≈ 1
    @test flow[5,2] ≈ 1
    @test flow[3,6] ≈ 1
    @test flow[4,6] ≈ 1
    @test flow[1,4] ≈ 1
    @test flow[2,4] ≈ 0
    @test flow[2,3] ≈ 1
    @test flow[1,3] ≈ 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1, nout] for nout in outneighbors(g, n1)) ≈ sum(flow[nin, n1] for nin in inneighbors(g, n1))
    end

    #--- Same test without edge_demands & positive cost should result in null flow
    flow = mincost_flow(g, node_demand, capacity, cost, o, source_nodes=[5], sink_nodes=[6])

    @test flow[5,1] ≈ 0
    @test flow[5,2] ≈ 0
    @test flow[3,6] ≈ 0
    @test flow[4,6] ≈ 0
    @test flow[1,4] ≈ 0
    @test flow[2,4] ≈ 0
    @test flow[2,3] ≈ 0
    @test flow[1,3] ≈ 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    #run same test again with nodal demands
    node_demand=spzeros(6)
    node_demand[5]=-2
    node_demand[6]=2
    #--- Same test without edge_demands
    flow = mincost_flow(g, node_demand, capacity, cost, o, source_nodes=[5],sink_nodes=[6])
    @test flow[5,1] ≈ 1
    @test flow[5,2] ≈ 1
    @test flow[3,6] ≈ 1
    @test flow[4,6] ≈ 1
    @test flow[1,4] ≈ 1
    @test flow[2,4] ≈ 0
    @test flow[2,3] ≈ 1
    @test flow[1,3] ≈ 0
    @test sum(diag(flow)) ≈ 0

    # flow conservation property
    for n1 = 1:4
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end

    # no demand => null flow
    d2 = spzeros(6,6)
    flow = mincost_flow(g, spzeros(Graphs.nv(g)), capacity, cost, o, edge_demand=d2, source_nodes=[5], sink_nodes=[6])
    for i in 1:6
        for j in 1:6
            @test flow[i,j] ≈ 0.0
        end
    end

    # graph without sink or source
    g = Graphs.DiGraph(6)
    # create circuit
    for s in 1:5
        Graphs.add_edge!(g, s, s+1)
    end
    Graphs.add_edge!(g, 6, 1)
    Graphs.add_edge!(g, 2, 5) # shortcut

    capacity = ones(6,6)
    demand = spzeros(6,6)
    demand[1,2] = 1
    costs = ones(6,6)
    flow = mincost_flow(g,spzeros(Graphs.nv(g)), capacity, costs, o, edge_demand=demand)
    active_flows = [(1,2), (2,5), (5,6),(6,1)]
    for s in 1:6
        for t in 1:6
            if (s,t) in active_flows
                @test flow[s,t] ≈ 1.
            else
                @test flow[s,t] ≈ 0.
            end
        end
    end
    # flow conservation property
    for n1 = 1:6
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end
    # higher short-circuit cost
    costs[2,5] = 10.
    flow = mincost_flow(g,spzeros(Graphs.nv(g)), capacity, costs, o, edge_demand=demand)
    active_flows = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1)]
    for s in 1:6
        for t in 1:6
            if (s,t) in active_flows
                @test flow[s,t] ≈ 1.
            else
                @test flow[s,t] ≈ 0.
            end
        end
    end
    # flow conservation property
    for n1 = 1:6
        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
    end
end


text:
# Depth-first visit / traversal

#################################################
#
#  Depth-first visit
#
#################################################

type DepthFirst <: AbstractGraphVisitAlgorithm
end

function depth_first_visit_impl!{V,E}(
    graph::AbstractGraph{V,E},      # the graph
    stack,                          # an (initialized) stack of vertex
    vertexcolormap::Vector{Int},    # an (initialized) color-map to indicate status of vertices
    edgecolormap::Vector{Int},      # an (initialized) color-map to indicate status of edges
    visitor::AbstractGraphVisitor)  # the visitor

    while !isempty(stack)
        u, uegs, tstate = pop!(stack)
        found_new_vertex = false

        while !done(uegs, tstate) && !found_new_vertex
            v_edge, tstate = next(uegs, tstate)
            v = v_edge.target
            v_color = vertexcolormap[vertex_index(v, graph)]
            e_color = edgecolormap[edge_index(v_edge, graph)]
            examine_neighbor!(visitor, u, v, v_color, e_color)

            if e_color == 0
                edgecolormap[edge_index(v_edge, graph)] = 1
            end

            if v_color == 0
                found_new_vertex = true
                vertexcolormap[vertex_index(v, graph)] = 1
                if !discover_vertex!(visitor, v)
                    return
                end
                push!(stack, (u, uegs, tstate))

                open_vertex!(visitor, v)
                vegs = out_edges(v, graph)
                push!(stack, (v, vegs, start(vegs)))
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vertexcolormap[vertex_index(u, graph)] = 2
        end
    end
end

function traverse_graph{V,G <: AbstractGraph}(
    graph::G,
    alg::DepthFirst,
    s::V,
    visitor::AbstractGraphVisitor;
    vertexcolormap = zeros(Int, num_vertices(graph)),
    edgecolormap = zeros(Int, num_edges(graph)))

    @graph_requires graph incidence_list vertex_map

    vertexcolormap[vertex_index(s, graph)] = 1
    if !discover_vertex!(visitor, s)
        return
    end

    segs = out_edges(s, graph)
    sstate = start(segs)
    stack = [(s, segs, sstate)]

    depth_first_visit_impl!(graph, stack, vertexcolormap, edgecolormap, visitor)
end

#################################################
#
#  Useful applications
#
#################################################

# Test whether a graph is cyclic

type DFSCyclicTestVisitor <: AbstractGraphVisitor
    found_cycle::Bool

    DFSCyclicTestVisitor() = new(false)
end

function examine_neighbor!{V}(
    vis::DFSCyclicTestVisitor,
    u::V,
    v::V,
    vcolor::Int,
    ecolor::Int)

    if vcolor == 1 && ecolor == 0
        vis.found_cycle = true
    end
end

discover_vertex!(vis::DFSCyclicTestVisitor, v) = !vis.found_cycle

function test_cyclic_by_dfs{G <: AbstractGraph}(graph::G)
    @graph_requires graph vertex_list incidence_list vertex_map

    cmap = zeros(Int, num_vertices(graph))
    visitor = DFSCyclicTestVisitor()

    for s in vertices(graph)
        if cmap[vertex_index(s, graph)] == 0
            traverse_graph(graph, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end

        if visitor.found_cycle
            return true
        end
    end
    return false
end

# Topological sort using DFS

type TopologicalSortVisitor{V} <: AbstractGraphVisitor
    vertices::Vector{V}

    function TopologicalSortVisitor(n::Int)
        vs = Array(Int, 0)
        sizehint!(vs, n)
        new(vs)
    end
end

function examine_neighbor!{V}(visitor::TopologicalSortVisitor{V}, u::V, v::V, vcolor::Int, ecolor::Int)
    if vcolor == 1 && ecolor == 0
        throw(ArgumentError("The input graph contains at least one loop."))
    end
end

function close_vertex!{V}(visitor::TopologicalSortVisitor{V}, v::V)
    push!(visitor.vertices, v)
end

function topological_sort_by_dfs{V}(graph::AbstractGraph{V})
    @graph_requires graph vertex_list incidence_list vertex_map

    cmap = zeros(Int, num_vertices(graph))
    visitor = TopologicalSortVisitor{V}(num_vertices(graph))

    for s in vertices(graph)
        if cmap[vertex_index(s, graph)] == 0
            traverse_graph(graph, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end
    end

    reverse(visitor.vertices)
end


text:
A = randn(3,5)
b = randn(3)

x_pg = Variable(5)
prob_pg = problem(ls(A*x_pg - b) + 1e-3*norm(x_pg, 1))
sol_pg = solve(prob_pg, ForwardBackward())

x_zfpr = Variable(5)
prob_zfpr = problem(ls(A*x_zfpr - b) + 1e-3*norm(x_zfpr, 1))
sol_zfpr = solve(prob_zfpr, ZeroFPR())

x_pnc = Variable(5)
prob_pnc = problem(ls(A*x_pnc - b) + 1e-3*norm(x_pnc, 1))
sol_pnc = solve(prob_pnc, PANOC())


text:
module aoc_15

using Test

export read_input, test_part1, solve_part1, test_part2, solve_part2

function read_input(filename)
    lines = readlines(filename)
    numbers = split(lines[1],",")
    return [parse(Int,x) for x in numbers]
end

function test_part1()
    starting = [0,3,6]
    @test solve_part1(starting) == 436
end

function solve_part1(starting = read_input("input.txt"), finish = 2020)
    memory = Dict{Int,Int}()
    for i = 1:length(starting)
        memory[ starting[ i ] ] = i
    end
    last = starting[end]
    age = 0
    for i = length( starting ) : finish - 1
        last_turn = get( memory, last, i )
        memory[ last ] = i
        age = i - last_turn
        last = age
    end
    return age
end

# Part 2

function test_part2()
#=     Given 0,3,6, the 30000000th number spoken is 175594.
Given 1,3,2, the 30000000th number spoken is 2578.
Given 2,1,3, the 30000000th number spoken is 3544142.
Given 1,2,3, the 30000000th number spoken is 261214.
Given 2,3,1, the 30000000th number spoken is 6895259.
Given 3,2,1, the 30000000th number spoken is 18.
Given 3,1,2, the 30000000th number spoken is 362. =#
    @test solve_part2([0,3,6]) == 175594
    @test solve_part2([1,3,2]) == 2578
    @test solve_part2([2,1,3]) == 3544142
    @test solve_part2([1,2,3]) == 261214
    @test solve_part2([2,3,1]) == 6895259
    @test solve_part2([3,2,1]) == 18
    @test solve_part2([3,1,2]) == 362
end

function solve_part2(starting = read_input("input.txt"))
    return solve_part1(starting, 30000000)
end

end


text:
module Extra_Test

using FactCheck, JNeuron

myimport=input(string(dirname(Base.source_path()),"/../examples/data/cell2.asc"));
blank_neuron=instantiate(myimport);
set_nsegs!(blank_neuron);

myneuron=add(blank_neuron,(HH(),Passive()));
myneuron2=add(blank_neuron,HH());

mynetwork=Network(myneuron,15.0);

electrode1_1=Extracellular([0.0,0.0,0.0]);
add!(mynetwork,electrode1_1);

electrode2_1=Extracellular(JNeuron.Point(),[0.0,0.0,0.0])
add!(mynetwork,electrode2_1);

electrode3_1=Extracellular(JNeuron.Mixed(),[0.0,0.0,0.0])
add!(mynetwork,electrode3_1);

mynetwork2=Network((deepcopy(myneuron),myneuron2),15.0)

electrode1_2=Extracellular([0.0,0.0,0.0]);
add!(mynetwork2,electrode1_2);

electrode2_2=Extracellular(JNeuron.Point(),[0.0,0.0,0.0])
add!(mynetwork2,electrode2_2);

electrode3_2=Extracellular(JNeuron.Mixed(),[0.0,0.0,0.0])
add!(mynetwork2,electrode3_2);

run!(mynetwork,true);
run!(mynetwork2,true);

facts() do

    for i=1:3
        @fact length(mynetwork.extra[i].v) --> 601
        @fact length(mynetwork.extra[i].coeffs) --> 1
        @fact length(mynetwork.extra[i].coeffs[1].c) --> 726
        @fact length(mynetwork.extra[i].coeffs[1].inds) --> 726
        @fact mynetwork.extra[i].v[end] --> roughly(0.0, atol=1e-4)
    end

    for i=1:3
        @fact length(mynetwork2.extra[i].coeffs) --> 2
        @fact length(mynetwork2.extra[i].v) --> 601
        @fact mynetwork2.extra[i].v[end] --> roughly(0.0, atol=1e-4)
        for j=1:2
            @fact length(mynetwork2.extra[i].coeffs[j].c) --> 726
            @fact length(mynetwork2.extra[i].coeffs[j].inds) --> 726
        end
    end

end

mynetwork3=Network(deepcopy(myneuron),15.0)
mystim=Stim(5.0,1,1,924,5.0,5.125)
add!(mynetwork3,mystim);
i3=JNeuron.runc(mynetwork3,true);

mye=Extracellular([500.0,125.0,0.0])
myv=JNeuron.nete(myneuron,i3,mye,100);
(myv3, mys)=JNeuron.extrap(myv,50.0);

facts() do

    @fact size(i3) --> (925,601)
    @fact size(myv) --> (582,100)
    @fact size(myv3) --> (2001,)

end

mye=Extracellular([500.0,125.0,0.0])
mye2=Extracellular([500.0,75.0,0.0])
myv=JNeuron.nete(myneuron,i3,[mye,mye2],100);
(myv3, mys)=JNeuron.extrap(myv,50.0);

facts() do

    @fact size(myv) --> (582,100,2)
    @fact size(myv3) --> (2001,2)

end

end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
"""
`MyDictionary` is a custom dictionary type that illustrates the minimal
functionality required to implement a dictionary.
"""
struct MyDictionary{S,T} <: Dictionary{S,T}
    n       ::  Int
end

# We provide a constructor that will choose a default length and default domain type.
MyDictionary(n = 10) = MyDictionary{Float64}(n)
# For convenience, we automatically choose a codomaintype if only a domaintype is given.
MyDictionary{S}(n) where {S} = MyDictionary{S,S}(n)

# The `length` function is implemented in terms of size.
# The size of a dictionary is always a tuple, just like for arrays.
size(dict::MyDictionary) = (dict.n,)

# By implementing `similar`, we support promotion of the domain type, as well
# as resizing the dictionary.
# If resizing is not possible, one can add a statement such as
# "@assert n == length(dict)", which will cause an error upon any user attempt to resize.
# Same for promotion.
similar(dict::MyDictionary, ::Type{S}, n::Int) where {S} = MyDictionary{S}(n)

support(dict::MyDictionary) = UnitInterval{domaintype(dict)}()

unsafe_eval_element(dict::MyDictionary, idx, x) = x^(idx-1)


text:
using WiltonInts84

abstract type RetardedPotential{T} <: Operator end
Base.eltype(::RetardedPotential{T}) where {T} = T
scalartype(A::RetardedPotential{T}) where {T} = T

mutable struct EmptyRP{T} <: RetardedPotential{T}
    speed_of_light::T
end
Base.eltype(::EmptyRP) = Int
quaddata(op::EmptyRP, xs...) = nothing
quadrule(op::EmptyRP, xs...) = nothing
momintegrals!(z, op::EmptyRP, xs...) = nothing

function allocatestorage(op::RetardedPotential, testST, basisST,
	::Type{Val{:densestorage}},
	::Type{LongDelays{:ignore}})

    tfs = spatialbasis(testST)
    bfs = spatialbasis(basisST)

    M = numfunctions(tfs)
    N = numfunctions(bfs)

    K0 = zeros(Int, M, N)
    K1 = zeros(Int, M, N)

    function store(v,m,n,k)
        K0[m,n] = (K0[m,n] == 0) ? K0[m,n] : min(K0[m,n],k)
        K1[m,n] = max(K1[m,n],k)
    end

    aux = EmptyRP(op.speed_of_light)
    print("Allocating memory for convolution operator: ")
    assemble!(aux, testST, basisST, store)
    println("\nAllocated memory for convolution operator.")

    kmax = maximum(K1);
	Z = zeros(eltype(op), M, N, kmax)
    store1(v,m,n,k) = (Z[m,n,k] += v)
    return ()->MatrixConvolution(Z), store1
end

# function allocatestorage(op::RetardedPotential, testST, basisST,
# 	::Type{Val{:bandedstorage}},
# 	::Type{LongDelays{:ignore}})

#     tfs = spatialbasis(testST)
#     bfs = spatialbasis(basisST)

#     M = numfunctions(tfs)
#     N = numfunctions(bfs)

#     K0 = fill(typemax(Int), M, N)
#     K1 = zeros(Int, M, N)

#     function store(v,m,n,k)
#         K0[m,n] = min(K0[m,n],k)
#         K1[m,n] = max(K1[m,n],k)
#     end

#     aux = EmptyRP(op.speed_of_light)
#     print("Allocating memory for convolution operator: ")
#     assemble!(aux, testST, basisST, store)
#     println("\nAllocated memory for convolution operator.")

# 	maxk1 = maximum(K1)
# 	bandwidth = maximum(K1 .- K0 .+ 1)
# 	data = zeros(eltype(op), bandwidth, M, N)
# 	Z = SparseND.Banded3D(K0, data, maxk1)
#     store1(v,m,n,k) = (Z[m,n,k] += v)
#     return ()->Z, store1
# end

struct Storage{T} end

function allocatestorage(op::RetardedPotential, testST, basisST,
	::Type{Val{:bandedstorage}},
    ::Type{LongDelays{:compress}})

    @info "Allocating mem for RP op compressing the static tail..."

	T = eltype(op)

    tfs = spatialbasis(testST)
    bfs = spatialbasis(basisST)

	Δt = timestep(temporalbasis(basisST))
	Nt = numfunctions(temporalbasis(basisST))

	tbf = convolve(temporalbasis(testST), temporalbasis(basisST))
	has_tail = !all(tbf.polys[end].data .== 0)

    M = numfunctions(tfs)
    N = numfunctions(bfs)

    K0 = fill(typemax(Int), M, N)
    K1 = zeros(Int, M, N)

    function store_alloc(v,m,n,k)
        K0[m,n] = min(K0[m,n],k)
        K1[m,n] = max(K1[m,n],k)
    end

    op_alloc = EmptyRP(op.speed_of_light)
	tbf_trunc = truncatetail(tbf)
	δ = timebasisdelta(Δt, Nt)
    print("Allocating memory for convolution operator: ")
    assemble!(op_alloc, tfs⊗δ, bfs⊗tbf_trunc, store_alloc)
    println("\nAllocated memory for convolution operator.")

	bandwidth = maximum(K1 .- K0 .+ 1)
	data = zeros(T, bandwidth, M, N)
    tail = zeros(T, M, N)
    # kmax = maximum(K1)
    len = has_tail ? Nt : maximum(K1)
	Z = ConvOp(data, K0, K1, tail, len)

    function store1(v,m,n,k)
		if Z.k0[m,n] ≤ k ≤ Z.k1[m,n]
			Z.data[k - Z.k0[m,n] + 1,m,n] += v
		elseif k == Z.k1[m,n]+1
			Z.tail[m,n] += v
		end
	end

    return ()->Z, store1
end

function assemble!(op::LinearCombinationOfOperators, tfs::SpaceTimeBasis, bfs::SpaceTimeBasis, store,
    threading=Threading{:multi})

    for (a,A) in zip(op.coeffs, op.ops)
        store1(v,m,n,k) = store(a*v,m,n,k)
        assemble!(A, tfs, bfs, store1)
    end
end

function assemble!(op::RetardedPotential, testST, trialST, store,
    threading=Threading{:multi})

	Y, S = spatialbasis(testST), temporalbasis(testST)

	P = Threads.nthreads()
	@assert P >= 1
	splits = [round(Int,s) for s in range(0, stop=numfunctions(Y), length=P+1)]

	# @info "Starting assembly with $P threads:"
	# Threads.@threads for i in 1:P
	# 	lo, hi = splits[i]+1, splits[i+1]
	# 	lo <= hi || continue
	# 	Y_p = subset(Y, lo:hi)
	# 	store1 = (v,m,n,k) -> store(v,lo+m-1,n,k)
	# 	assemble_chunk!(op, Y_p ⊗ S, trialST, store1)
	# end

	assemble_chunk!(op, testST, trialST, store)
end

function assemble_chunk!(op::RetardedPotential, testST, trialST, store)

	myid = Threads.threadid()

    testspace  = spatialbasis(testST)
    trialspace = spatialbasis(trialST)
    timebasisfunction = convolve(temporalbasis(testST), temporalbasis(trialST))

	testels, testad   = assemblydata(testspace)
	trialels, trialad = assemblydata(trialspace)

	speedoflight = op.speed_of_light
	Δt = timestep(timebasisfunction)
	ct, hs = boundingbox(vertices(geometry(trialspace)))
	diam = 2 * sqrt(3) * hs
	#kmax = ceil(Int, diam/speedoflight/timestep(timebasisfunction)) + (numintervals(timebasisfunction)-1)
	kmax = ceil(Int, (numintervals(timebasisfunction)-1) + diam/speedoflight/Δt)+1
	kmax = max(kmax, numfunctions(timebasisfunction))
    timead = temporalassemblydata(timebasisfunction, kmax=kmax)

	Δt = timestep(timebasisfunction)
	ΔR = Δt * speedoflight

	Nt = numfunctions(timebasisfunction)
	tmax = (Nt-1) * Δt

    U = refspace(testspace)
    V = refspace(trialspace)
    W = refspace(timebasisfunction)

    qd = quaddata(op, U, V, W, testels, trialels, nothing)

    udim = numfunctions(U)
    vdim = numfunctions(V)
    wdim = numfunctions(W)
    z = zeros(eltype(op), udim, vdim, wdim)

	@show length(testels) length(trialels)

    myid == 1 && print("dots out of 10: ")
    todo, done, pctg = length(testels), 0, 0
    for p in eachindex(testels)
        τ = testels[p]
        for q in eachindex(trialels)
            σ = trialels[q]
	        for r in rings(τ,σ,ΔR)
				r > numfunctions(timebasisfunction) && continue
	            ι = ring(r,ΔR)

	            # compute interactions between reference shape functions
	            fill!(z, 0)
	            qr = quadrule(op, U, V, W, p, τ, q, σ, r, ι, qd)
                momintegrals!(z, op, U, V, W, τ, σ, ι, qr)

		        # assemble in the global matrix
		        for i in 1 : udim
		            for j in 1 : vdim
		                for d in 1 : wdim

		                    v = z[i,j,d]

                            for (m,a) in testad[p,i]
                                av = a*v
                                for (n,b) in trialad[q,j]
									abv = b*av
									tad_rd = timead[r,d]
                                    for (k,c) in tad_rd
										store(c*abv, m, n, k)
									end # next κ
		                        end # next ν
		                    end # next μ
		                end # next d
		            end # next j
		        end #next i
		    end # next r
		end # next q

        done += 1
        new_pctg = round(Int, done / todo * 100)
        if myid == 1 && new_pctg > pctg + 9
            print(".")
            pctg = new_pctg
        end
    end # next p

    println("")
end

struct WiltonInts84Strat{T,V,W}
	outer_quad_points::T
    binomials::V
    workspace::W
end

function momintegrals!(z, op, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)

    XW = qr.outer_quad_points
    for p in 1 : length(XW)
        x = XW[p].point
        w = XW[p].weight
		innerintegrals!(z, op, x, g, f, T, τ, σ, ι, qr, w)
    end

end


text:
"""
    Learn Dash server

Use `run_Dash.jl` to run the Dash server.
For learning or debugging purpose,

```julia-repl
julia> include("learn_Dash.jl")
julia> db.@enter your_function_name()
# Enter compiled mode by `C` (capital C), otherwise it is too slow.
1|debug> C
```
"""

# Use run_Dash.jl to run Dash server
#
# For debugging,
#

# Deprecatd, use run_Dash.jl
#
# Use version 1.5.3 and run in command line
#
#   julia --project=@. learn_Dash.jl
#
# this will enable hot reloading, meaning the server refresh the page when
# learn_Dash.jl is updated.
#
# Run in interactive mode is fine by
#
#   include("learn_Dash.jl")
#
# but needs to shutdown and restart the server every time the file updated.

# https://dash-julia.plotly.com
# add Dash DashCoreComponents DashHtmlComponents DashTable

# https://dash-julia.plotly.com/getting-started
using Dash, DashHtmlComponents, DashCoreComponents, DashBootstrapComponents
using DataFrames, CSV, HTTP
using PlotlyJS #, RDatasets
import JSON3

import Debugger
db = Debugger

using Logging
logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

# https://dash-julia.plotly.com/getting-started

"""
    dataset(dir, name)

Replace the dataset function from RDatasets.  It stopped working.
"""
function dataset(dir, name)
    url = "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/$(dir)/$(name).csv"
    body = HTTP.get(url).body
    csv = CSV.File(body)
    df = DataFrame(csv)
end

function urldownload(url)
    body = HTTP.get(url).body
    csv = CSV.File(body)
    df = DataFrame(csv)
end

function hello_dash()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        h1 = html_h1("Hello Dash",
            style = Dict("color" => "#7FDBFF", "textAlign" => "center")
        )
        d = html_div("Dash: A web application")
        dcc = dcc_graph(
            id = "example-graph-1",
            figure = (
                data = [
                    (
                        x = ["giraffes", "orangutans", "monkeys"],
                        y = [20, 14, 23], type = "bar", name = "SF"
                    ),
                    (
                        x = ["giraffes", "orangutans", "monkeys"],
                        y = [12, 18, 29], type = "bar", name = "Montreal"
                    )
                ],
                layout = (
                    title = "Dash Data Visualization", barmode="group",
                    plot_bgcolor = "#111111",
                    paper_bgcolor = "#114411",
                    font = Dict("color" => "#7FDBFF")
                )
            )
        )

        # The original example uses comma to separate the 3 tags, and return
        # the value that way.

        return (h1, d, dcc)
    end

    run_server(app, "0.0.0.0", debug=true)
end

function generate_table(df, max_rows=10)
    t1 = html_table([
        html_thead(html_tr([html_th(col) for col in names(df)])),
        html_tbody([
            html_tr([html_td(df[r,c]) for c in names(df)]) for r = 1:min(nrow(df),max_rows)
        ])
    ])
end

function dash_table()
    url = "https://raw.githubusercontent.com/plotly/datasets/master/2011_us_ag_exports.csv"
    #csv = urldownload("https://raw.githubusercontent.com/plotly/datasets/master/2011_us_ag_exports.csv")
    csv = CSV.File(HTTP.get(url).body)
    @info typeof(csv)
    df = DataFrame(csv)

    #app = dash()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        html_h4("US Agriculture Exports (2011)"),
        generate_table(df, 10)
    end

    run_server(app, "0.0.0.0", debug=true)
end

"""
    dash_plotly()

Demonstrate simple Plotly plot in Dash.  Use `Layout` to specify figure size.
Use my own `dataset` function to retrieve iris data from R.

[More About Visualization](https://dash-julia.plotly.com/getting-started)
"""
function dash_plotly()
    iris = dataset("datasets", "iris")
    layout = Layout(title="Iris", width=800, height=600)

    plt = Plot(
        iris, layout,
        x=Symbol("Sepal.Length"), y=Symbol("Sepal.Width"), mode="markers",
        marker_size=8, group=:Species,
    )
    @info typeof(plt)
    @info typeof(plt.layout)
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        html_h4("Iris Sepal Length vs Sepal Width"),
        dcc_graph(
            id = "example-graph-3",
            figure = plt,
        )
    end

    run_server(app, "0.0.0.0", debug=true)

    return plt
end

function dash_markdown()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    markdown_text = """
    ### Dash and Markdown

    Dash apps can be written in Markdown.
    Dash uses the [CommonMark](http://commonmark.org/)
    specification of Markdown.
    Check out their [60 Second Markdown Tutorial](http://commonmark.org/help/)
    if this is your first introduction to Markdown!

    """

    app.layout = html_div() do
        dcc_markdown(markdown_text)
    end

    run_server(app, "0.0.0.0", debug=true)
end

function dash_core()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    dropdown_options = [
        Dict("label" => "New York City", "value" => "NYC"),
        Dict("label" => "Montreal", "value" => "MTL"),
        Dict("label" => "San Francisco", "value" => "SF"),
    ]
    app.layout = html_div(style = Dict("columnCount" => 2)) do
        html_label("Dropdown"),
        dcc_dropdown(options = dropdown_options, value = "MTL"),
        html_label("Multi-Select Dropdown"),
        dcc_dropdown(
            options = dropdown_options,
            value = ["MTL", "SF"],
            multi = true,
        ),
        html_label("Radio Items"),
        dcc_radioitems(options = dropdown_options, value = "MTL"),
        html_label("Checkboxes"),
        dcc_checklist(options = dropdown_options, value = ["MTL", "SF"]),
        html_label("Text Input"),
        dcc_input(value = "MTL", type = "text"),
        html_label("Slider"),
        dcc_slider(
            min = 0,
            max = 9,
            marks = Dict([i => (i == 1 ? "Label $(i)" : "$(i)") for i = 1:6]),
            value = 5,
        )
    end

    run_server(app, "0.0.0.0", debug = true)
end

"""
    dash_basic_callback()

[Basic Dash Callback](https://dash-julia.plotly.com/basic-callbacks)
"""
function dash_basic_callback()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        dcc_input(id = "input-3", value = "initial value",
            type = "text",
        ),
        html_div(id="output-1")
    end

    callback!(
        app,
        Output("output-1", "children"),
        Input("input-3", "value")) do input_value
        "You have entered $(input_value)"
    end

    run_server(app, "0.0.0.0", debug = true)
end

"""

[Dash App Layout With Figure and Slider](https://dash-julia.plotly.com/basic-callbacks)
"""
function dash_layout_figure_slider()
    df1 = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/gapminderDataFiveYear.csv"))
    years = unique(df1[!, :year])
    @show years
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])
    app.layout = html_div() do
        dcc_graph(id = "graph-1"),
        dcc_slider(
            id = "year-slider-1",
            min = minimum(years),
            max = maximum(years),
            marks = Dict([Symbol(v) => Symbol(v) for v in years]),
            value = minimum(years),
            step = nothing,
        )
    end

    callback!(
        app,
        Output("graph-1", "figure"),    # figure is an attribute of graph-1
        Input("year-slider-1", "value") # value is an attribut of year-slider-1
        ) do selected_year
        pt = Plot(
            df1[df1.year .== selected_year, :],
            Layout(
                xaxis_type = "log",
                xaxis_title = "GDP Per Capita",
                yaxis_title = "Life Expectancy",
                legend_x = 0,
                legend_y = 1,
                hovermode = "closest",
                transition_duration = 500
            ),
            x = :gdpPercap,
            y = :lifeExp,
            text = :country,
            group = :continent,
            mode = "markers",
            marker_size = 15,
            marker_line_color = "white"
        )

        return pt
    end

    run_server(app, "0.0.0.0", debug=true)
end

"""

[Dash App With Multiple Inputs](https://dash-julia.plotly.com/basic-callbacks)
"""
function dash_app_multiple_inputs()
    df = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/country_indicators.csv"))
    #@show size(df)
    #@show names(df)
    dropmissing!(df)
    #display(df)

    # [!, ] returns a reference rather than a copy
    # https://stackoverflow.com/questions/60900001/what-is-the-meaning-of-the-exclamation-mark-in-indexing-a-julia-dataframe
    available_indicators = unique(df[!, "Indicator Name"])
    years = unique(df[!, "Year"])

    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        div_1 = html_div(
            children = [
                dcc_dropdown(
                    id = "xaxis-column",
                    options = [(label=i, value=i) for i in available_indicators],
                    value = available_indicators[3] #"Fertility rate, total (births per woman)"
                ),
                dcc_radioitems(
                    id = "xaxis-type",
                    options = [(label=v, value=v) for v in ["linear", "log"]],
                    value = "linear"
                )
            ],
            style = (width="48%", display="inline-block")
        )
        div_2 = html_div(
            children=[
                dcc_dropdown(
                    id = "yaxis-column",
                    options = [(label=v, value=v) for v in available_indicators],
                    value = available_indicators[4]
                ),
                dcc_radioitems(
                    id = "yaxis-type",
                    options = [(label = v, value = v) for v in ["linear", "log"]],
                    value = "linear",
                )
            ],
            style = (width="48%", display="inline-block", float = "right")

        )
        graph = dcc_graph(id = "indicator-graphic")
        slider = dcc_slider(
            id = "year-slider-2",
            min = minimum(years), max = maximum(years),
            marks = Dict([Symbol(v) => Symbol(v) for v in years]),
            value = minimum(years),
            step = nothing,
        )

        return div_1,div_2, graph, slider
    end

    callback!(
        app,
        Output("indicator-graphic", "figure"),
        Input("xaxis-column", "value"),
        Input("yaxis-column", "value"),
        Input("xaxis-type", "value"),
        Input("yaxis-type", "value"),
        Input("year-slider-2", "value"),
    ) do xaxis_column_name, yaxis_column_name, xaxis_type, yaxis_type, year_value
        dff = df[df.Year .== year_value, :]

        pt = Plot(
            dff[dff[!, Symbol("Indicator Name")] .== xaxis_column_name, :Value],
            dff[dff[!, Symbol("Indicator Name")] .== yaxis_column_name, :Value],
            Layout(
                xaxis_type = xaxis_type == "linear" ? "linear" : "log",
                xaxis_title = xaxis_column_name,
                yaxis_type = yaxis_type == "linear" ? "linear" : "log",
                yaxis_title = yaxis_column_name,
                hovermode = "closest",
            ),
            kind = "scatter",
            text = dff[
                dff[!, Symbol("Indicator Name")] .== yaxis_column_name,
                Symbol("Country Name")
            ],
            mode = "markers",
            marker_size = 15,
            marker_opacity = 0.5,
            marker_line_width = 1,
            marker_line_color = "white"
        )

        return pt
    end

    run_server(app, "0.0.0.0", debug=true)
end
"""

[Dash App With Multiple Outputs](https://dash-julia.plotly.com/basic-callbacks)
"""
function dash_app_multiple_outputs()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        dcc_input(id = "input-1", value = "1", type="text"),
        html_tr(
            (html_td("x^2 ="), html_td(id="square"))
        ),
        html_tr(
            (html_td("x^3 ="), html_td(id="cube"))
        )
    end

    callback!(
        app,
        Output("square", "children"),
        Output("cube", "children"),
        Input("input-1", "value")
    ) do x
        @show x
        try
            x = parse(Int64, x)
            return (x^2, x^3)
        catch e
            if isa(e, ArgumentError)
                return ("", "")
            end

            # unknown error
            @error e
        end

    end

    run_server(app, "0.0.0.0", debug=true)
end

"""
[Dash App With Chained Callbacks]
"""
function dash_app_chained_callbacks()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    all_options = Dict(
        "America" => ["New York City", "San Francisco", "Cincinnati"],
        "Canada"  => ["Montreal", "Toronto", "Ottawa"],
    )
    app.layout = html_div() do
        html_div(
            children = [
                dcc_radioitems(
                    id = "countries-radio",
                    options = [(label = x, value = x) for x in keys(all_options)],
                    value = "America", # keys(all_options)[1]
                ),
                html_hr(),
                dcc_radioitems(id = "cities-radio"),
                html_hr(),
                html_div(id = "display-selected-values"),
            ],
        )
    end

    callback!(
        app,
        Output("cities-radio", "options"),
        Input("countries-radio", "value")
    ) do selected_country
        cities = all_options[selected_country]
        Ans = [(label = x, value = x) for x in cities]

        return Ans
    end
    callback!(
        app,
        Output("cities-radio", "value"),
        Input("cities-radio", "options")
    ) do available_options
        #@show available_options
        #return available_options[1].value
        return available_options[1][:value]
    end

    callback!(
        app,
        Output("display-selected-values", "children"),
        Input("countries-radio", "value"),
        Input("cities-radio", "value"),
    ) do selected_country, selected_city
        if isnothing(selected_city)
            return "No city selected"
        else
            return "$(selected_city) is a city in $(selected_country)"
        end
    end

    run_server(app, "0.0.0.0", debug=true)
end

function dash_app_state()
    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        input_1 = dcc_input(id = "input-1-state", type = "text", value = "Montreal")
        input_2 = dcc_input(id = "input-2-state", type = "text", value = "Canada")
        button = html_button(id = "submit-button-state", children = "submit", n_clicks = 0)
        out   = html_div(id = "output-state")
        return input_1, input_2, button, out
    end

    callback!(
        app,
        Output("output-state", "children"),
        Input("submit-button-state", "n_clicks"),
        State("input-1-state", "value"),
        State("input-2-state", "value"),
    ) do clicks, input_1, input_2
        Ans = "The button has been clicked $(clicks) times, inputs are $(input_1) and $(input_2)"

        return Ans
    end

    run_server(app, "0.0.0.0", debug=true)
end

"""
[Interactive Visualizations](https://dash-julia.plotly.com/interactive-graphing)
"""
function dash_intro_visual()
    df = DataFrame(
        x = [1,2,1,2],
        y = [1,2,3,4],
        customdata = [1,2,3,4],
        fruit = ["apple", "apple", "orange", "orange"],
    )

    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])
    app.layout = html_div() do
        fig = dcc_graph(id = "basic-interactions", figure = (
            data = [
                (
                    x = [1,2,3,4],
                    y = [4,1,3,5],
                    text = ["a", "b", "c", "d"],
                    customdata = ["c.a", "c.b", "c.c", "c.d"],
                    name = "Trace-1",
                    mode = "markers",
                    marker = (size=12,)
                ),
                (
                    x = [1,2,3,4],
                    y = [9,4,1,4],
                    customdata = ["c.w", "c.x", "c.y", "c.z"],
                    name = "Trace-2",
                    mode = "markers",
                    marker = (size=12,)
                )
            ],
            layout = (clickmode = "event+select",)
        ))
        tag = html_div(
            children = [
                html_div(
                    children = [
                        dcc_markdown("
                        **Hover Data**

                        Mouse over values in the graph.
                        "),
                        html_pre(id = "hover-data"),
                    ],
                ),
                html_div(
                    children = [
                        dcc_markdown("
                        **Click Data**

                        Click on points in the graph.
                        "),
                        html_pre(id = "click-data"),
                    ],
                ),
                html_div(
                    children = [
                        dcc_markdown("
                        **Selection Data**

                        Choose the lasso or rectangle tool in the graph's menu
                        bar and then select points in the graph.

                        Note that if `layout.clickmode = 'event+select'`, selection data also
                        accumulates (or un-accumulates) selected data if you hold down the shift
                        button while clicking.
                        "),
                        html_pre(id = "selected-data"),
                    ],
                ),
                html_div(
                    children = [
                        dcc_markdown("
                        **Zoom and Relayout Data**

                        Click and drag on the graph to zoom or click on the zoom
                        buttons in the graph's menu bar.
                        Clicking on legend items will also fire
                        this event.
                        "),
                        html_pre(id = "relayout-data"),
                    ],

                ),
            ]
        )

        return fig, tag
    end

    # Hover
    callback!(
        app,
        Output("hover-data", "children"),
        Input("basic-interactions", "hoverData"),
    ) do hover_data
        @show typeof(hover_data)
        @show hover_data

        return JSON3.write(hover_data)
    end

    # Click
    callback!(
        app,
        Output("click-data", "children"),
        Input("basic-interactions", "clickData"),
    ) do click_data
        @show click_data

        return JSON3.write(click_data)
    end

    # Select
    callback!(
        app,
        Output("selected-data", "children"),
        Input("basic-interactions", "selectedData"),
    ) do selected_data
        @show selected_data

        return JSON3.write(selected_data)
    end

    # Re-layout
    callback!(
        app,
        Output("relayout-data", "children"),
        Input("basic-interactions", "relayoutData"),
    ) do relayout_data
        @show relayout_data
        return JSON3.write(relayout_data)
    end

    #return df
    run_server(app, "0.0.0.0", debug=true)
end

"""
[Update Graphs on Hover](https://dash-julia.plotly.com/interactive-graphing)
"""
function dash_update_graphs_on_hover()
    println("Random command to study Debugger")

    df = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/country_indicators.csv"))

    dropmissing!(df)
    #display(df)
    #@show names(df)

    available_indicators = unique(df[:, "Indicator Name"])
    years = unique(df[!, "Year"])

    #app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])
    app = dash(external_stylesheets=[dbc_themes.BOOTSTRAP])
    app.layout = html_div() do
        tag_1 = html_div(
            children = [
                html_div( # x
                    children=[
                        dcc_dropdown( # x-column
                            id = "crossfilter-xaxis-column",
                            options = [(label = x, value = x) for x in available_indicators],
                            value = available_indicators[3],
                        ),
                        dcc_radioitems( # x-axis
                            id = "crossfilter-xaxis-type",
                            options = [(label = x, value = x) for x in ["linear", "log"]],
                            value = "linear",
                        )
                    ],
                    style = (width = "49%", display = "inline-block")
                ),
                html_div( # y
                    children = [
                        dcc_dropdown( # y-column
                            id = "crossfilter-yaxis-column",
                            options = [(label = x, value = x) for x in available_indicators],
                            value = available_indicators[4],
                        ),
                        dcc_radioitems( # y-axis
                            id = "crossfilter-yaxis-type",
                            options = [(label = x, value = x) for x in ["linear", "log"]],
                            value = "linear",
                        ),
                    ],
                    style = (
                        width = "49%",
                        float = "right",
                        display = "inline-block",
                    ),
                ),
            ],
            style = (
                borderBottom = "thin lightgrey solid",
                backgroundColor = "rgb(250,250,250)",
                padding = "10px 5px",
            ),
        )
        tag_2 = html_div(
            children = [
                dcc_graph(id = "crossfilter-indicator-scatter"),
                dcc_slider(
                    id = "crossfilter-year-slider",
                    min = minimum(years),
                    max = maximum(years),
                    marks = Dict([Symbol(x) => Symbol(x) for x in years]),
                    value = minimum(years),
                    step = nothing,
                )
            ],
            style = (
                width = "49%",
                display = "inline-block",
            ),
        )
        tag_3 = html_div(
            children = [
                dcc_graph(id = "x-time-series"),
                dcc_graph(id = "y-time-series"),
            ],
            style = (
                width = "49%", display = "inline-block",
            #    float="right"
            ),
        )

        return tag_1, tag_2, tag_3
    end

    callback!(
        app,
        Output("crossfilter-indicator-scatter", "figure"),
        Input("crossfilter-xaxis-column", "value"),
        Input("crossfilter-yaxis-column", "value"),
        Input("crossfilter-xaxis-type", "value"),
        Input("crossfilter-yaxis-type", "value"),
        Input("crossfilter-year-slider", "value"),
    ) do xaxis_column_name, yaxis_column_name, xaxis_type, yaxis_type, year_slider_value
        dff = df[df.Year .== year_slider_value, :] # filtered
        pt = Plot(
            dff[dff[!, Symbol("Indicator Name")] .== xaxis_column_name, :Value],
            dff[dff[!, Symbol("Indicator Name")] .== yaxis_column_name, :Value],
            Layout(
                xaxis_title = xaxis_column_name,
                xaxis_type = xaxis_type == "linear" ? "linear" : "log",
                yaxis_title = yaxis_column_name,
                yaxis_type = yaxis_type == "linear" ? "linear" : "log",
                hovermode = "closest",
                height = 450,
            ),
            kind = "scatter",
            text = dff[
                dff[!, Symbol("Indicator Name")] .== yaxis_column_name,
                Symbol("Country Name"),
            ],
            customdata = dff[
                dff[!, Symbol("Indicator Name")] .== yaxis_column_name,
                Symbol("Country Name"),
            ],
            mode = "markers",
            marker_size = 15,
            marker_opacity = 0.5,
            marker_line_width = 0.5,
            marker_line_color = "white",
        )

        return pt
    end

    callback!(
        app,
        Output("x-time-series", "figure"),
        Input("crossfilter-indicator-scatter", "hoverData"),
        Input("crossfilter-xaxis-column", "value"),
        Input("crossfilter-xaxis-type", "value"),
    ) do hover_data, xaxis_column_name, axis_type
        #@show hover_data
        country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
        dff = df[df[:, Symbol("Country Name")] .== country_name, :]
        dff = dff[dff[:, Symbol("Indicator Name")] .== xaxis_column_name, :]
        title = "<b>$(country_name)</b><br>$(xaxis_column_name)"

        return create_time_series(dff, axis_type, title)
    end

    callback!(
        app,
        Output("y-time-series", "figure"),
        Input("crossfilter-indicator-scatter", "hoverData"),
        Input("crossfilter-yaxis-column", "value"),
        Input("crossfilter-yaxis-type", "value"),
    ) do hover_data, yaxis_column_name, axis_type
        #@show hover_data
        country_name = isnothing(hover_data) ? "" : hover_data.points[1].customdata
        dff = df[df[:, Symbol("Country Name")] .== country_name, :]
        dff = dff[dff[:, Symbol("Indicator Name")] .== yaxis_column_name, :]
        #title = "<b>$(country_name)</b><br>$(yaxis_column_name)"

        return create_time_series(dff, axis_type, yaxis_column_name)
    end

    #return df
    run_server(app, "0.0.0.0", debug=true)
end

"""

Utility function for `dash_update_graphs_on_hover()`.
"""
function create_time_series(df, axis_type, title)
    pt = Plot(
        df[:, :Year],
        df[:, :Value],
        Layout(
            yaxis_type = axis_type == "linear" ? "linear" : "log",
            xaxis_showgrid = false,
            annotations = [attr(
                x = 0,
                y = 0.85,
                xanchor = "left",
                yanchor = "bottom",
                xref = "paper",
                yref = "paper",
                showarrow = false,
                align = "left",
                bgcolor = "rgba(255,255,255,0.5)",
                text = title,
            )],
            height = 225,
        ),
        mode = "lines+markers",
    )

    return pt
end

"""
[Update Graphs on Hover](https://dash-julia.plotly.com/interactive-graphing)
"""
function dash_generic_crossfilter_recipe()
    data = [("Col $(i)" => rand(30)) for i = 1:6] # array of pairs
    db = Dict(data)
    df = DataFrame(db)

    app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

    app.layout = html_div() do
        tag_1 = html_div(
            dcc_graph(id = "graph_1"),
            style = (width = "32%", display = "inline-block"),
        )
        tag_2 = html_div(
            dcc_graph(id = "graph_2"),
            style = (width = "32%", display = "inline-block"),
        )
        tag_3 = html_div(
            dcc_graph(id = "graph_3"),
            style = (width = "32%", display = "inline-block"),
        )

        return tag_1, tag_2, tag_3
    end

    callback!(
        app,
        Output("graph_1", "figure"),
        Output("graph_2", "figure"),
        Output("graph_3", "figure"),
        Input("graph_1", "selectedData"),
        Input("graph_2", "selectedData"),
        Input("graph_3", "selectedData"),
    ) do s1, s2, s3
        selectedpoints = 1:size(df)[1]

        for selected_data in [s1, s2, s3]
            if selected_data != nothing
                # JS is 0-offset?
                selectedpoints = [p[:customdata] for p in selected_data.points]
                @info selectedpoints
            end
        end

        fig_1 = create_figure(df, "Col 1", "Col 2", selectedpoints, s1)
        fig_2 = create_figure(df, "Col 3", "Col 4", selectedpoints, s2)
        fig_3 = create_figure(df, "Col 5", "Col 6", selectedpoints, s3)

        return fig_1, fig_2, fig_3
    end

    run_server(app, "0.0.0.0", debug=true)
end

"""

Utility function for `dash_generic_crossfilter_recipe()`
"""
function create_figure(df, x_col, y_col, selectedpoints, local_selectedpoints)

    # This outcome of this `if` is not used.
    # Get selection bounds
    if local_selectedpoints != nothing
        @info "local_selectedpoint is something"
        # ranges = (x = Any[0.29, 1.09], y = Any[4.59, 23.78])
        ranges = local_selectedpoints[:range]
        @info ranges
        selection_bounds = Dict(
            "x0" => ranges[:x][1],
            "x1" => ranges[:x][2],
            "y0" => ranges[:y][1],
            "x0" => ranges[:y][2],
        )
    else
        @info "local_selectedpoint is nothing"
        selection_bounds = Dict(
            "x0" => minimum(df[:, x_col][1]),
            "x1" => minimum(df[:, x_col][2]),
            "y0" => minimum(df[:, y_col][1]),
            "x0" => minimum(df[:, x_col][2]),
        )
    end

    pt = Plot(
        df,

        # use this group
        Symbol(x_col),
        Symbol(y_col),
        # or this group
        #x = df[:, x_col],
        #y = df[:, y_col],

        mode = "markers+text",
        marker_size = 20,
        text = 1:size(df)[1],
        customdata = 1:size(df)[1],
        selectedpoints = selectedpoints .- 1, # Index in JS use 0-offset
        unselected = (
            marker = (opacity = 0.3, textfont = (color = "rgba(0,0,0,0)"))
        )
    )

    return pt
end

# Use run_Dash.jl to run
#
# https://dash-julia.plotly.com/getting-started
#hello_dash()
#dash_table()
#dash_plotly() # RDatasets does not work
#dash_markdown()
#dash_core()
#dash_basic_callback()
#dash_layout_figure_slider()
#Ans = dash_app_multiple_inputs()
#Ans = dash_app_multiple_outputs()
#Ans = dash_app_chained_callbacks()
#Ans = dash_app_state()
#Ans = dash_intro_visual()
#Ans = dash_update_graphs_on_hover()

nothing


text:
################### spin system.

function evalitpproxysys(qs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysParamsType1{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ
    #κs_β = x.κs_β

    @assert length(d) == length(qs)

    out = zero(Complex{T})

    for i = 1:length(qs)
        r = u_rad - d[i]

        for k = 1:length(qs[i])

            out += qs[i][k](r, κs_λ[i])
        end
    end

    ## slower possibly due to r = u_rad - d[i] being evaluated every time qs is called.
    #out = sum( sum(qs[i][k](u_rad - d[i], κs_λ[i]) for k = 1:length(qs[i])) for i = 1:length(qs) )

    return out
end

function evalitpproxysys(qs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysParamsType1{T}, κs_α::Vector{Vector{T}})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ
    #κs_β = x.κs_β

    @assert length(d) == length(qs)

    out = zero(Complex{T})

    for i = 1:length(qs)
        r = u_rad - d[i]

        for k = 1:length(qs[i])

            out += κs_α[i][k]*qs[i][k](r, κs_λ[i])
        end
    end

    ## slower possibly due to r = u_rad - d[i] being evaluated every time qs is called.
    #out = sum( sum(qs[i][k](u_rad - d[i], κs_λ[i]) for k = 1:length(qs[i])) for i = 1:length(qs) )

    return out
end

function evalitpproxysys(qs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysParamsType2{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ
    #κs_β = x.κs_β

    @assert length(d) == length(qs)

    out = zero(Complex{T})

    #u_rad = 2*π*u
    for i = 1:length(qs)

        for k = 1:length(qs[i])
            r = u_rad - d[i][k]

            #out += qs[i][k](r, κs_λ[i][k], κs_β[i])
            out += qs[i][k](r, κs_λ[i][k])
        end
    end

    return out
end

function evalκitpproxysys(κ_α::Vector{Vector{T}}, qs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysParamsType1{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ
    #κs_β = x.κs_β

    @assert length(d) == length(qs)

    out = zero(Complex{T})

    #u_rad = 2*π*u
    for i = 1:length(qs)
        r = u_rad - d[i]

        for k = 1:length(qs[i])
            #out += κ_α[i][k]*qs[i][k](r, κs_λ[i], κs_β[i])
            out += κ_α[i][k]*qs[i][k](r, κs_λ[i])
        end
    end

    return out
end

function evalκitpproxysys(κ_α::Vector{Vector{T}}, qs::Vector{Vector{Function}},
    u_rad::T, x::SpinSysParamsType2{T})::Complex{T} where T

    d = x.d
    κs_λ = x.κs_λ

    @assert length(d) == length(qs)

    out = zero(Complex{T})

    #u_rad = 2*π*u
    for i = 1:length(qs)

        for k = 1:length(qs[i])
            r = u_rad - d[i][k]

            #out += κ_α[i][k]*qs[i][k](r, κs_λ[i][k], κs_β[i])
            out += κ_α[i][k]*qs[i][k](r, κs_λ[i][k])
        end
    end

    return out
end

###################### singlets.
function evalsinglets(u_rad::T, d::Vector{T}, αs_singlets::Vector{T}, Ωs_singlets,
    βs_singlets, λ0::T, λ_multipliers::Vector{T}) where T <: Real

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for i = 1:length(αs_singlets)
        τ = u_rad - d[i]

        λ = λ0*λ_multipliers[i]
        Ω = Ωs_singlets[i]
        out += αs_singlets[i]*exp(im*βs_singlets[i])/(λ+im*(τ-Ω))
    end
    return out
end

function evalsinglets(u_rad::T, d::Vector{T},
    αs_singlets::Vector{T}, Ωs_singlets,
    βs_singlets, λ0::T, λ_multipliers::Vector{T},
    κ_α_singlets::Vector{T}) where T <: Real

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for i = 1:length(αs_singlets)
        τ = u_rad - d[i]

        λ = λ0*λ_multipliers[i]
        Ω = Ωs_singlets[i]
        out += κ_α_singlets[i]*αs_singlets[i]*exp(im*βs_singlets[i])/(λ+im*(τ-Ω))
    end
    return out
end

###################### front end.
# with proxy.
function evalitpproxymixture(u_rad, As::Vector{SHType{T}},
    Bs::Vector{FIDModelType{T,SST}};
    w::Vector{T} = ones(T, length(As)))::Complex{T} where {T <: Real,SST}

    #u_rad = 2*π*u

    out = zero(Complex{T})

    for n = 1:length(As)
        out += w[n]*evalitpproxycompound(u_rad, As[n], Bs[n])
    end

    return out
end

# with proxy.
function evalitpproxycompound(u_rad, A::SHType{T}, B::FIDModelType{T,SST})::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out_sys = evalitpproxysys(B.qs, u_rad, B.ss_params)

    out_singlets = evalsinglets(u_rad, B.d_singlets, A.αs_singlets, A.Ωs_singlets,
    B.β_singlets, B.λ0, B.κs_λ_singlets)

    return out_sys + out_singlets
end

# with κ-proxy. refactor/remove this later.
function evalitpproxymixture(u_rad, As::Vector{SHType{T}},
    Es::Vector{καFIDModelType{T,SST}};
    w::Vector{T} = ones(T, length(As)))::Complex{T} where {T <: Real,SST}

    #u_rad = 2*π*u

    out = zero(Complex{T})

    for n = 1:length(Es)
        out += w[n]*evalitpproxycompound(u_rad, As[n], Es[n])
    end

    return out
end

# with κ-proxy.
function evalitpproxycompound(u_rad, A::SHType{T}, E::καFIDModelType{T,SST})::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out_sys = evalitpproxysys(E.core.qs, u_rad, E.core.ss_params, E.κs_α)

    out_singlets = evalsinglets(u_rad, E.core.d_singlets, A.αs_singlets, A.Ωs_singlets,
    E.core.β_singlets, E.core.λ0, E.core.κs_λ_singlets, E.κs_α_singlets)

    return out_sys + out_singlets
end

function findsimplecoherences(c2::Vector{T};
    atol = 1e-6)::Vector{Int} where T

    c2a = collect( abs.(c2[i]) for i = 1:length(c2))
    sum_c2a = collect(sum(c2a[i]) for i = 1:length(c2a))
    flags = isapprox.(sum_c2a, 1.0, atol = atol)

    return collect(1:length(c2))[flags]
end


text:
mutable struct GFSimulateState
    trace::GFTrace
    score::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFSimulateState(params::Dict{Symbol,Any})
    GFSimulateState(GFTrace(), 0., AddressVisitor(), params)
end

get_args_change(state::GFSimulateState) = nothing
get_addr_change(state::GFSimulateState, addr) = nothing
set_ret_change!(state::GFSimulateState, value) = begin end

function addr(state::GFSimulateState, dist::Distribution{T}, args, addr, arg_change) where {T}
    visit!(state.visitor, addr)
    retval::T = random(dist, args...)
    score = logpdf(dist, retval, args...)
    call = CallRecord(score, retval, args)
    state.trace = assoc_primitive_call(state.trace, addr, call)
    state.trace.has_choices = true
    state.score += score
    retval
end

function addr(state::GFSimulateState, gen::Generator{T,U}, args, addr, arg_change) where {T,U}
    visit!(state.visitor, addr)
    trace::U = simulate(gen, args)
    call::CallRecord = get_call_record(trace)
    state.trace = assoc_subtrace(state.trace, addr, trace)
    state.trace.has_choices = state.trace.has_choices || has_choices(trace)
    state.score += call.score
    call.retval::T
end

splice(state::GFSimulateState, gf::GenFunction, args::Tuple) = exec(gf, state, args)

function simulate(gen::GenFunction, args)
    state = GFSimulateState(gen.params)
    retval = exec(gen, state, args)
    # TODO add return type annotation for gen functions
    call = CallRecord{Any}(state.score, retval, args)
    state.trace.call = call
    state.trace
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Simpler
using Base.Test

d = Dict("a" => [Dict("x"=>1,"y"=>2)])
@test """Dict("a"=>Dict{String,Int64}[Dict("x"=>1,"y"=>2)])""" == string(d)
@test "{a: [{x: 1, y: 2}]}" == simple(d, 1)


text:
using Base.Test

using Surprise
using Persa

#Testing
ds = Persa.createdummydataset()

holdout = Persa.HoldOut(ds, 0.9)

(ds_train, ds_test) = Persa.get(holdout)

model = Surprise.IRSVD(ds_train)

Persa.train!(model, ds_train)

@assert Persa.aval(model, ds_test).mae >= 0.0


text:
export FiniteOperator

immutable FiniteOperator{AT<:AbstractMatrix,T<:Number,DS,RS} <: Operator{T}
    matrix::AT
    domainspace::DS
    rangespace::RS
end

FiniteOperator(M::AbstractMatrix,ds::Space,rs::Space) =
    FiniteOperator{typeof(M),eltype(M),typeof(ds),typeof(rs)}(M,ds,rs)

FiniteOperator(M::AbstractMatrix) = FiniteOperator(M,ℓ⁰,ℓ⁰)

domainspace(F::FiniteOperator) = F.domainspace
rangespace(F::FiniteOperator) = F.rangespace

function getindex(F::FiniteOperator,k::Integer,j::Integer)
    if k ≤ size(F.matrix,1) && j ≤ size(F.matrix,2)
        F.matrix[k,j]
    elseif k ≤ size(F,1) && j ≤ size(F,2)
        zero(eltype(F))
    else
        throw(BoundsError())
    end
end

function getindex(F::FiniteOperator,k::Integer)
    @assert size(F,1) == 1
    if k ≤ length(F.matrix)
        F.matrix[k]
    else
        zero(eltype(F))
    end
end

function Base.convert{AT<:BandedMatrix,T}(::Type{BandedMatrix},S::SubOperator{T,FiniteOperator{AT,T}})
    kr,jr=parentindexes(S)
    if last(kr[1]) ≤ size(S.matrix,1) &&
        last(jr[2]) ≤ size(S.matrix,2)
        matrix[kr,jr]
    else
        default_copy(S)
    end
end

bandinds(T::FiniteOperator) = (1-size(T.matrix,1),size(T.matrix,2)-1)
bandinds{AT<:BandedMatrix}(T::FiniteOperator{AT}) = bandinds(T.matrix)

Base.maximum(K::FiniteOperator) = maximum(K.matrix)


text:
# This file was generated, do not modify it. # hide
ŷ = inverse_transform(box, ẑ)


text:
"""
    unbroadcast(x,dx)

Bring dx to x's size via unbroadcasting (reduction). This is needed
when defining gradients of multi-argument broadcasting functions where
the arguments and the result may be of different sizes.

"""
function unbroadcast(x, dx)
    if size(x)==size(dx)
        return dx
    elseif isa(value(x),Number)
        return sum(dx)
    elseif isa(value(dx),Number)
        return fill!(similar(value(x)),dx)
    else
        d = []
        for i=1:ndims(dx)
            size(x,i) == size(dx,i) > 1 && continue
            size(x,i) != 1 && throw(DimensionMismatch())
            push!(d,i)
        end
        length(d)==1 && (d=d[1])
        return reshape(sum(dx, dims=d), size(x))
    end
end

# This fixes unbroadcast when x isa Broadcasted
Base.size(x::Base.Broadcast.Broadcasted,i...)=length.(axes(x,i...))

# Broadcasting dimensions:
# y = f(x...) where f is a broadcasting operation.
# size(y) = broadcast_shape(x...)
# ndims(y) = max ndims(x)
# size(y,i) = max size(x,i)
# size(x,i) = 1 or size(y,i) for all x and i<=ndims(x)
# if ndims(x) < ndims(y) the extra dimensions of x are treated as 1


text:
__precompile__(false)

module LogManager

export stringColor
export info
export debug
export warn
#export error

RESET_COLOR = "\x1b[39m\x1b[49m"
FONT_COLOR = RESET_COLOR

function stringColor(s...;color=:normal, save=false)
  global FONT_COLOR
  colorstr=FONT_COLOR

  if color == :yellow colorstr="\x1b[93m"
  elseif color == :darkyellow colorstr="\x1b[33m"
  elseif color == :red colorstr="\x1b[91m"
  elseif color == :darkred colorstr="\x1b[31m"
  elseif color == :cyan colorstr="\x1b[96m"
  elseif color == :darkcyan colorstr="\x1b[36m"
  elseif color == :magenta colorstr="\x1b[95m"
  elseif color == :darkmagenta colorstr="\x1b[35m"
  elseif color == :reset colorstr=RESET_COLOR
  end

  if save FONT_COLOR=colorstr end

  return string(colorstr,s...,FONT_COLOR)
end

setFontColor(s...) = replace(string(s...),RESET_COLOR=>FONT_COLOR)

info(s...) = println(string(stringColor("";color=:yellow,save=true),setFontColor(s...),stringColor("";color=:reset,save=true)))
debug(s...) = println("\x1b[36mDEBUG: ",string(stringColor("";color=:cyan,save=true),setFontColor(s...),stringColor("";color=:reset,save=true)))
warn(s...) = println("\x1b[35mWARNING: ",string(stringColor("";color=:magenta,save=true),setFontColor(s...),stringColor("";color=:reset,save=true)))
error(s...) = println("\x1b[41mERROR\x1b[49m: ",string(stringColor("";color=:red,save=true),setFontColor(s...),stringColor("";color=:reset,save=true)))

end #LogManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
export Polynomial

# Invariant:
# a and x might be empty: meaning it is the zero polynomial
# a does not contain any zeros
# x is increasing in the monomial order (i.e. grlex)
struct Polynomial{C, T} <: AbstractPolynomial{T}
    a::Vector{T}
    x::MonomialVector{C}

    function Polynomial{C, T}(a::Vector{T}, x::MonomialVector{C}) where {C, T}
        length(a) == length(x) || throw(ArgumentError("There should be as many coefficient than monomials"))
        p = new{C, T}(a, x)
        _remove_zeros!(p)
        return p
    end
end

iscomm(::Type{Polynomial{C, T}}) where {C, T} = C

Base.broadcastable(p::Polynomial) = Ref(p)
MA.mutable_copy(p::Polynomial{C, T}) where {C, T} = Polynomial{C, T}(MA.mutable_copy(p.a), MA.mutable_copy(p.x))
Base.copy(p::Polynomial) = MA.mutable_copy(p)
Base.zero(::Type{Polynomial{C, T}}) where {C, T} = Polynomial(T[], MonomialVector{C}())
Base.one(::Type{Polynomial{C, T}}) where {C, T} = Polynomial([one(T)], MonomialVector{C}(PolyVar{C}[], [Int[]]))
Base.zero(p::Polynomial{C, T}) where {C, T} = Polynomial(T[], emptymonovec(copy(_vars(p))))
Base.one(p::Polynomial{C, T}) where {C, T} = Polynomial([one(T)], MonomialVector(copy(_vars(p)), [zeros(Int, nvariables(p))]))

Polynomial{C, T}(a::AbstractVector, x::MonomialVector) where {C, T} = Polynomial{C, T}(Vector{T}(a), x)
Polynomial{C, T}(a::AbstractVector, X::DMonoVec) where {C, T} = Polynomial{C, T}(monovec(a, X)...)
Polynomial{C}(a::Vector{T}, x) where {C, T} = Polynomial{C, T}(a, x)
Polynomial(af::Union{Function, Vector}, x::DMonoVec{C}) where {C} = Polynomial{C}(af, x)

# TODO Remove with MP v0.2.8
Polynomial{C, T}(p::Polynomial{C, T}) where {C, T} = p

Base.convert(::Type{Polynomial{C, T}}, p::Polynomial{C, T}) where {C, T} = p
function Base.convert(::Type{Polynomial{C, T}},
                      p::Polynomial{C, S}) where {C, S, T}
    return Polynomial{C}(convert(Vector{T}, p.a), p.x)
end
#function convert(::Type{Polynomial{C, T}},
#                 p::AbstractPolynomialLike) where {C, T}
#    return convert(Polynomial{C, T}, polynomial(p, T))
#end
function Base.convert(::Type{Polynomial{C, T}}, t::Term{C}) where {C, T}
    return Polynomial{C, T}(T[t.α], [t.x])
end
function Base.convert(::Type{Polynomial{C, T}}, m::DMonomialLike{C}) where {C, T}
    return Polynomial(convert(Term{C, T}, m))
end
function MP.convertconstant(::Type{Polynomial{C, T}}, α) where {C, T}
    return Polynomial(convert(Term{C, T}, α))
end

Polynomial{C}(p::Union{Polynomial{C}, Term{C}, Monomial{C}, PolyVar{C}}) where {C} = Polynomial(p)
Polynomial{C}(α) where {C} = Polynomial(Term{C}(α))

Polynomial(p::Polynomial) = p
Polynomial(t::Term{C, T}) where {C, T} = Polynomial{C, T}([t.α], [t.x])
Polynomial(x::Union{PolyVar{C}, Monomial{C}}) where {C} = Polynomial(Term{C}(x))

#Base.convert(::Type{TermContainer{C, T}}, p::Polynomial{C}) where {C, T} = Polynomial{C, T}(p)

function Polynomial{C, T}(f::Function, x::MonomialVector{C}) where {C, T}
    a = T[f(i) for i in 1:length(x)]
    Polynomial{C, T}(a, x)
end
function Polynomial{C, T}(f::Function, x::AbstractVector) where {C, T}
    σ, X = sortmonovec(x)
    a = T[f(i) for i in σ]
    Polynomial{C, T}(a, X)
end
Polynomial{C}(f::Function, x) where {C} = Polynomial{C, Base.promote_op(f, Int)}(f, x)

#Base.convert(::Type{PolyType{C}}, p::TermContainer{C}) where {C} = p

# needed to build [p Q; Q p] where p is a polynomial and Q is a matpolynomial in Julia v0.5
#Base.convert(::Type{TermType{C}}, p::TermContainer{C}) where {C} = p
#Base.convert(::Type{TermType{C, T}}, p::TermContainer{C, T}) where {C, T} = p

Base.length(p::Polynomial) = length(p.a)
Base.isempty(p::Polynomial) = isempty(p.a)
Base.iterate(p::Polynomial) = isempty(p) ? nothing : (p[1], 1)
function Base.iterate(p::Polynomial, state::Int)
    state < length(p) ? (p[state+1], state+1) : nothing
end
#eltype(::Type{Polynomial{C, T}}) where {C, T} = T
Base.getindex(p::Polynomial, I::Int) = Term(p.a[I[1]], p.x[I[1]])

#Base.transpose(p::Polynomial) = Polynomial(map(transpose, p.a), p.x) # FIXME invalid age range update

struct TermIterator{C, T} <: AbstractVector{Term{C, T}}
    p::Polynomial{C, T}
end
Base.firstindex(p::TermIterator) = firstindex(p.p.a)
Base.lastindex(p::TermIterator) = lastindex(p.p.a)
Base.length(p::TermIterator) = length(p.p.a)
Base.size(p::TermIterator) = (length(p),)
Base.isempty(p::TermIterator) = isempty(p.p.a)
Base.iterate(p::TermIterator) = isempty(p) ? nothing : (p[1], 1)
function Base.iterate(p::TermIterator, state::Int)
    state < length(p) ? (p[state+1], state+1) : nothing
end

Base.getindex(p::TermIterator, I::Int) = Term(p.p.a[I[1]], p.p.x[I[1]])

MP.terms(p::Polynomial) = TermIterator(p)
MP.coefficients(p::Polynomial) = p.a
MP.monomials(p::Polynomial) = p.x
_vars(p::Polynomial) = _vars(p.x)

MP.extdegree(p::Polynomial) = extdegree(p.x)
MP.mindegree(p::Polynomial) = mindegree(p.x)
MP.maxdegree(p::Polynomial) = maxdegree(p.x)

MP.leadingcoefficient(p::Polynomial{C, T}) where {C, T} = iszero(p) ? zero(T) : first(p.a)
MP.leadingmonomial(p::Polynomial) = iszero(p) ? constantmonomial(p) : first(p.x)
MP.leadingterm(p::Polynomial) = iszero(p) ? zeroterm(p) : first(terms(p))

function MP.removeleadingterm(p::Polynomial)
    Polynomial(p.a[2:end], p.x[2:end])
end
function MA.mutable_operate!(::typeof(MP.removeleadingterm), p::Polynomial)
    deleteat!(p.a, 1)
    deleteat!(p.x, 1)
    return p
end
function MP.removemonomials(p::Polynomial, x::MonomialVector)
    # use the fact that monomials are sorted to do this O(n) instead of O(n^2)
    j = 1
    I = Int[]
    for (i,t) in enumerate(p)
        while j <= length(x) && x[j] > t.x
            j += 1
        end
        if j > length(x) || x[j] != t.x
            push!(I, i)
        end
    end
    Polynomial(p.a[I], p.x[I])
end
MP.removemonomials(p::Polynomial, x::Vector) = removemonomials(p, MonomialVector(x))

function _remove_zeros!(p::Polynomial)
    zeroidx = Int[]
    for (i, α) in enumerate(p.a)
        if iszero(α)
            push!(zeroidx, i)
        end
    end
    if !isempty(zeroidx)
        deleteat!(p.a, zeroidx)
        deleteat!(p.x.Z, zeroidx)
    end
end

function removedups_to!(a::Vector{T}, Z::Vector{Vector{Int}},
                        adup::Vector{T}, Zdup::Vector{Vector{Int}}) where T
    σ = sortperm(Zdup, rev=true, lt=grlex)
    i = 0
    j = 1
    while j <= length(adup)
        k = σ[j]
        if j == 1 || Zdup[k] != Zdup[σ[j-1]]
            push!(Z, Zdup[k])
            push!(a, adup[k])
            i += 1
        else
            a[i] = MA.operate!(+, a[i], adup[k])
        end
        j += 1
    end
end
function polynomialclean(vars::Vector{PolyVar{C}}, adup::Vector{T}, Zdup::Vector{Vector{Int}}) where {C, T}
    Z = Vector{Int}[]
    a = T[]
    removedups_to!(a, Z, adup, Zdup)
    Polynomial{C, T}(a, MonomialVector{C}(vars, Z))
end
function polynomialclean_to!(p::Polynomial{C, T}, vars::Vector{PolyVar{C}}, adup::Vector{T}, Zdup::Vector{Vector{Int}}) where {C, T}
    Future.copy!(p.x.vars, vars)
    empty!(p.a)
    empty!(p.x.Z)
    removedups_to!(p.a, p.x.Z, adup, Zdup)
    _remove_zeros!(p)
    return p
end

MP.polynomial!(a::Vector, x::DMonoVec, ::MP.ListState) = Polynomial(a, x)
MP.polynomial(a::AbstractVector, x::DMonoVec, s::MP.ListState) = MP.polynomial!(collect(a), MA.mutable_copy(x), s)

#MP.polynomial(f::Function, x::AbstractVector) = Polynomial(f, x)
#MP.polynomial(ts::AbstractVector{Term{C, T}}) where {C, T} = Polynomial(coefficient.(ts), monomial.(ts)) # FIXME invalid age range update

# i < j
function trimap(i, j, n)
    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1
end
MP.polynomial(Q::AbstractMatrix{T}, mv::MonomialVector) where T = MP.polynomial(Q, mv, Base.promote_op(+, T, T))
function MP.polynomial(Q::AbstractMatrix, mv::MonomialVector{C}, ::Type{T}) where {C, T}
    if isempty(Q)
        zero(Polynomial{C, T})
    else
        n = length(mv)
        if C
            N = trimap(n, n, n)
            Z = Vector{Vector{Int}}(undef, N)
            a = Vector{T}(undef, N)
            for i in 1:n
                for j in i:n
                    k = trimap(i, j, n)
                    Z[k] = mv.Z[i] + mv.Z[j]
                    if i == j
                        a[k] = Q[i, j]
                    else
                        a[k] = Q[i, j] + Q[j, i]
                    end
                end
            end
            v = _vars(mv)
        else
            N = n^2
            x = Vector{Monomial{C}}(undef, N)
            a = Vector{T}(undef, N)
            offset = 0
            for i in 1:n
                # for j in 1:n wouldn't be cache friendly for Q
                for j in i:n
                    k = trimap(i, j, n)
                    q = Q[i, j]
                    x[offset+k] = mv[i] * mv[j]
                    a[offset+k] = q
                    if i != j
                        offset += 1
                        x[offset+k] = mv[j] * mv[i]
                        a[offset+k] = q
                    end
                end
            end
            a, X = monovec(a, x)
            v = _vars(X)
            Z = X.Z
        end
        polynomialclean(v, a, Z)
    end
end

function MA.mutable_operate!(::typeof(zero), p::Polynomial)
    empty!(p.a)
    empty!(p.x.Z)
    return p
end
function MA.mutable_operate!(::typeof(one), p::Polynomial{C, T}) where {C, T}
    if isempty(p.a)
        push!(p.a, one(T))
        push!(p.x.Z, zeros(Int, length(p.x.vars)))
    else
        resize!(p.a, 1)
        MA.mutable_operate!(one, p.a[1])
        resize!(p.x.Z, 1)
        MA.mutable_operate!(zero, p.x.Z[1])
    end
    return p
end


text:
import ImageView
using Images, Gtk.ShortNames
using TestImages

gui = imshow_gui((400, 300), (2, 2); name="canvasgrid")
c = gui["canvas"]
imshow(c[1,1], testimage("lighthouse"))
imshow(c[1,2], testimage("mountainstream"))
imshow(c[2,1], testimage("moonsurface"))
imshow(c[2,2], testimage("mandrill"))
Gtk.showall(gui["window"])
sleep(0.01)


text:
# RunAllSimulations
println("Spinning up model code ... ")
include("../src/new_SimulationFunctions.jl")

# Establish simulation parameters
TSTART = 0.0;
Ts = 1.0; # 0.1
TSTOP = 250.0; # 100 ; 500

# Load data dictionary
data_dictionary = DataFile(TSTART,TSTOP,Ts)

 # # # #
#println(data_dictionary["CONTROL_PARAMETER_ARRAY"])

# update control parameter array here for now
GAIN_MASTER = 0.3
num_regs = 19;
control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"]
for i in 1:num_regs
  control_parameter_array[i,1]=GAIN_MASTER # Gain for transcriptional activity
end
# TGFB3_inhibition: TGFB3_Ab inhibition; TGFB3_binding
control_parameter_array[1,1]=30.0 # GAIN # _tgfb3_Ab
# NFATc_inhibition: vivit inhibition; Calcineurin_activation_NFATc
control_parameter_array[2,1]=3.0 # GAIN # NFATc_inhibition: vivit inhibition; Calcineurin_activation_NFATc
# GSK3_P_inhibition_Binding_APC_AXIN_forw: GSK3_P inhibition; Binding_APC_AXIN_forw
control_parameter_array[3,1]=0.5 # GAIN # GSK3_P_inhibition_Binding_APC_AXIN_forw
control_parameter_array[3,2]=1.0 # ORDER# GSK3_P_inhibition_Binding_APC_AXIN_forw

# Active_TCF4_inhibition_Ecad: Ecadherin_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[4,1]=0.1 # GAIN # Active_TCF4_inhibition_Ecad
control_parameter_array[4,1]=0.9 # GAIN # Active_TCF4_inhibition_Ecad
# Active_TCF4_inhibition_APC: APC_AXIN_Bcatenin inhibition; TCF4_Bcatenin_complex_formation
control_parameter_array[5,1]=0.2 # GAIN # Active_TCF4_inhibition_APC
control_parameter_array[5,2]=1.0 # ORDER # Active_TCF4_inhibition_APC

# activation_Ecadherin_by_Active_NFATc: Active_NFATc induction induction_gene_Ecadherin
control_parameter_array[6,1]=0.3 # GAIN # activation_Ecadherin_by_Active_NFATc
# repression_Ecadherin_by_pSmad: pSmad repression induction_gene_Ecadherin
control_parameter_array[7,1]=5.0 # GAIN # REP_ECAD_PSMAD
control_parameter_array[7,2]=2.0 # ORDER # REP_ECAD_PSMAD
# repression_Ecadherin_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_Ecadherin
control_parameter_array[8,1]=0.5 # GAIN_SNAIL_SLUG_ECAD #2.0
# repression_Ecadherin_by_Active_LEF1: Active_LEF1 repression induction_gene_Ecadherin
control_parameter_array[9,1]=0.6# GAIN_ACTIVE_LEF1_ECAD #2.0 0.5
control_parameter_array[9,2]=1.0# ORDER_ACTIVE_LEF1_ECAD
#  control_parameter_array[6,1]=0.3 # GAIN_BCAT_TGFB3
# activation_Ecadherin_by_ecad_virus: ecad_virus induction induction_gene_Ecadherin
control_parameter_array[10,1]=5.0 # ecad_virus ***
# activation_TGFB3_by_TCF4: Active_TCF4 induction induction_gene_TGFB3
control_parameter_array[11,1]=0.15 # GAIN_TCF4_TGFB3
control_parameter_array[11,2]=1.5 # GAIN_TCF4_TGFB3

# repression_TGFB3_by_YREG1: YREG1 repression induction_gene_TGFB3
control_parameter_array[12,1]=1.0 # YREG_TGFB3 #1.0
control_parameter_array[12,2]=1.0 # YREG_TGFB3
# repression_YREG1_by_SNAIL_SLUG: SNAIL_SLUG repression induction_gene_YREG1
control_parameter_array[13,1]=0.1 # GAIN_SNAIL_YREG1 # 0.1
control_parameter_array[13,2]=2.0 # ORDER_SNAIL_YREG1 # 2.0
# activation_LEF1_by_pSmad: pSmad induction induction_gene_LEF1
control_parameter_array[14,1]=0.8# GAIN # _PSMAD_LEF1 # 0.5 # 0.8
control_parameter_array[14,2]=2.0# ORDER# _PSMAD_LEF1       # 1.0, 2.0
# repression_LEF1_by_YREG1: YREG1 repression induction_gene_LEF1
control_parameter_array[15,1]=5.0 # GAIN YREG1_LEF1 # 0.1
control_parameter_array[15,2]=0.5 # ORDER YREG1_LEF1 # 2.0
# activation_vimentin_by_Active_LEF1: Active_LEF1 induction induction_gene_vimentin
control_parameter_array[16,1]=0.3# GAIN_LEF1_VIM
control_parameter_array[16,2]=2.0# ORDER_LEF1_VIM
# activation_SNAIL_SLUG_by_AP1_SP1_P: AP1_SP1_P induction induction_gene_SNAIL_SLUG
control_parameter_array[18,1]=1.0# SS by AP1 SP1
# activation_SNAIL_SLUG_by_ss_virus: ss_virus induction induction_gene_SNAIL_SLUG
control_parameter_array[19,1]=5.0# SS_VIRUS ***

# # # #

data_dictionary["CONTROL_PARAMETER_ARRAY"]=control_parameter_array
#println(data_dictionary["CONTROL_PARAMETER_ARRAY"])

println("Running to steady state ... ")

# # # STEADY STATE
#@show data_dictionary["RATE_CONSTANT_ARRAY"][141] # here

println(now())
# WARNING data_dicionary will be edited within EstimateSteadyState unless copy/deepcopy used!
XSS = EstimateSteadyState(deepcopy(data_dictionary))
println(now())
initial_condition_array = vec(XSS);
#println(initial_condition_array) ###
#@show data_dictionary["RATE_CONSTANT_ARRAY"][141]

writedlm("../output/xss.dat", XSS)

data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
xss_data_dict = deepcopy(data_dictionary)

# # # # TGFB SIMULATION
#
# println("Running TGFB simulation ... ")
#
# # not required for first simulation
# #data_dictionary= deepcopy(xss_data_dict) # reinitialize dictionary with steady state dictionary
#
# (TSIM,X,ICA) =  SolveModelAddTGFB(TSTART,TSTOP,Ts,data_dictionary)
#
# writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/TGFbSimulation/TSIM.dat", TSIM)
# writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/TGFbSimulation/X.dat", X)

# # # VEGF SIMULATION
## WARNING, SIMULATION DATAFILE MAY BE LEAKING
# println("Running VEGF simulation ... ")
#
# data_dictionary = deepcopy(xss_data_dict) # reinitialize dictionary with steady state dictionary
#
# (TSIM,X,ICA) =  SolveModelAddVEGF(TSTART,TSTOP,Ts,data_dictionary)
#
# writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/VEGFSimulation/TSIM.dat", TSIM)
# writedlm("/Users/dbassen/work/varner_gould_emt_model/emt_grn_bassen_model/Results/VEGFSimulation/X.dat", X)
#
#@show xss_data_dict["RATE_CONSTANT_ARRAY"][141]

simulations_name = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_mapkin.jl","edit_file_ss.jl","edit_file_ssecad.jl","edit_file_ss_tgfb3ab.jl","edit_file_ss_dnlef1.jl","edit_file_tgfb3.jl","edit_file_dnsmad.jl","edit_file_tgfb3_dnlef1.jl","edit_file_tgfb_tgfb3ab.jl"]
#simulations_name = ["edit_file_tgfb.jl","edit_file_vegf.jl","edit_file_vegf_vi.jl","edit_file_tgfb_vi.jl","edit_file_tgfb_vegf.jl","edit_file_vegf_vi.jl","edit_file_tgfb_vegf_vi.jl","edit_file_yreg_ox.jl","edit_file_yreg_ko.jl","edit_file_yreg_ox_tgfb.jl","edit_file_yreg_ko_tgfb.jl"] # Robustness sims
include("spin_up_edits.jl")
#include("robustness_spin_up_edits.jl")
sim_data = [] # Keep data in RAM
i_simdat = 1
for sim in simulations
  name = simulations_name[i_simdat]
  msg = string("Running ",name," simulation ... ")
  println(msg)
  data_dictionary_active = deepcopy(xss_data_dict)
  (TSIM,X) = SolveModelGeneric_ram(TSTART,TSTOP,Ts,data_dictionary_active,sim)

  outdir="../output/Fig2Simulations/"
  tpath= string(outdir,"TSIM",name[10:end-3],".dat")
  xpath= string(outdir,"X",name[10:end-3],".dat")
  push!(sim_data,X)
  writedlm(tpath, TSIM)
  writedlm(xpath, X)
  i_simdat = i_simdat + 1
end


text:
module Crux
    CRUX_WARNINGS = true

    set_crux_warnings(val::Bool) = CRUX_WARNINGS = val

    using Random
    using Distributions
    using POMDPs
    using POMDPModelTools:render
    using Parameters
    using TensorBoardLogger
    using Flux
    using Zygote
    import Zygote: ignore
    using Flux.Optimise: train!
    using CUDA
    using LinearAlgebra
    using ValueHistories
    using DataStructures
    using Plots
    using ColorSchemes
    import Images: save
    using Statistics
    using Base.Iterators: partition

    extra_functions = Dict()
    function set_function(key, val)
        extra_functions[key] = val
    end

    export device, cpucall, gpucall, mdcall, bslice
    include("devices.jl")

    export AbstractSpace, DiscreteSpace, ContinuousSpace, type, dim, state_space
    include("spaces.jl")

    export ObjectCategorical, ConstantLayer, whiten, to2D, tovec,
           LinearDecaySchedule, MultitaskDecaySchedule
    include("utils.jl")

    export MinHeap, inverse_query, mdp_data, PriorityParams, ExperienceBuffer, buffer_like, minibatch,
           clear!, trim!, isprioritized, dim, episodes, update_priorities!, uniform_sample!,
           prioritized_sample!, capacity, normalize!, extra_columns, get_last_N_indices
    include("experience_buffer.jl")

    export TrainingParams, batch_train!
    include("training.jl")

    export PolicyParams, NetworkPolicy, polyak_average!, ContinuousNetwork, DiscreteNetwork,
           DoubleNetwork, ActorCritic, GaussianPolicy, SquashedGaussianPolicy,
           DistributionPolicy, MixedPolicy, ϵGreedyPolicy,
           GaussianNoiseExplorationPolicy, FirstExplorePolicy,
           entropy, logpdf, action_space, exploration, layers, actor, critic
    include("policies.jl")

    export Sampler, initial_observation, terminate_episode!, step!, steps!,
           episodes!, fillto!, metric_by_key, metrics_by_key, undiscounted_return,
           discounted_return, failure, fill_gae!, fill_returns!, trim!
    include("sampler.jl")

    export elapsed, LoggerParams, aggregate_info, log_performance,
           log_discounted_return, log_undiscounted_return, log_failure,
           log_exploration, log_metric_by_key, log_metrics_by_key, log_validation_error
    include("logging.jl")

    export smooth, readtb, plot_learning, episode_frames, gif, percentile,
           find_crossing, plot_jumpstart, directories, plot_peak_performance,
           plot_learning, plot_cumulative_rewards, plot_steps_to_threshold
    include("analysis.jl")

    export DenseSN, ConvSN
    include("extras/spectral_normalization.jl")

    export GAN_BCELoss, GAN_LSLoss, GAN_HingeLoss, GAN_WLossGP, GAN_WLoss, Lᴰ, Lᴳ
    include("extras/gans.jl")

    export gradient_penalty
    include("extras/gradient_penalty.jl")

    export DiagonalFisherRegularizer, add_fisher_information_diagonal!, update_fisher!
    include("extras/fisher_information.jl")

    export BatchRegularizer, value_regularization, action_regularization, action_value_regularization
    include("extras/batch_regularization.jl")

    export OrthogonalRegularizer
    include("extras/orthogonal_regularization.jl")

    export log_multitask_performances!, continual_learning
    include("extras/multitask_learning.jl")

    export OnPolicySolver, OffPolicySolver, BatchSolver
    include("model_free/on_policy.jl")
    include("model_free/off_policy.jl")
    include("model_free/batch.jl")

    export REINFORCE, A2C, PPO, LagrangePPO, DQN, DDPG, TD3, SAC
    include("model_free/rl/reinforce.jl")
    include("model_free/rl/a2c.jl")
    include("model_free/rl/ppo.jl")
    include("model_free/rl/dqn.jl")
    include("model_free/rl/ddpg.jl")
    include("model_free/rl/td3.jl")
    include("model_free/rl/sac.jl")

    export TIER, LatentConditionedNetwork
    include("model_free/rl/tier.jl")

    export OnPolicyGAIL, OffPolicyGAIL, BC, AdVIL, SQIL, AdRIL, ASAF
    export mse_action_loss, logpdf_bc_loss, mse_value_loss
    include("model_free/il/bc.jl")
    include("model_free/il/AdVIL.jl")
    include("model_free/il/on_policy_gail.jl")
    include("model_free/il/off_policy_gail.jl")
    include("model_free/il/asaf.jl")
    include("model_free/il/sqil.jl")
    include("model_free/il/AdRIL.jl")

    export CQL, BatchSAC
    include("model_free/batch/cql.jl")
    include("model_free/batch/sac.jl")

    export AdversarialOffPolicySolver, RARL, RARL_DQN, RARL_TD3, ISARL, ISARL_DQN, ISARL_TD3
    include("model_free/adversarial/adv_off_policy.jl")
    include("model_free/adversarial/rarl.jl")
    include("model_free/adversarial/isarl.jl")
end


text:
using Processing
@Processing.version "2DOpenGL"
@Processing.load

coordSystem(-1, 1, 1, -1) # user-defined coordinate system (x-axis min = -1,
						  # x-axis max = 1, y-axis min = -1, y-axis max = 1,
						  # xcent = 0, ycent = 0)
						  #
						  # I can just think easier this way...
						  #
						  # default coordinate system is in units of pixels with
						  # top-left corner of drawing surface labeled (0,0),
						  # just like Processing's default coordinate system
smooth()

ellipse(0,0,0.3,0.5)
ellipse(0,0,0.2,0.2)
triangle(.3, .75, .58, .20, .86, .75)
triangle(.3, .6, .58, .40, .86, .7)
point(0, 0)
fill(1,0,0,1)
quad(-.3, -.75, -.58, -.20, -.86, -.75, -.2, -.4)
noFill()
rect(-.6, -.6, .2, .5)
line(.5, -.4, .7, -.5)
line(.7, -.4, .5, -.5)
strokeWeight(1)
stroke(0.9, 0, 0, 1)
arc(-.6, .6, .3, 0.4pi, pi)
animate()

x = 0
while true
	fill(sin(x)/2+0.5,cos(x)/2+0.5,(sin(x)*cos(x))/2+0.5,1)
	ellipse(0,0,0.2,0.2)
	point(0, 0)
	animate()
	x = x + 0.01
end


text:
using Test: @test, @testset
import Pastebin

@testset "pastebin" begin
    @test Pastebin.content(Pastebin.getraw("Lt9PtFc4")) == "Hola mundo"
    @test Pastebin.content(Pastebin.getraw("gdQAnG2f")) == "Batido de mamey"
    @test Pastebin.success(Pastebin.getraw("gdQAnG2f"))
    @test Pastebin.parsekey("https://pastebin.com/gdQAnG2f") == "gdQAnG2f"

    client = Pastebin.Client(ENV["PASTEBIN_DEVKEY"])
    data = "hola $(randn())"
    response = Pastebin.paste(client, "A pastebin test", data; expire="1H")
    key = Pastebin.parsekey(response)
    @test Pastebin.content(Pastebin.getraw(key)) == data
end