text:
export kmeans

using mlpack._Internal.io

import mlpack_jll
const kmeansLibrary = mlpack_jll.libmlpack_julia_kmeans

# Call the C binding of the mlpack kmeans binding.
function kmeans_mlpackMain()
  success = ccall((:kmeans, kmeansLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module kmeans_internal
  import ..kmeansLibrary

end # module

"""
    kmeans(clusters, input; [algorithm, allow_empty_clusters, in_place, initial_centroids, kill_empty_clusters, labels_only, max_iterations, percentage, refined_start, samplings, seed, verbose])

This program performs K-Means clustering on the given dataset.  It can return
the learned cluster assignments, and the centroids of the clusters.  Empty
clusters are not allowed by default; when a cluster becomes empty, the point
furthest from the centroid of the cluster with maximum variance is taken to fill
that cluster.

Optionally, the Bradley and Fayyad approach ("Refining initial points for
k-means clustering", 1998) can be used to select initial points by specifying
the `refined_start` parameter.  This approach works by taking random samplings
of the dataset; to specify the number of samplings, the `samplings` parameter is
used, and to specify the percentage of the dataset to be used in each sample,
the `percentage` parameter is used (it should be a value between 0.0 and 1.0).

There are several options available for the algorithm used for each Lloyd
iteration, specified with the `algorithm`  option.  The standard O(kN) approach
can be used ('naive').  Other options include the Pelleg-Moore tree-based
algorithm ('pelleg-moore'), Elkan's triangle-inequality based algorithm
('elkan'), Hamerly's modification to Elkan's algorithm ('hamerly'), the
dual-tree k-means algorithm ('dualtree'), and the dual-tree k-means algorithm
using the cover tree ('dualtree-covertree').

The behavior for when an empty cluster is encountered can be modified with the
`allow_empty_clusters` option.  When this option is specified and there is a
cluster owning no points at the end of an iteration, that cluster's centroid
will simply remain in its position from the previous iteration. If the
`kill_empty_clusters` option is specified, then when a cluster owns no points at
the end of an iteration, the cluster centroid is simply filled with DBL_MAX,
killing it and effectively reducing k for the rest of the computation.  Note
that the default option when neither empty cluster option is specified can be
time-consuming to calculate; therefore, specifying either of these parameters
will often accelerate runtime.

Initial clustering assignments may be specified using the `initial_centroids`
parameter, and the maximum number of iterations may be specified with the
`max_iterations` parameter.

As an example, to use Hamerly's algorithm to perform k-means clustering with
k=10 on the dataset `data`, saving the centroids to `centroids` and the
assignments for each point to `assignments`, the following command could be
used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> centroids, assignments = kmeans(10, data)
```

To run k-means on that same dataset with initial centroids specified in
`initial` with a maximum of 500 iterations, storing the output centroids in
`final` the following command may be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> initial = CSV.read("initial.csv")
julia> final, _ = kmeans(10, data; initial_centroids=initial,
            max_iterations=500)
```

# Arguments

 - `clusters::Int`: Number of clusters to find (0 autodetects from initial
      centroids).
 - `input::Array{Float64, 2}`: Input dataset to perform clustering on.
 - `algorithm::String`: Algorithm to use for the Lloyd iteration ('naive',
      'pelleg-moore', 'elkan', 'hamerly', 'dualtree', or 'dualtree-covertree').
      Default value `naive`.

 - `allow_empty_clusters::Bool`: Allow empty clusters to be persist.
      Default value `false`.

 - `in_place::Bool`: If specified, a column containing the learned cluster
      assignments will be added to the input dataset file.  In this case,
      --output_file is overridden. (Do not use in Python.)  Default value
      `false`.

 - `initial_centroids::Array{Float64, 2}`: Start with the specified
      initial centroids.
 - `kill_empty_clusters::Bool`: Remove empty clusters when they occur.
      Default value `false`.

 - `labels_only::Bool`: Only output labels into output file.  Default
      value `false`.

 - `max_iterations::Int`: Maximum number of iterations before k-means
      terminates.  Default value `1000`.

 - `percentage::Float64`: Percentage of dataset to use for each refined
      start sampling (use when --refined_start is specified).  Default value
      `0.02`.

 - `refined_start::Bool`: Use the refined initial point strategy by
      Bradley and Fayyad to choose initial points.  Default value `false`.

 - `samplings::Int`: Number of samplings to perform for refined start (use
      when --refined_start is specified).  Default value `100`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `centroid::Array{Float64, 2}`: If specified, the centroids of each
      cluster will  be written to the given file.
 - `output::Array{Float64, 2}`: Matrix to store output labels or labeled
      data to.

"""
function kmeans(clusters::Int,
                input;
                algorithm::Union{String, Missing} = missing,
                allow_empty_clusters::Union{Bool, Missing} = missing,
                in_place::Union{Bool, Missing} = missing,
                initial_centroids = missing,
                kill_empty_clusters::Union{Bool, Missing} = missing,
                labels_only::Union{Bool, Missing} = missing,
                max_iterations::Union{Int, Missing} = missing,
                percentage::Union{Float64, Missing} = missing,
                refined_start::Union{Bool, Missing} = missing,
                samplings::Union{Int, Missing} = missing,
                seed::Union{Int, Missing} = missing,
                verbose::Union{Bool, Missing} = missing,
                points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, kmeansLibrary), Nothing, ());

  IORestoreSettings("K-Means Clustering")

  # Process each input argument before calling mlpackMain().
  IOSetParam("clusters", clusters)
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(algorithm)
    IOSetParam("algorithm", convert(String, algorithm))
  end
  if !ismissing(allow_empty_clusters)
    IOSetParam("allow_empty_clusters", convert(Bool, allow_empty_clusters))
  end
  if !ismissing(in_place)
    IOSetParam("in_place", convert(Bool, in_place))
  end
  if !ismissing(initial_centroids)
    IOSetParamMat("initial_centroids", initial_centroids, points_are_rows)
  end
  if !ismissing(kill_empty_clusters)
    IOSetParam("kill_empty_clusters", convert(Bool, kill_empty_clusters))
  end
  if !ismissing(labels_only)
    IOSetParam("labels_only", convert(Bool, labels_only))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(percentage)
    IOSetParam("percentage", convert(Float64, percentage))
  end
  if !ismissing(refined_start)
    IOSetParam("refined_start", convert(Bool, refined_start))
  end
  if !ismissing(samplings)
    IOSetParam("samplings", convert(Int, samplings))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("centroid")
  IOSetPassed("output")
  # Call the program.
  kmeans_mlpackMain()

  return IOGetParamMat("centroid", points_are_rows),
         IOGetParamMat("output", points_are_rows)
end


text:
########## Theese methods belong to configuration folder ############

using DemeNet: Certificate, Intent, Signer, Deme
import DemeNet: serialize, deserialize

using PeaceVote: BraidChain
using PeaceVote.BraidChains: getrecord, record!, readbootrecord, writebootrecord, Ledger

# Directory could be obtained with a method.

function binary(x)
    io = IOBuffer()
    serialize(io,x)
    return take!(io)
end

# A dity hack
# function serialize(ledger::Ledger,config::Certificate{PeaceFounderConfig})
#     fname = "PeaceFounder.toml"
#     bytes = binary(config)
#     record!(ledger,fname,bytes,1)
#     write(ledger.dir * "/" * fname, bytes)
# end

serialize(ledger::Ledger,config::Certificate{PeaceFounderConfig}) = writebootrecord(ledger,config)
#     fname = "PeaceFounder.toml"
#     bytes = binary(config)
#     record!(ledger,fname,bytes,1)
#     write(ledger.dir * "/" * fname, bytes)
# end

#serialize(ledger,config,"PeaceFounder.toml")

#serialize(ledger::Ledger,config::PeaceFounderConfig) = serialize(ledger,config,"PeaceFounder.toml")

function serialize(deme::Deme,config::PeaceFounderConfig) # #writebootrecord(ledger,config)
    uuid = deme.spec.uuid
    fname = datadir(uuid) * "/PeaceFounder.toml"
    mkpath(dirname(fname))

    io = IOBuffer()
    serialize(io,config)
    bytes = take!(io)

    write(fname,bytes)
    # open(fname, "w") do io
    #     return serialize(io,config)
    # end
end

function readbytes(fname::AbstractString)
    file = open(fname,"r")
    data = UInt8[]
    while !eof(file)
        push!(data,read(file,UInt8))
    end
    close(file)
    return data
end

### Did this one worked?
function deserialize(deme::Deme,config::Type{PeaceFounderConfig})
    fname = datadir(deme.spec.uuid) * "/PeaceFounder.toml"

    bytes = readbytes(fname)
    return deserialize(IOBuffer(bytes),PeaceFounderConfig)

    # open(fname, "r") do io
    #     return cert = deserialize(io,PeaceFounderConfig)
    # end
end

function deserialize(deme::Deme,config::Type{Certificate{PeaceFounderConfig}})
    fname = datadir(deme.spec.uuid) * "/PeaceFounder.toml"

    bytes = readbytes(fname)
    return deserialize(IOBuffer(bytes),Certificate{PeaceFounderConfig})
    # open(fname, "r") do io
    #     @show readavailable(io)

    #     return cert = deserialize(io,Certificate{PeaceFounderConfig})
    # end
end

### I could unite theese two methods into one
deserialize(ledger::Ledger,type::Type{Certificate{PeaceFounderConfig}}) = readbootrecord(ledger,type)
#     rec =
#     #rec = getrecord(ledger,"PeaceFounder.toml") # I did not update the ledger I guess
#     #println(String(copy(rec.data)))
#     deserialize(rec,type)
# end

deserialize(ledger::Ledger,type::Type{PeaceFounderConfig}) = readbootrecord(ledger,type)

# function deserialize(ledger::Ledger,type::Type{PeaceFounderConfig})
#     rec = readbootrecord(ledger)
#     #rec = getrecord(ledger,"PeaceFounder.toml")
#     deserialize(rec,type)
# end

function deserialize(chain::BraidChain,::Type{PeaceFounderConfig})
    sc = deserialize(chain.ledger,Certificate{PeaceFounderConfig})
    #@show sc #Is there a bug?
    #@show chain.deme.notary
    intent = Intent(sc,chain.deme.notary)
    @assert intent.reference==chain.deme.spec.maintainer
    return intent.document
end

#serialize(deme::BraidChain,config::PeaceFounderConfig) = serialize(deme.ledger,config)

# function certify(chain::BraidChain,signer::Signer)
#     @assert chain.deme.spec.maintainer==signer.id "You are not eligible to certify PeaceFounder.toml for this deme"
#     #sc = deserialize(chain.ledger,PeaceFounderConfig)
#     sc = readbootrecord(chain.ledger,PeaceFounderConfig)
#     cert = Certificate(sc,signer)
#     #serialize(chain.ledger,cert)
#     writebootrecord(chain.ledger,cert)
# end

ext: jl

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
using Test
using SafeTestsets
using Plots

unicodeplots()

include(joinpath(@__DIR__, "generate_example_tests.jl"))

include(joinpath(@__DIR__, "download_dumps.jl"))

# Note: comment outer @testset to stop after first @safetestset failure
@time @testset verbose = true "Krotov.jl Package" begin

    print("\n* Example 1 (examples/simple_state_to_state.jl):")
    @time @safetestset "Example 1" begin
        include(joinpath("examples", "simple_state_to_state.jl"))
    end

    print("\n* Example 2 (examples/rho_3states.jl):")
    @time @safetestset "Example 2" begin
        include(joinpath("examples", "rho_3states.jl"))
    end

    print("\n* Example 3 (examples/state_to_state_parametrizations.jl):")
    @time @safetestset "Example 3" begin
        include(joinpath("examples", "state_to_state_parametrizations.jl"))
    end

    print("\n* Example 4 (examples/perfect_entanglers.jl):")
    @time @safetestset "Example 4" begin
        include(joinpath("examples", "perfect_entanglers.jl"))
    end

    println("")

end


text:
# Short script to test that the 2 x 3 = 6 different solvers (compared
# in Fig. 1 of the manuscript) work on a medium-size simulated data
# set.
#
# Note that this code was tested in Julia 0.6.2 and may not work in
# later versions of Julia (particularly versions Julia 1.0 and later).
using Distributions
using Mosek
using JuMP
include("../code/datasim.jl");
include("../code/likelihood.jl");
include("../code/primaldual.jl");
include("../code/mixSQP.jl");

# Initialize the pseudorandom number generator.
srand(1);

# Generate a data set.
@printf "Creating data set.\n"
n = round(Int,5e4);
z = normtmixdatasim(n);

# Compute the likelihood matrix.
@printf "Computing likelihood matrix.\n"
sd = autoselectmixsd(z,nv = 20);
L  = normlikmatrix(z,sd = sd);

# Fit the mixture model using the different formulations
# (simplex-constrained, non-negatively-constrained and dual) and
# different solvers (IP and SQP).
@printf "Fitting model...\n"
@printf "simplex, IP:  "
@time x_simplexip  = simplexIP(L);
@printf "non-neg, IP:  "
@time x_nonnegip   = nonnegIP(L);
@printf "dual,    IP:  "
@time x_dualip     = dualIP(L);
@printf "simplex, SQP: "
@time x_simplexsqp = simplexSQP(L);
@printf "non-neg, SQP: \n"
@time x_nonnegsqp = sqp_box(L);
@printf "dual,    SQP: \n"
@time x_dualsqp = sqp_dual(L);

# Check quality of the solutions.
f_simplexip  = mixobjective(L,x_simplexip);
f_nonnegip   = mixobjective(L,x_simplexip);
f_dualip     = mixobjective(L,x_dualip);
f_simplexsqp = mixobjective(L,x_simplexsqp);
f_nonnegsqp  = mixobjective(L,x_nonnegsqp);
f_dualsqp    = mixobjective(L,x_dualsqp);
f_best = minimum([f_simplexip f_nonnegip f_dualip f_simplexsqp f_dualsqp]);
@printf "Difference between objective at given solution, and best solution:\n"
@printf "simplex, IP:  %0.2e\n" f_simplexip - f_best
@printf "non-neg, IP:  %0.2e\n" f_nonnegip - f_best
@printf "dual,    IP:  %0.2e\n" f_dualip - f_best
@printf "simplex, SQP: %0.2e\n" f_simplexsqp - f_best
@printf "dual,    SQP: %0.2e\n" f_dualsqp - f_best

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
# faster than builtin, same behavior

@inline fast_minmax{T<:AbstractFloat}(a::T, b::T) = ifelse(a<b,(a,b),(b,a))
@inline fast_maxmin{T<:AbstractFloat}(a::T, b::T) = ifelse(a>b,(a,b),(b,a))

@inline minmax{T<:AbstractFloat}(a::T, b::T) = ifelse(a<b, (a,b), ifelse(isnan(b),(a,b),(b,a)))
@inline maxmin{T<:AbstractFloat}(a::T, b::T) = ifelse(a>b, (a,b), ifelse(isnan(b),(a,b),(b,a)))

function fast_minmax{T<:AbstractFloat}(a::T, b::T, c::T)
    b,c = fast_minmax(b,c)
    a,c = fast_minmax(a,c)
    a,b = fast_minmax(a,b)
    a,b,c
end
function fast_maxmin{T<:AbstractFloat}(a::T, b::T, c::T)
    b,c = fast_maxmin(b,c)
    a,c = fast_maxmin(a,c)
    a,b = fast_maxmin(a,b)
    a,b,c
end

function minmax{T<:AbstractFloat}(a::T, b::T, c::T)
    if isnan(/4a+b/4-c/4)
        b,c = minmax(b,c)
        a,c = minmax(a,c)
        a,b = minmax(a,b)
    else
        b,c = fast_minmax(b,c)
        a,c = fast_minmax(a,c)
        a,b = fast_minmax(a,b)
    end
    a,b,c
end

function maxmin{T<:AbstractFloat}(a::T, b::T, c::T)
    if isnan(a/4+b/4-c/4)
        b,c = maxmin(b,c)
        a,c = maxmin(a,c)
        a,b = maxmin(a,b)
    else
        b,c = fast_maxmin(b,c)
        a,c = fast_maxmin(a,c)
        a,b = fast_maxmin(a,b)
    end
    a,b,c
end

#=
function minmax{T<:AbstractFloat}(a::T, b::T, c::T)
    b,c = minmax(b,c)
    a,c = minmax(a,c)
    a,b = minmax(a,b)
    a,b,c
end
function maxmin{T<:AbstractFloat}(a::T, b::T, c::T)
    b,c = maxmin(b,c)
    a,c = maxmin(a,c)
    a,b = maxmin(a,b)
    a,b,c
end
=#

#=
function maxmin{T<:AbstractFloat}(a::T, b::T, c::T)
    b,c = ifelse(b<c,(c,b),(b,c))
    a,c = ifelse(a<c,(c,a),(a,c))
    a,b = ifelse(a<b,(b,a),(a,b))
    a,b,c
end
=#

function minmax{T<:AbstractFloat}(a::T,b::T,c::T,d::T)
    a,b = ifelse(a>b,(b,a),(a,b)); c,d = ifelse(c>d,(d,c),(c,d))
    a,c = ifelse(a>c,(c,a),(a,c)); b,d = ifelse(b>d,(d,b),(b,d))
    b,c = ifelse(b>c,(c,b),(b,c))
    a,b,c,d
end

function maxmin{T<:AbstractFloat}(a::T,b::T,c::T,d::T)
    a,b = ifelse(a<b,(b,a),(a,b)); c,d = ifelse(c<d,(d,c),(c,d))
    a,c = ifelse(a<c,(c,a),(a,c)); b,d = ifelse(b<d,(d,b),(b,d))
    b,c = ifelse(b<c,(c,b),(b,c))
    a,b,c,d
end

function minmax{T<:AbstractFloat}(a::T,b::T,c::T,d::T,e::T)
    a,b = ifelse(a>b,(b,a),(a,b)); d,e = ifelse(d>e,(e,d),(d,e))
    c,e = ifelse(c>e,(e,c),(c,e))
    c,d = ifelse(c>d,(d,c),(c,d)); b,e = ifelse(b>e,(e,b),(b,e))
    a,d = ifelse(a>d,(d,a),(a,d))
    a,c = ifelse(c>a,(c,a),(a,c)); b,d = ifelse(b>d,(d,b),(b,d))
    b,c = ifelse(b>c,(c,b),(b,c))
    a,b,c,d,e
end

function maxmin{T<:AbstractFloat}(a::T,b::T,c::T,d::T,e::T)
    a,b = ifelse(a<b,(b,a),(a,b)); d,e = ifelse(d<e,(e,d),(d,e))
    c,e = ifelse(c<e,(e,c),(c,e))
    c,d = ifelse(c<d,(d,c),(c,d)); b,e = ifelse(b<e,(e,b),(b,e))
    a,d = ifelse(a<d,(d,a),(a,d))
    a,c = ifelse(c<a,(c,a),(a,c)); b,d = ifelse(b<d,(d,b),(b,d))
    b,c = ifelse(b<c,(c,b),(b,c))
    a,b,c,d,e
end

#=
http://pages.ripco.net/~jgamble/nw.html

6
[[1,2],[4,5]]
[[0,2],[3,5]]
[[0,1],[3,4],[2,5]]
[[0,3],[1,4]]
[[2,4],[1,3]]
[[2,3]]

7
[[1,2],[3,4],[5,6]]
[[0,2],[3,5],[4,6]]
[[0,1],[4,5],[2,6]]
[[0,4],[1,5]]
[[0,3],[2,5]]
[[1,3],[2,4]]
[[2,3]]

8
[[0,1],[2,3],[4,5],[6,7]]
[[0,2],[1,3],[4,6],[5,7]]
[[1,2],[5,6],[0,4],[3,7]]
[[1,5],[2,6]]
[[1,4],[3,6]]
[[2,4],[3,5]]
[[3,4]]
=#


text:
function verifyBaseNodeAttributes!( mpSim::MPsim )

    missingAttributes = Dict{String, Vector{String}}()
    missingValues = Vector{NTuple{3, String}}()

    # Verify the consistency of the attributes in each base node.
    for nodeName in keys( mpSim.baseNodeList )
        node = mpSim.baseNodeList[nodeName]

        for attribute in keys( node.requirements )
            value = node.requirements[attribute]
            if haskey( missingAttributes, attribute )
                push!( missingAttributes[attribute], nodeName )
            elseif !haskey( mpSim.attributeList, attribute )
                missingAttributes[attribute] = [nodeName]
            elseif !isAttributeValuePossible( mpSim.attributeList[attribute],
                value )
                push!( missingValues, (nodeName, attribute, value) )
            end  # if haskey( missingAttributes, attribute )
        end  # for attribute in node.requirements
    end  # for nodeName in keys( mpSim.baseNodeList )

    # Report missing attributes in base node requirements.
    warnString = map( collect( keys( missingAttributes ) ) ) do attribute
        return string( "Unknown attribute '", attribute, "' in base node(s) '",
            join( missingAttributes[attribute], "', '", "', and '" ), "'" )
    end  # map( ... ) do attribute

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

    # Report missing values in base node requirements.
    warnString = map( missingValues ) do missingVal
        return string( "Unknown value '", missingVal[3], "' for attribute '",
            missingVal[2], "' in base node '", missingVal[1], "'" )
    end  # map( missingValues ) do missingVal

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyBaseNodeAttributes!( mpSim )

function verifyBaseNodeAttrition!( mpSim::MPsim )

    # Verify consistency of attrition schemes in base nodes.
    missingAttrition = Dict{String, Vector{String}}()

    for nodeName in keys( mpSim.baseNodeList )
        node = mpSim.baseNodeList[nodeName]
        attrition = node.attrition

        if !haskey( mpSim.attritionSchemes, attrition )
            if haskey( missingAttrition, attrition )
                push!( missingAttrition[attrition], node.name )
            else
                missingAttrition[attrition] = [node.name]
            end  # if haskey( missingAttrition, attrition )
        end  # if !haskey( mpSim.attritionSchemes, attrition )
    end  # for node in mpSim.baseNodeList

    # Report missing attrition schemes in base nodes.
    warnString = map( collect( keys( missingAttrition ) ) ) do attrition
        return string( "Unknown attrition scheme '", attrition,
            "' in base nodes '", join( missingAttrition[attrition], "', '",
            "', and '" ), "'" )
    end  # map( ... ) do attrition

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyBaseNodeAttrition!( mpSim )

function verifyTransitionBaseNodes!( mpSim::MPsim )

    # Verify consistency of base nodes in each transition.
    missingSourceNodes = Dict{String, Vector{String}}()
    missingTargetNodes = Dict{String, Vector{String}}()

    for node in keys( mpSim.transitionsBySource )
        if !haskey( mpSim.baseNodeList, node )
            missingSourceNodes[node] = map( transition -> transition.name,
                mpSim.transitionsBySource[node] )
        end  # if !haskey( mpSim.baseNodeList, node )
    end  # for node in keys( mpSim.transitionsBySource )

    for node in keys( mpSim.transitionsByTarget )
        if ( node != "OUT" ) && !haskey( mpSim.baseNodeList, node )
            missingTargetNodes[node] = map( transition -> transition.name,
                mpSim.transitionsByTarget[node] )
        end  # if !haskey( mpSim.baseNodeList, node )
    end  # for node in keys( mpSim.transitionsByTarget )

    # Report missing source/target nodes in compound node composition.
    warnString = map( collect( keys( missingSourceNodes ) ) ) do node
        return string( "Unknown source node '", node, "' in transition(s) '",
            join( missingSourceNodes[node], "', '", "', and '" ), "'" )
    end  # warnString = map( ... ) do node

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

    warnString = map( collect( keys( missingTargetNodes ) ) ) do node
        return string( "Unknown target node '", node, "' in transition(s) '",
            join( missingTargetNodes[node], "', '", "', and '" ), "'" )
    end  # warnString = map( ... ) do node

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyTransitionBaseNodes!( mpSim )

function verifyTransitionConditions!( mpSim::MPsim )

    # Verify consistency of attributes in each transition.
    missingAttributes = Dict{String, Vector{String}}()
    missingValues = Vector{NTuple{3, String}}()

    for transitionName in keys( mpSim.transitionsByName ),
        transition in mpSim.transitionsByName[transitionName],
        condition in transition.extraConditions

        attribute = condition.attribute
        isTime = lowercase( attribute ) ∈ timeAttributes

        if !isTime && !haskey( mpSim.attributeList, attribute )
            if haskey( missingAttributes, attribute )
                push!( missingAttributes[attribute], transitionName )
            else
                missingAttributes[attribute] = [transitionName]
            end  # if haskey( missingAttributes, attribute )
        elseif !isTime
            attribute = mpSim.attributeList[attribute]

            if ( condition.operator == Base.:(==) ) &&
                !isAttributeValuePossible( attribute, condition.value )
                !haskey( mpSim.attributeList, condition.value )
                push!( missingValues,
                    (transitionName, attribute.name, condition.value) )
            elseif condition.operator == Base.:∈
                for value in filter(
                    val -> !isAttributeValuePossible( attribute, val ),
                    condition.value )
                    push!( missingValues, (transitionName, attribute.name,
                        value) )
                end  # for value in filter( ... )
            end  # if ( condition.operator == Base.:(==) ) && ...
        end  # if !isTime && ...
    end  # for transitionName in mpSim.transitionsByName, ...

    # Report missing attributes in base node requirements.
    warnString = map( collect( keys( missingAttributes ) ) ) do attribute
        return string( "Unknown attribute '", attribute, "' in base node(s) '",
            join( missingAttributes[attribute], "', '", "', and '" ), "'" )
    end  # map( ... ) do attribute

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

    # Report missing values in base node requirements.
    warnString = map( missingValues ) do missingVal
        return string( "Unknown value '", missingVal[3], "' for attribute '",
            missingVal[2], "' in base node '", missingVal[1], "'" )
    end  # map( missingValues ) do missingVal

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyTransitionConditions!( mpSim )

function verifyTransitionChanges!( mpSim::MPsim )

    missingAttributes = Dict{String, Vector{String}}()
    missingAttrVals = Dict{String, Vector{String}}()

    for name in keys( mpSim.transitionsByName ),
        transition in mpSim.transitionsByName[name],
        attrName in keys( transition.extraChanges )
        # Check if the attribute exists.
        if !haskey( mpSim.attributeList, attrName )
            if haskey( missingAttributes, attrName )
                push!( missingAttributes[attrName], name )
            else
                missingAttributes[attrName] = [name]
            end  # if haskey( missingAttributes, attrName )

            continue
        end  # if !haskey( mpSim.attributeList, attrName )

        # check if the new value is a possible value for the attribute.
        attribute = mpSim.attributeList[attrName]
        val = transition.extraChanges[attrName]

        if val ∉ attribute.possibleValues
            if haskey( missingAttrVals, attrName )
                push!( missingAttrVals[attrName], val )
            else
                missingAttrVals[attrName] = [val]
            end  # if haskey( missingAttrVals, attrName )
        end  # if transition.extraChanges[attrName] ∉ attribute.possibleValues
    end  # for name in keys( mpSim.transitionsByName ), ...

    # Report missing attributes in extra changes of transition.
    warnString = map( collect( keys( missingAttributes ) ) ) do name
        return string( "Unknown attribute '", name,
        "' to be changed in transitions '",
        join( unique( missingAttributes[name] ), "', '", "', and '" ), "'." )
    end  # map( collect( keys( missingAttributes ) ) ) do name

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

    # Report missing attribute values in extra changes of transition.
    warnString = map( collect( keys( missingAttrVals ) ) ) do name
        return string( "Unknown value(s) '",
            join( missingAttrVals[name], "', '", "', and '" ),
            "' in attribute '", name, "'." )
    end  # map( collect( keys( missingAttributes ) ) ) do name

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyTransitionChanges!( mpSim )

function verifyCompoundNodeComponents!( mpSim::MPsim )

    # Verify consistency of base nodes in each compound node.
    missingNodes = Dict{String, Vector{String}}()

    for nodeName in keys( mpSim.compoundNodeList )
        node = mpSim.compoundNodeList[nodeName]

        for baseNode in node.baseNodeList
            if haskey( missingNodes, baseNode )
                push!( missingNodes[baseNode], nodeName )
            elseif !haskey( mpSim.baseNodeList, baseNode )
                missingNodes[baseNode] = [nodeName]
            end  # if haskey( missingNodes, baseNode )
        end  # for baseNode in node.baseNodeList
    end  # for nodeName in keys( mpSim.compoundNodeList )

    # Report missing base nodes in compound node composition.
    warnString = map( collect( keys( missingNodes ) ) ) do node
        return string( "Unknown base node '", node, "' in compound node(s) '",
            join( missingNodes[node], "', '", "', and '" ), "'" )
    end  # warnString = map( ... ) do node

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyCompoundNodeComponents!( mpSim )

function verifyRecruitmentTargets!( mpSim::MPsim )

    # Verify target nodes of recruitment schemes.
    missingNodes = Dict{String, Vector{String}}()

    for nodeName in keys( mpSim.recruitmentByTarget )
        if !haskey( mpSim.baseNodeList, nodeName )
            missingNodes[nodeName] = map( recruitment -> recruitment.name,
                mpSim.recruitmentByTarget[nodeName] )
        end  # if !haskey( mpSim.baseNodeList, nodeName )
    end  # for nodeName in keys( mpSim.recruitmentByTarget )

    # Report missing target nodes of recruitment schemes.
    warnString = map( collect( keys( missingNodes ) ) ) do node
        return string( "Unknown target node '", node,
            "' in recruitment scheme(s) '",
            join( missingNodes[node], "', '", "', and '" ), "'" )
    end  # map( ... ) do node

    if !isempty( warnString )
        @warn join( warnString, "\n" )
        mpSim.isConsistent = false
    end  # if !isempty( warnString )

end  # verifyRecruitmentTargets!( mpSim )

function setSimAttributes!( mpSim::MPsim, attributes::Vector{Attribute},
    wipeConfig::Bool )

    attributeNames = getfield.( attributes, :name )

    if length( attributeNames ) != length( unique( attributeNames ) )
        @warn "Multiple attributes with same name in attribute list, not making any changes."
        return false
    end  # if length( attributeNames ) != length( unique( attributeNames ) )

    if wipeConfig
        clearSimulationAttributes!( mpSim )
    end  # wipeConfig

    for attribute in attributes
        mpSim.attributeList[attribute.name] = attribute
    end  # for attribute in attributes

    mpSim.isStale = true
    return true

end  # setSimAttributes!( mpSim, attributes, wipeConfig )

function setSimBaseNodes!( mpSim::MPsim, nodes::Vector{BaseNode},
    wipeConfig::Bool )

    nodeNames = getfield.( nodes, :name )

    if length( nodeNames ) != length( unique( nodeNames ) )
        @warn "Multiple nodes with same name in base node list, not making any changes."
        return false
    end  # if length( nodeNames ) != length( unique( nodeNames ) )

    if wipeConfig
        clearSimulationBaseNodes!( mpSim )
    end  # if wipeConfig

    for node in nodes
        mpSim.baseNodeList[node.name] = node
    end  # for node in nodes

    mpSim.isStale = true
    return true

end  # setSimBaseNodes!( mpSim, nodes, wipeConfig )

function setSimCompoundNodes!( mpSim::MPsim, nodes::Vector{CompoundNode},
    wipeConfig::Bool )

    nodeNames = getfield.( nodes, :name )

    if length( nodeNames ) != length( unique( nodeNames ) )
        @warn "Multiple nodes with same name in base node list, not making any changes."
        return false
    end  # if length( nodeNames ) != length( unique( nodeNames ) )

    if wipeConfig
        clearSimulationCompoundNodes!( mpSim )
    end  # if wipeConfig

    for node in nodes
        mpSim.compoundNodeList[node.name] = node
    end  # for node in nodes

    mpSim.isStale = true
    return true

end  # setSimCompoundNodes!( mpSim, nodes, wipeConfig )

function setSimRecruitment!( mpSim::MPsim, recruitmentList::Vector{Recruitment},
    wipeConfig::Bool )

    tmpRecruitment = filter( recruitment ->
        recruitment.targetNode ∉ ["", "dummy"], recruitmentList )

    if isempty( tmpRecruitment )
        return false
    end  # if isempty( tmpRecruitment )

    if wipeConfig
        clearSimulationRecruitment!( mpSim )
    end  # if wipeConfig

    # Add the recruitment schemes to the appropriate lists.
    for recruitment in tmpRecruitment
        if haskey( mpSim.recruitmentByName, recruitment.name )
            push!( mpSim.recruitmentByName[recruitment.name], recruitment )
        else
            mpSim.recruitmentByName[recruitment.name] = [recruitment]
        end  # if haskey( mpSim.recruitmentByName, recruitment.name )

        if haskey( mpSim.recruitmentByTarget, recruitment.targetNode )
            push!( mpSim.recruitmentByTarget[recruitment.targetNode],
                recruitment )
        else
            mpSim.recruitmentByTarget[recruitment.targetNode] =
                [recruitment]
        end  # if haskey( mpSim.recruitmentByTarget, recruitment.targetNode )
    end  # for recruitment in tmpRecruitment

    mpSim.isStale = true
    return true

end  # setSimRecruitment!( mpSim, recruitmentList, wipeConfig )

function setSimTransitions!( mpSim::MPsim, transitions::Vector{Transition},
    wipeConfig::Bool )

    tmpTransitions = filter( transitions ) do transition
        return ( lowercase( transition.sourceNode ) ∉ ["", "dummy"] ) &&
            ( transition.isOutTransition ? true :
            lowercase( transition.targetNode ) ∉ ["", "dummy"] )
    end  # filter( transitions ) do transition

    if isempty( tmpTransitions )
        return false
    end  # if isempty( tmpTransitions )

    if wipeConfig
        clearSimulationTransitions!( mpSim )
    end  # if wipeConfig

    # Add the transitions to the appropriate lists.
    for transition in tmpTransitions
        if haskey( mpSim.transitionsByName, transition.name )
            push!( mpSim.transitionsByName[transition.name], transition )
        else
            mpSim.transitionsByName[transition.name] = [transition]
        end  # if haskey( mpSim.transitionsByName, transitions.name )

        if haskey( mpSim.transitionsBySource, transition.sourceNode )
            push!( mpSim.transitionsBySource[transition.sourceNode],
                transition )
        else
            mpSim.transitionsBySource[transition.sourceNode] = [transition]
        end  # if haskey( mpSim.transitionsBySource, transitions.sourceNode )

        if transition.isOutTransition
            push!( mpSim.transitionsByTarget["OUT"], transition )
        elseif haskey( mpSim.transitionsByTarget, transition.targetNode )
            push!( mpSim.transitionsByTarget[transition.targetNode],
                transition )
        else
            mpSim.transitionsByTarget[transition.targetNode] =
                [transition]
        end  # if transition.isOutTransition
    end  # for transition in tmpTransitions

    mpSim.isStale = true
    return true

end  # setSimTransitions!( mpSim, transitions, wipeConfig )

function setSimTransitionTypeOrder!( mpSim::MPsim,
    transOrder::Dict{String, Int}, wipeConfig::Bool )

    if isempty( transOrder )
        return false
    end  # if isempty( transOrder )

    if wipeConfig
        clearSimulationTransitionTypeOrder!( mpSim )
    end   # if wipeConfig

    for name in keys( transOrder )
        mpSim.transitionTypeOrder[name] = transOrder[name]
    end  # for name in keys( transOrder )

    return true

end  # setSimTransitionTypeOrder!( mpSim, transOrder, wipeConfig )

function setSimBaseNodeOrder!( mpSim::MPsim,
    nodeOrder::Dict{String, Int}, wipeConfig::Bool )

    if isempty( nodeOrder )
        return false
    end  # if isempty( nodeOrder )

    if wipeConfig
        clearSimulationBaseNodeOrder!( mpSim )
    end   # if wipeConfig

    for name in keys( nodeOrder )
        mpSim.baseNodeOrder[name] = nodeOrder[name]
    end  # for name in keys( nodeOrder )

    return true

end  # setSimBaseNodeOrder!( mpSim, nodeOrder, wipeConfig )

function setSimAttrition!( mpSim::MPsim, attritionList::Vector{Attrition},
    wipeConfig::Bool )

    attritionNames = getfield.( attritionList, :name )
    nDefaults = count( map( attrition -> lowercase( attrition ) ∈
        ["", "default"], attritionNames ) )

    if ( length( attritionNames ) != length( unique( attritionNames ) ) ) ||
        ( nDefaults > 1 )
        @warn "Multiple attrition schemes with same name in attribute list, not making any changes."
        return false
    end  # if length( attritionNames ) != length( unique( attritionNames ) )

    if wipeConfig
        clearSimulationAttrition!( mpSim )
    end  # if wipeConfig

    for attrition in attritionList
        attritionName = lowercase( attrition.name ) ∈ ["", "default"] ?
            "default" : attrition.name
        mpSim.attritionSchemes[attritionName] = deepcopy( attrition )
    end  # for attrition in attritionList

    mpSim.isStale = true
    return true

end  # setSimAttrition!( mpSim, attritionList, wipeConfig )

function validateDatabaseAge!( mpSim::MPsim )

    if haskey( SQLite.tables( mpSim.simDB ), :name ) &&
        ( mpSim.persDBname ∈ SQLite.tables( mpSim.simDB )[:name] )
        mpSim.isOldDB = "startState" ∈ SQLite.columns( mpSim.simDB,
            mpSim.transDBname )[:name]
    else
        mpSim.isOldDB = false
    end  # haskey( SQLite.tables( mpSim.simDB ), :name ) && ...

    if mpSim.isOldDB
        @warn "Old style simulation results database. This style will be deprecated in a future version."
    end  # if mpSim.isOldDB

    mpSim.sNode = mpSim.isOldDB ? "startState" : "sourceNode"
    mpSim.tNode = mpSim.isOldDB ? "endState" : "targetNode"
    mpSim.valName = mpSim.isOldDB ? "strValue" : "value"

end  # validateDatabaseAge!( mpSim )


text:
@testset "Grid Test" begin
    @test ALifeBenchmark._get_grid_coordinates((1, 1)) == (1, 1)
    @test ALifeBenchmark._get_grid_coordinates((0.5, 0.5)) == (1, 1)
    @test ALifeBenchmark._get_grid_coordinates((6.5, 3)) == (7, 3)
    @test ALifeBenchmark._get_grid_coordinates((2.5, 3.4)) == (3, 4)

    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    @test ALifeBenchmark._get_cell_in_direction(model, (0.5, 0.5), 0) == (2,  1)
    @test ALifeBenchmark._get_cell_in_direction(model, (0.5, 0.5), 90) == (1, 2)
    @test ALifeBenchmark._get_cell_in_direction(model, (0.5, 0.5), 180) == (2, 1)
    @test ALifeBenchmark._get_cell_in_direction(model, (0.5, 0.5), 270) == (1, 2)

    @test ALifeBenchmark._get_cell_in_direction(model, (1.5, 0.5), 0) == (1,  1)
    @test ALifeBenchmark._get_cell_in_direction(model, (1.5, 0.5), 90) == (2, 2)
    @test ALifeBenchmark._get_cell_in_direction(model, (1.5, 0.5), 180) == (1, 1)
    @test ALifeBenchmark._get_cell_in_direction(model, (1.5, 0.5), 270) == (2, 2)

    organism_11.direction = 0
    @test ALifeBenchmark._get_organism_in_front(model, organism_11) == organism_21
    @test ALifeBenchmark._get_organism_behind(model, organism_11) == organism_21

    organism_22.direction = 90
    @test ALifeBenchmark._get_organism_in_front(model, organism_22) == organism_21
    @test ALifeBenchmark._get_organism_behind(model, organism_22) == organism_21
end

@testset "Cross-Over Test" begin
    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    organism_11.direction = 0

    ALifeBenchmark.perform!(ALifeBenchmark.CrossOver(), model, organism_11, 1)
    @test model.grid[2, 1] != organism_21
    @test model.grid[2, 1].genotype == "00"

    organism_12.genotype = "01"
    organism_22.genotype = "11"
    organism_12.direction = 0
    ALifeBenchmark.perform!(ALifeBenchmark.CrossOver(), model, organism_12, 1)
    @test model.grid[2, 2] != organism_22
    @test model.grid[2, 2].genotype == "011"

    organism_11.genotype = "011"
    organism_12.genotype = "0011"
    organism_11.direction = 270
    ALifeBenchmark.perform!(ALifeBenchmark.CrossOver(), model, organism_11, 1)
    @test model.grid[1, 2] != organism_12
    @test model.grid[1, 2].genotype in ["00011", "011", "01011", "011"]

    organism_11.genotype = "011"
    model.grid[1, 2].genotype = "0011"
    organism_11.direction = 270
    ALifeBenchmark.perform!(ALifeBenchmark.CrossOver(), model, organism_11, 1)
    @test model.grid[1, 2].genotype in ["00011", "011", "01011", "011"]

    organism_11.genotype = "011"
    model.grid[1, 2].genotype = "0011"
    organism_11.direction = 270
    ALifeBenchmark.perform!(ALifeBenchmark.CrossOver(), model, organism_11, 1)
    @test model.grid[1, 2].genotype in ["00011", "011", "01011", "011"]
end

@testset "Fight Test" begin
    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    organism_11.direction = 0

    ALifeBenchmark.perform!(ALifeBenchmark.Fight(), model, organism_11, 1)
    @test model.grid[2, 1] === nothing
end

@testset "Turning Test" begin
    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    organism_11.direction = 0

    ALifeBenchmark.perform!(ALifeBenchmark.TurnClockwise(), model, organism_11, 20/360)
    @test organism_11.direction == 3
    ALifeBenchmark.perform!(ALifeBenchmark.TurnClockwise(), model, organism_11, 350/360)
    @test organism_11.direction == 61
    ALifeBenchmark.perform!(ALifeBenchmark.TurnAntiClockwise(), model, organism_11, 40/360)
    @test organism_11.direction == 54
    ALifeBenchmark.perform!(ALifeBenchmark.TurnAntiClockwise(), model, organism_11, 100/360)
    @test organism_11.direction == 37
end

@testset "Moving Test" begin
    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    organism_11.direction = 0

    @test organism_11.coordinates == (0.5, 0.5)

    ALifeBenchmark.perform!(ALifeBenchmark.MoveForward(), model, organism_11, 0.1)
    @test organism_11.coordinates == (0.6, 0.5)

    ALifeBenchmark.perform!(ALifeBenchmark.MoveForward(), model, organism_11, 0.5)
    @test organism_11.coordinates == (0.6, 0.5)

    ALifeBenchmark.kill!(model, organism_21)
    ALifeBenchmark.perform!(ALifeBenchmark.MoveForward(), model, organism_11, 0.5)
    @test organism_11.coordinates == (1.1, 0.5)
    @test model.grid[1, 1] === nothing
    @test model.grid[2, 1] === organism_11
end

@testset "Complete Action Test" begin
    Random.seed!(100)

    # build network
    node1 = ALifeBenchmark.Node("1", [], [], [], [])
    node2 = ALifeBenchmark.Node("2", [], [], [], [])
    node3 = ALifeBenchmark.Node("3", [], [], [], [])
    node4 = ALifeBenchmark.Node("4", [], [], [], [])
    node5 = ALifeBenchmark.Node("5", [], [], [], [])
    node6 = ALifeBenchmark.Node("1112", [], [], [], [])  # move forward
    node7 = ALifeBenchmark.Node("11001", [], [], [], []) # turn clockwise

    node1.type = ALifeBenchmark.InputNode()
    node2.type = ALifeBenchmark.InputNode()
    node6.type = ALifeBenchmark.OutputNode()
    node7.type = ALifeBenchmark.OutputNode()

    push!(node1.out_excitatory, node3)
    push!(node1.out_inhibitory, node4)

    push!(node2.out_excitatory, node3)
    push!(node2.out_excitatory, node4)

    push!(node3.in_excitatory, node1)
    push!(node3.in_excitatory, node2)
    push!(node3.out_excitatory, node5)

    push!(node4.in_inhibitory, node1)
    push!(node4.in_excitatory, node2)
    push!(node4.out_excitatory, node7)

    push!(node5.in_excitatory, node3)
    push!(node5.out_excitatory, node6)
    push!(node5.out_excitatory, node7)

    push!(node6.in_excitatory, node5)
    push!(node6.out_excitatory, node7)

    push!(node7.in_excitatory, node4)
    push!(node7.in_excitatory, node5)
    push!(node7.in_excitatory, node6)

    network = ALifeBenchmark.Network()
    ALifeBenchmark.update_inputs_outputs!(network, [node1])
    ALifeBenchmark.fire!(network, [0, 0])
    ALifeBenchmark.fire!(network, [0, 0])
    ALifeBenchmark.fire!(network, [0, 0])
    ALifeBenchmark.fire!(network, [0, 0])
    ALifeBenchmark.fire!(network, [0, 0])

    # build Geb model

    model = ALifeBenchmark.GebModel(size=2)

    organism_11 = model.grid[1, 1]
    organism_12 = model.grid[1, 2]
    organism_21 = model.grid[2, 1]
    organism_22 = model.grid[2, 2]

    organism_11.network = network
    organism_11.direction = 0

    @assert organism_11.coordinates == (0.5, 0.5)

    ALifeBenchmark.perform!(model, organism_11)

    @test organism_11.coordinates == (0.6233837738869451, 0.5)
    @test organism_11.direction == 15
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# Tests data_import.jl

using Test
using GeophysicalModelGenerator

# should throw an error with a 2D dataset
Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,-50km);
Data1           =   Depth*2;                # some data
Vx1,Vy1,Vz1     =   Data1*3,Data1*4,Data1*5
Data_set2D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon, Velocity=(Vx1,Vy1,Vz1)))
@test_throws ErrorException CrossSection(Data_set2D, Depth_level=-10)

# Create 3D volume with some fake data
Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
Data            =   Depth*2;                # some data
Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)))

# Create 3D volume with some fake data
Lon,Lat,Depth           =   LonLatDepthGrid(10:20,30:40,(0:-25:-300)km);
Data                    =   Depth*2;                # some data
Vx,Vy,Vz                =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
Data_set3D_reverse      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)))

# Create cross-sections in various directions (no interpolation which is default)
test_cross      =   CrossSection(Data_set3D, Depth_level=-100km)
@test test_cross.fields[1][41]==-200km
@test test_cross.fields[2][31]==18
@test test_cross.fields[3][1][30]==-600
@test test_cross.fields[3][2][30]==-800
@test test_cross.fields[3][3][30]==-1000

# throw error if outside bounds
@test_throws ErrorException CrossSection(Data_set3D, Depth_level=100km)

test_cross      =   CrossSection(Data_set3D, Lon_level=15)
@test test_cross.fields[1][41]==-450km
@test test_cross.fields[2][31]==15
@test test_cross.fields[3][1][30]==-1500
@test test_cross.fields[3][2][30]==-2000
@test test_cross.fields[3][3][30]==-2500

test_cross      =   CrossSection(Data_set3D, Lat_level=35)
@test test_cross.fields[1][41]==-450km
@test test_cross.fields[2][31]==18
@test test_cross.fields[3][1][30]==-1500
@test test_cross.fields[3][2][30]==-2000
@test test_cross.fields[3][3][30]==-2500

# Create cross-sections with interpolation in various directions
test_cross      =   CrossSection(Data_set3D, Depth_level=-100km, dims=(50,100), Interpolate=true)
@test size(test_cross.fields[1])    ==  (50,100,1)
@test size(test_cross.fields[3][2]) ==  (50,100,1)

test_cross      =   CrossSection(Data_set3D, Lon_level=15, dims=(50,100), Interpolate=true)
@test size(test_cross.fields[3][2])==(1,50,100)
@test Write_Paraview(test_cross, "profile_test")==nothing

test_cross      =   CrossSection(Data_set3D, Lat_level=35, dims=(50,100), Interpolate=true)
@test size(test_cross.fields[3][2])==(50,1,100)

# Diagonal cross-section
test_cross      =   CrossSection(Data_set3D, Start=(10,30), End=(20,40), dims=(50,100), Interpolate=true)
@test size(test_cross.fields[3][2])==(50,100,1)
@test Write_Paraview(test_cross, "profile_test")==nothing

#test_cross_rev  =   CrossSection(Data_set3D_reverse, Start=(10,30), End=(20,40), dims=(50,100), Interpolate=true)
#@test size(test_cross_rev.fields[3][2])==(50,100,1)
#@test Write_Paraview(test_cross_rev, "profile_test_rev")[1]=="profile_test_rev.vts"

# Extract sub-volume

# with interpolation
Data_sub_Interp = ExtractSubvolume(Data_set3D,Lon_level=(10,15), Lat_level=(30,32), Interpolate=true, dims=(51,21,32))
@test Data_sub_Interp.fields[1][11]==-600km
@test size(Data_sub_Interp.lat)==(51,21,32)

# no interpolation
Data_sub_NoInterp = ExtractSubvolume(Data_set3D,Lon_level=(10,15), Lat_level=(30,32), Interpolate=false, dims=(51,21,32))
@test Data_sub_NoInterp.fields[1][11]==-600km
@test size(Data_sub_NoInterp.lat)==(6,3,13)

# Extract subset of cross-section
test_cross      =   CrossSection(Data_set3D, Lat_level=35, dims=(50,100), Interpolate=true)
Data_sub_cross  =   ExtractSubvolume(test_cross, Depth_level=(-100km,0km), Interpolate=false)
@test Data_sub_cross.fields[1][11]==-200.00000000000003km
@test size(Data_sub_cross.lat)==(50,1,34)

# compute the mean velocity per depth in a 3D dataset and subtract the mean from the given velocities
Data_pert   =   SubtractHorizontalMean(ustrip(Data))    # 3D, no units
@test Data_pert[10] == 0.0

Data_pert   =   SubtractHorizontalMean(Data)            # 3D with units
@test Data_pert[10] == 0.0km

Data_pert   =   SubtractHorizontalMean(Data, Percentage=true)            # 3D with units
@test Data_pert[1000] == 0.0

Data2D = Data[:,1,:];
Data_pert   =   SubtractHorizontalMean(Data2D, Percentage=true)         # 2D version with units [dp the same along a vertical profile]

Data_set2D  =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon,Pertdata=Data_pert ,Velocity=(Vx,Vy,Vz)))
@test Data_set2D.fields[3][10,8,1] == 0

# Create surface ("Moho")
Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,-40km);
Depth           =   Depth + Lon*km;     # some fake topography on Moho
Data_Moho       =   GeoData(Lon,Lat,Depth,(MohoDepth=Depth,LonData=Lon))

# Test intersecting a surface with 2D or 3D data sets
Above       =   AboveSurface(Data_set3D, Data_Moho);            # 3D regular ordering
@test Above[1,1,12]==true
@test Above[1,1,11]==false

Above       =   AboveSurface(Data_set3D_reverse, Data_Moho);    #  3D reverse depth ordering
@test Above[1,1,2]==true
@test Above[1,1,3]==false

Above       =   AboveSurface(Data_sub_cross, Data_Moho);        # 2D cross-section
@test Above[end]==true
@test Above[1]==false

# Test VoteMaps
Data_VoteMap = VoteMap(Data_set3D, "Depthdata<-560",dims=(10,10,10))
@test Data_VoteMap.fields[:VoteMap][101]==0
@test Data_VoteMap.fields[:VoteMap][100]==1

Data_VoteMap = VoteMap(Data_set3D_reverse, "Depthdata<-560",dims=(10,10,10))
@test Data_VoteMap.fields[:VoteMap][101]==0
@test Data_VoteMap.fields[:VoteMap][100]==1

# Combine 2 datasets
Data_VoteMap = VoteMap([Data_set3D_reverse, Data_set3D], ["Depthdata<-560","LonData>19"],dims=(10,10,10))
@test Data_VoteMap.fields[:VoteMap][10,9,1]==2
@test Data_VoteMap.fields[:VoteMap][9 ,9,1]==1
@test Data_VoteMap.fields[:VoteMap][9 ,9,2]==0

# Test rotation routines
X,Y,Z   =   LonLatDepthGrid(10:20,30:40,-50:-10);
Data_C  =   ParaviewData(X,Y,Z,(Depth=Z,))
RotateTranslateScale!(Data_C, Rotate=30);

@test Data_C.x.val[10] ≈ 20.964101615137753
@test Data_C.y.val[10] ≈ 32.66987298107781
@test Data_C.z.val[20] == -50

RotateTranslateScale!(Data_C, Scale=10, Rotate=10, Translate=(1,2,3));
@test Data_C.x.val[10] ≈ 213.78115820908607
@test Data_C.y.val[10] ≈ 339.4092822315127
@test Data_C.z.val[20] == -497.0

# Test


text:
@testset verbose = true "fixedpoint" begin
    @testset "supnorm" begin
        @test supnorm(ones(3), ones(3)) == 0.0
        @test supnorm(ones(3, 3, 3), ones(3, 3, 3)) == 0.0
        @test supnorm([1, 2, 3], [4, 4, 4]) == 3
    end

    @testset "fixedpoint" begin
        init = ones(10)
        cache = (ones(10) .* 2, ones(10) .* 3)
        function iterate!(x1, x0, β)
            return x1 .= β .* x0
        end
        # test a fixed point which converges
        res = fixedpoint((x1, x0) -> iterate!(x1, x0, 0.5), init, cache)
        @test res.converged
        @test supnorm(res.value, zeros(10)) <= 1e-10

        # test that the cache works
        res = fixedpoint((x1, x0) -> iterate!(x1, x0, 0.2), init)
        @test res.converged
        @test supnorm(res.value, zeros(10)) <= 1e-10

        # test a fixed point which doesn't converge
        res = fixedpoint((x1, x0) -> iterate!(x1, x0, 1.2), init, cache; maxiter=50)
        @test !res.converged
        @test res.iter == 50
    end
end


text:
"""
    forward_score(c::CRF, x::Array)

Compute the Normalization / partition function
or the Forward Algorithm score - `Z`
"""
function forward_score(c::CRF, x, init_α)
    forward_var = log_sum_exp(c.W .+ x[1]' .+ init_α)

    for i in 2:length(x)
        forward_var = log_sum_exp((c.W .+ x[i]') .+ forward_var')
    end

    return log_sum_exp(c.W[:, c.n + 2] .+ forward_var')[1]
end

"""
    score_sequence(c::CRF, xs, label_seq)

Calculating the score of the desired `label_seq` against sequence `xs`.
Not exponentiated as required for negative log likelihood,
thereby preventing operation.

`label_seq`<:Array/ CuArray
eltype(label_seq) = Flux.OneHotVector
"""
function score_sequence(c::CRF, x, label_seq)
    score = preds_first(c, label_seq[1]) + onecold(label_seq[1], x[1])

    for i in 2:length(label_seq)
        score += preds_single(c, label_seq[i], label_seq[i-1]) +
                    onecold(label_seq[i], x[i])
    end

    return score + preds_last(c, label_seq[end])
end

# REGULARIZATION TERM
crf_loss(c::CRF, x, label_seq, init_α) = forward_score(c, x, init_α) - score_sequence(c, x, label_seq)


text:
@time using Plots, LaTeXStrings
@time using CSV, DataFrames
@time using Statistics

default(alpha = 0.2)

last_seed = 30

directory = "D:/trash 0.20/"; cd(directory); println(pwd())
scenario_list = ["00", "0V", "M0", "MV"]
# scenario_list = ["00"]
scenario_color = Dict([
    "00" => :black,
    "0V" => colorant"#C00000",
    "M0" => colorant"#0070C0",
    "MV" => colorant"#7030A0"]
)
scenario_shape = Dict([
    "00" => :circle,
    "0V" => :utriangle,
    "M0" => :square,
    "MV" => :star]
)

l = @layout [[I; R; Rt] RV]

plot_tevl_I = plot(ylabel = L"I(T)")
plot_tevl_R = plot(ylabel = L"R(T)")
plot_tevl_Rt = plot(xlabel = L"T", ylabel = L"R_T")
plot_smry_RV = plot(xlabel = L"\#R", ylabel = L"\#V", size = (500,500))

plot_trng_I = plot(ylabel = L"I(T)")
plot_trng_R = plot(ylabel = L"R(T)")
plot_trng_Rt = plot(xlabel = L"T", ylabel = L"R_T")

tevl_I = zeros(200, last_seed)
tevl_R = zeros(200, last_seed)
tevl_Rt = zeros(200, last_seed)
smry_0V = zeros(last_seed, 2)
smry_MV = zeros(last_seed, 2)

for scenario_code ∈ scenario_list, seed_number ∈ 1:last_seed

temp_tevl = CSV.read("./$scenario_code//$(lpad(seed_number, 4, '0')) tevl.csv", DataFrame)
if length(temp_tevl.I) != 200 continue end
tevl_I[:, seed_number] = temp_tevl.I
tevl_R[:, seed_number] = temp_tevl.R
tevl_Rt[:, seed_number] = temp_tevl.Rt

plot!(plot_tevl_I,
    temp_tevl.I,
    color = scenario_color[scenario_code],
    label = ""
)
plot!(plot_tevl_R,
    temp_tevl.R,
    color = scenario_color[scenario_code],
    label = ""
)
plot!(plot_tevl_Rt,
    temp_tevl.Rt,
    color = scenario_color[scenario_code],
    label = ""
)
if scenario_code == "0V" || scenario_code == "MV"
    temp_smry = CSV.read("./$scenario_code//$(lpad(seed_number, 4, '0')) smry.csv", DataFrame)
    scatter!(plot_smry_RV,
        temp_smry.R/(10^5), temp_smry.V/(10^5),
        color = scenario_color[scenario_code],
        markershape = scenario_shape[scenario_code],
        label = "",
        markerstrokewidth = 0,
        alpha = 0.8
    )
    if scenario_code == "0V"
        smry_0V[seed_number, 1] = temp_smry.R[1]/(10^5)
        smry_0V[seed_number, 2] = temp_smry.V[1]/(10^5)
    end
    if scenario_code == "MV"
        smry_MV[seed_number, 1] = temp_smry.R[1]/(10^5)
        smry_MV[seed_number, 2] = temp_smry.V[1]/(10^5)
    end
end

if seed_number == last_seed
    mean_I = mean(tevl_I, dims = 2)
    std_I = std(tevl_I, dims = 2)
    plot!(plot_trng_I, mean_I + std_I, fillrange = mean_I - std_I,
    color = scenario_color[scenario_code], linealpha = 0.2,
    label = "")

    mean_R = mean(tevl_R, dims = 2)
    std_R = std(tevl_R, dims = 2)
    plot!(plot_trng_R, mean_R + std_R, fillrange = mean_R - std_R,
    color = scenario_color[scenario_code], linealpha = 0.2,
    label = "")

    mean_Rt = mean(tevl_Rt, dims = 2)
    std_Rt = std(tevl_Rt, dims = 2)
    plot!(plot_trng_Rt, mean_Rt + std_Rt, fillrange = mean_Rt - std_Rt,
    color = scenario_color[scenario_code], linealpha = 0.2,
    label = "")
end
end

plot_smry_RV; png(plot_smry_RV, "plot_smry_RV")

plot(
    plot_tevl_I,
    plot_tevl_R,
    plot_tevl_Rt,
    layout = (3,1), size = (600,800)
); png("plot_tevl")

plot(
    plot_trng_I,
    plot_trng_R,
    plot_trng_Rt,
    layout = (3,1), size = (600,800)
); png("plot_trng")

using StatsPlots
using Statistics
boxplot(smry_0V[:,2], color = colorant"#C00000", label = "0V", size = (400,400))
boxplot!(smry_MV[:,2], color = colorant"#7030A0", label = "MV")
png("box_V")

boxplot(smry_0V[:,1], color = colorant"#C00000", label = "0V", size = (400,400))
boxplot!(smry_MV[:,1], color = colorant"#7030A0", label = "MV")
png("box_R")

plot(rand(10), linealpha = 0)

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
# A Generator represents the rules and meta information of a generator.
abstract Generator

# A choice model represents a way for how godel numbers are sampled during
# a derivation with a generator. A choice model only returns godel numbers, nothing more.
abstract ChoiceModel

choicemodel(g::Generator) = g.choicemodel
function setchoicemodel!(g::Generator, cm::ChoiceModel)
	g.choicemodel = cm
	g # output g to give useful confirmation message
end

function show(io::IO, g::Generator)
	print(io, "data generator ", typeof(g), " with ", length(choicepointinfo(g)), " choice points using ", choicemodel(g))
end

# Return the choice point info associated with a generator. This includes both the choicepointinfo for
# the generator itself plus any subgenerators
# function choicepointinfo(g::Generator)
# 	cpi = g.choicepointinfo
# 	for subgen in g.subgens
# 		merge!(cpi, subgen.choicepointinfo)
# 	end
# 	cpi
# end
# no longer accumulate choice point info from subgenerators
choicepointinfo(g::Generator) = g.choicepointinfo

# A DerivationState is created for every unique generation/derivation process. It
# collects info about the godel numbers (choices) and choice points
# as they are generated/taken. It can also provide information to the choice model
# for more complex time/state/depth related choices.
abstract DerivationState

#
# Assumed interface to all subtypes of DerivationState. Override to implement more specific
# behavior.
#
generator(s::DerivationState) = s.generator

choicemodel(s::DerivationState) = choicemodel(generator(s))

# custom exception thrown when terminating generation
type GenerationTerminatedException <: Exception
	reason::AbstractString
	GenerationTerminatedException(reason::AbstractString) = new(reason)
end
Base.showerror(io::IO, e::GenerationTerminatedException) = print(io, "generation was terminated because ", e.reason);

# determines whether generation should continue or not
# the default behaviour is to terminate, by throwing a custom error, if the limit on the length of the godel sequence
# (i.e. number of choices made) would be exceeded
# this is useful in cleaning terminating when recursive generation would lead to a extremely large or infinite object
function checkterminationcriteria(s::DerivationState)
	if length(s.godelsequence) >= s.maxchoices
		throw(GenerationTerminatedException("number of the choices made exceeded $(s.maxchoices): specify a larger value of maxchoices as a parameter to generate"))
	end
end

function logchoicepoint(s::DerivationState, cpid::Integer, godelnumber::Real, trace::Dict)
	trace[:gdl] = godelnumber
	trace[:rul] = getcurrentrulename(s)
	push!(s.cmtrace, (cpid, trace))
	push!(s.godelsequence, godelnumber)
	pushcpseqnumber(s, length(s.godelsequence)) # sequence number is the new length of the godelsequence
end

# A default derivation state that generators can use unless they need
# something more specific.
type DefaultDerivationState <: DerivationState
	generator::Generator
	godelsequence::Vector{Real} 		# Can be integers or floats
	cmtrace::Vector{Tuple{Integer,Dict}}	# choice point plus trace info returned from the choice model
	maxchoices::Int # upper limit on the size of the Godel sequence
	maxseqreps::Int # upper limit on the length of sequences from sequence choice points
	maxruledepth::Int # upper limit on rule stack depth
	rulenamestack::Vector{Symbol} # stack of called rules - most recent at the end
	cpseqnumberstack::Vector{Vector{Int}} # corresponding to the rule name stack, the seq numbers of the encountered choice points
	# executiontreecoords::Vector{Int} # during the execution of a generator, this uniquely identifies the current rule in the execution tree
	# nextchildcoord::Int
	function DefaultDerivationState(g::Generator, maxchoices::Int, maxseqreps::Int, maxruledepth::Int)
		# new(g, cm, Vector{Real}[], Vector{Integer}[], maxchoices, maxseqreps, maxruledepth, Vector{Symbol}(), Vector{Int}(), 1)
		new(g, Vector{Real}[], Vector{Integer}[], maxchoices, maxseqreps, maxruledepth, Vector{Symbol}(), Vector{Vector{Int}}())
	end
end

#
# Abstract interface to all Generator sub-types. Override for specific behavior.
#
statetype(g::Generator) = DefaultDerivationState

subgenerator(g::Generator, index::Integer) = g.subgens[index]

# Return the meta information associated with a generator.
meta(g::Generator) = g.meta

#
# Operations on derivation state related to record rule execution trees and stacks
#
# recordstartofrulemethod and recordendofrulemethod are written into rule methods by generator macro
# getcurrentrecursiondepth and getrecursiondepth can be called by, for example, a choice model (derivation state is accessible via the choice context)
#
# Note this design permits record* methods to be reasonably fast, but get* methods are slower owing to need to count.  Could denormalise by
# maintaining, for example, a dict of counts during record operations, but this would make the record methods - called by all rules - slower, while
# the get*depth methods - used only by SOME choice models - faster.
#
# "rule name stack" are the names (actually the unique method names assigned by the generator macro) of the ancestors of the current rule, including
# the currently executing rule as the last entry
# "cp seq number stack" for each rule in the stack, stores a vector of the seq nos of the choice points encountered in this rule
# this enable the ancestors (in terms of the call stack, rather than chronologically) to be determined, without storing
# large data structure
#
# OBSOLETE - replaced with cpseqnumberstack which requires less momroy - but retaining these comments along with commented out code in case another use arises
# for these coords
# "execution tree coords" are a vector of numbers, e.g. [x, y, z] which identifies the current rule as the z^th 'child' rule executed by the y^th
# child of the x^th child of the (abstract) root (x will always be 1, in fact).  Thus the immediate parent (i.e. calling) rule of a rule with coords
# [x, y, z] will have the coords [x, y]; siblings have [x, y, i!=z], descendants have [x, y, z, w, ...].
# note: the execution tree coord is unique for each *instance* of a rule being executed; the same generator rule may be executed multiple times during
# the execution of a generator, and each time it is executed, it will have different coords
#

# called at the start of generator method for each rule
recordstartofrule(s::DerivationState, rulename::Symbol) = begin
	if length(s.rulenamestack) >= s.maxruledepth
		throw(GenerationTerminatedException("rule depth exceeds the limit of $(s.maxruledepth): specify a larger value of maxruledepth as a parameter to generate"))
	end
	push!(s.rulenamestack, rulename)
	push!(s.cpseqnumberstack, Int[])
	# push!(s.executiontreecoords, s.nextchildcoord)
	# s.nextchildcoord = 1
	# println("Starting rule: $(rulename); coords: $(s.executiontreecoords); stack: $(s.rulenamestack)")
end

# called at the end of generator method for each rule
recordendofrule(s::DerivationState) = begin
	pop!(s.rulenamestack)
	# s.nextchildcoord = pop!(s.executiontreecoords) + 1
	pop!(s.cpseqnumberstack)
end

# the rule name (actually the "internal" unique method name) of the currently executing rule
# note: the unique method name is (ironically) the same for all rule with the same name (which occurs with an implicit rule choice
# in the generator) - this is achieved by having a single "umbrella" method which decides which of the rules with the same name to call
getcurrentrulename(s::DerivationState) = s.rulenamestack[end]

# the "depth" at the current rule (including the rule itself)
getruledepth(s::DerivationState) = length(s.rulenamestack)

# the recursion depth of a given rule (in terms of its immediate ancestor rules with the same name)
getrecursiondepth(s::DerivationState, rulename::Symbol) = count(r->r==rulename, s.rulenamestack)

# the recursion depth of the current rule (in terms of its immediate ancestor rules with the same name)
getcurrentrecursiondepth(s::DerivationState) = getrecursiondepth(s, getcurrentrulename(s))

# record new cp seq number
pushcpseqnumber(s::DerivationState, seqnumber::Int) = push!(s.cpseqnumberstack[end], seqnumber)

# get immediate ancestor cp seq number as the most recent seq number in the stack
getimmediateancestorcpseqnumber(s::DerivationState) = begin
	ancestorcpseqnumber = nothing
	i = length(s.cpseqnumberstack)
	while (ancestorcpseqnumber == nothing) && (i >= 1)
		if !isempty(s.cpseqnumberstack[i])
			ancestorcpseqnumber = s.cpseqnumberstack[i][end]
		end
		i -=  1
	end
	ancestorcpseqnumber
end

# # tests if execution tree coords indicate ancestry
# isexecutiontreeancestor(ancestorcoords::Vector{Int}, descendantcoords::Vector{Int}) =
# 	(length(ancestorcoords) < length(descendantcoords)) && (ancestorcoords == descendantcoords[1:length(ancestorcoords)])
#
# # tests if execution tree coords indicate ancestry or same
# isexecutiontreeselforancestor(ancestorcoords::Vector{Int}, descendantcoords::Vector{Int}) =
# 		(length(ancestorcoords) <= length(descendantcoords)) && (ancestorcoords == descendantcoords[1:length(ancestorcoords)])
#
# # tests if execution tree coords indicate parenthood
# isexecutiontreeparent(parentcoords::Vector{Int}, childcoords::Vector{Int}) = parentcoords == childcoords[1:end-1]

#
# Core functions to use generators once they have been defined.
#

# default parameters that can be overriden by user parameters to generate etc.
const defaultgenerateparams = Dict(
		:state => nothing,
		:resetchoicemodelstate => true,
		:startrule => :start,
		:maxchoices => 10017,
		:maxseqreps => 4868,
		:maxruledepth => 11765)

# Generate an object from the generator using the startrule as entry point.
# Uses the default choice model and creates a new state object unless one is given.
function generate(g::Generator; kwparams...)
	params = merge(defaultgenerateparams, Dict(p[1]=>p[2] for p in kwparams)) # kwparams takes priority
	state = (params[:state]  == nothing) ? newstate(g, params[:maxchoices], params[:maxseqreps], params[:maxruledepth]) : params[:state]
	if params[:resetchoicemodelstate]
		resetstate!(choicemodel(g))
	end
	startfunc = methodforrulenamed(g, params[:startrule])
	# important: we evaluate in the module that owns the type and since this (rather than DataGenerators) will be where rule functions are defined
	# the correct eval function is set in the generator on creation
	result = g.evalfn(Expr(:call, startfunc, g, state))
	(result, state)
end

choose(g::Generator; kwparams...) = first(generate(g; kwparams...))

function choose(gt::DataType; kwparams...)
	if !(supertype(gt) == Generator)
		error("choose(::DataType) is currently supported outside of generator code only when the datatype is a generator type")
	end
	g = gt()
	first(generate(g; kwparams...))
end

# call generator and handle termination exception
# if termination occurs, return nothing as the object
function robustchoose(g::Generator; kwparams...)
	try
		return choose(g; kwparams...)
	catch e
		if isa(e, GenerationTerminatedException)
			return nothing
		else
			throw(e) # rethrow other types of error
		end
	end
end

# This enables a rule to be called indirectly by specifying its name as a symbol (plus any parameters)
# It is used by the type generators to workaround slow compilation times that appear to be a result of recursive type inference by "hiding" the call from the type inference algorithm
function indirectrulecall(g::Generator, s::DerivationState, rulename::Symbol, ruleparams...)
	if !haskey(g.rulemethodnames, rulename)
		throw(GenerationTerminatedException("no rule exists for requested name: $(rulename) in rule(...)"))
	end
	g.evalfn(Expr(:call, g.rulemethodnames[rulename], g, s, ruleparams...))
end

#
# Core choice point functions that have no short forms and should never be written directly in user code.
#

# choose a number (used by value choice points with numeric types: choose(Bool|Int|Float64|...))
function choosenumber(s::DerivationState, cpid, datatype, minval, maxval, paramsliteral)
	if !paramsliteral
		try
			minval = convert(datatype, minval)
		catch(InexactError)
			error("minimum value cannot be converted to a $(datatype) in choose_number")
		end
		try
			maxval = convert(datatype, maxval)
		catch(InexactError)
			error("maximum value cannot be converted to a $(datatype) in choose_number")
		end
		if maxval < minval
			error("maximum value is less than minimum in choose_number")
		end
	end
	querychoicemodel(s, :value, cpid, datatype, minval, maxval)
end

# choose a number of repetitions (used by sequence choice points from constructs reps,mult,plus)
function choosereps(s::DerivationState, cpid, minreps, maxreps, paramsliteral)
	if !paramsliteral
		try
			minreps = convert(Int, minreps)
		catch(InexactError)
			error("minimum repetitions is not an integer in choose_reps")
		end
		try
			maxreps = convert(Int, maxreps)
		catch(InexactError)
			error("maximum repetitions is not an integer in choose_reps")
		end
		if minreps < 0
			error("minimum repetitions is less than zero in choose_reps")
		end
		if maxreps < minreps
			error("maximum repetitions is less than minimum in choose_reps")
		end
	end
	reps = querychoicemodel(s, :sequence, cpid, Int, minreps, maxreps)
	if reps > s.maxseqreps
		throw(GenerationTerminatedException("a sequence choice point made more than $(reps) repetitions and so exceeded the limit of $(s.maxseqreps): specify a larger value of maxseqrep as a parameter to generate"))
		# warn("choice model specified $(reps) repetitions for a sequence choice point, but this is being reduced to the maximum of $(s.maxseqreps) - specify a larger value of maxseqreps as a parameter to generate if this is not the required behaviour")
		# reps = s.maxseqreps
	end
	reps
end

# The (implicit) rule choice is implemented with this function. The short form is that
# the user writes rulename() at least two times in the generator code.
function chooserule(s::DerivationState, cpid, numrules)
	querychoicemodel(s, :rule, cpid, Int, 1, numrules)
end

#
# Queries to the choice model
#

# ChoiceContext holds info about the current choice point such as its type and id number, and constraints such as the
# minimum and maximum values that should be returned.  This information should be used by the choice model to return
# an appropriate Godel number.
type ChoiceContext
	derivationstate::DerivationState
	cptype::Symbol
	cpid::UInt64
	datatype::DataType
	lowerbound::Real
	upperbound::Real
end

# All queries to choice model to retrieve Godel numbers made via this function
# sets up ChoiceContext that may be used by choice model
function querychoicemodel(s::DerivationState, cptype, cpid, datatype, lowerbound, upperbound)
	# check if generation should continue or not - function raises an error if not
	checkterminationcriteria(s)
	# choicecontext = ChoiceContext(s, cptype, cpid, datatype, lowerbound, upperbound, 0)
	choicecontext = ChoiceContext(s, cptype, cpid, datatype, lowerbound, upperbound)
	# TODO recursiondepth
	gn, trace = godelnumber(choicemodel(s), choicecontext)
	@assert lowerbound <= gn <= upperbound
	logchoicepoint(s, cpid, gn, trace)
	# convert the godel number to the specified type
	# this is most useful for choose_number to ensure returned value is of specified type
	# note that for Bool, convert treats 0 as false, and 1 as true
	convert(datatype, gn)
end

# Generate data from the subgenerator with given index.
subgen(g::Generator, s::DerivationState, index::Integer) = choose(subgenerator(g, index))

#
# Internal functions used by the core functions. They are not exported.
#

# Return the top-level method for thhe rule
# Since rule method have unique names to avoid issues with extending existing methods, need to look this up
function methodforrulenamed(g::Generator, rulename::Symbol)
	if !haskey(g.rulemethodnames, rulename)
		error("generator has no rule named $(rulename)")
	end
  g.rulemethodnames[rulename]
end

# Create a new derivation state object for a given generator and choice model.
function newstate(g::Generator, maxchoices::Int, maxseqreps::Int, maxruledepth::Int)
	st = statetype(g)
	st(g, maxchoices, maxseqreps, maxruledepth)
end


text:
export server, ServerActor

using Sockets

"""
    ServerActor{D, Address, Port}() where D

The `ServerActor` sends all `next!`/`error!`/`complete!` events to the local network listeners
with specified `Address` and `Port` parameters via `TCPSocket`.

See also: [`Actor`](@ref), [`server`](@ref)
"""
struct ServerActor{D, A, P} <: Actor{D}
    server  :: Sockets.TCPServer
    sockets :: Vector{Sockets.TCPSocket}

    ServerActor{D, A, P}() where { D, A, P } = begin
        self = new(listen(A, P), Vector{Sockets.TCPSocket}())

        @async begin
            while true
                listener = accept(self.server)
                push!(self.sockets, listener)
                filter!(socket -> isopen(socket), self.sockets)
            end
        end

        return self
    end
end

function on_next!(actor::ServerActor, data)
    filter!(socket -> __send_next(socket, data), actor.sockets)
end

function on_error!(actor::ServerActor, err)
    foreach(socket -> __send_error(socket, err), actor.sockets)
    close(actor.server)
end

function on_complete!(actor::ServerActor)
    foreach(socket -> __send_complete(socket), actor.sockets)
    close(actor.server)
end

function __send_next(socket, data::D) where D
    if isopen(socket)
        try
            return write(socket, 1) !== 0 && write(socket, data) !== 0
        catch _
            return false
        end
    end
    return false
end

function __send_next(socket, data::Vector{D}) where D
    if isopen(socket)
        try
            size = length(data)
            return write(socket, size) !== 0 && write(socket, data) !== 0
        catch _
            return false
        end
    end
    return false
end

function __send_error(socket, err)
    if isopen(socket)
        try
            return write(socket, -1) !== 0
        catch _
            return false
        end
    end
    return false
end

function __send_complete(socket)
    if isopen(socket)
        try
            return write(socket, 0) !== 0
        catch _
            return false
        end
    end
    return false
end

struct ServerActorFactory{Address, Port} <: AbstractActorFactory end

create_actor(::Type{L}, factory::ServerActorFactory{Address, Port}) where { L, Address, Port } = server(L, Address, Port)

"""
    server(port::Int)
    server(address::A, port::Int) where { A <: IPAddr }
    server(::Type{D}, port::Int)
    server(::Type{D}, address::A, port::Int) where { A <: IPAddr }

Creation operator for the `ServerActor` actor.

See also: [`AbstractActor`](@ref)
"""
server(port::Int)                                   = server(Sockets.localhost, port)
server(address::A, port::Int) where { A <: IPAddr } = ServerActorFactory{address, port}()

function server(::Type{D}, port::Int) where D
    return server(D, Sockets.localhost, port)
end

function server(::Type{D}, address::A, port::Int) where { D, A <: IPAddr }
    @assert isbits(zero(D)) "Network server actor supports only primitive data types"
    return ServerActor{D, address, port}()
end

function server(::Type{Vector{D}}, address::A, port::Int) where { D, A <: IPAddr }
    @assert isbits(zero(D)) "Network server actor supports only primitive data types"
    return ServerActor{Vector{D}, address, port}()
end


text:
#!/usr/bin/env julia --project=.

InFile1="test-dataset/scenario_1/case.con"
InFile2="test-dataset/scenario_1/case.inl"
InFile3="test-dataset/scenario_1/case.raw"
InFile4="test-dataset/scenario_1/case.rop"
TimeLimitInSeconds=600
ScoringMethod=2
NetworkModel="IEEE 14"

println("  $(InFile1)")
println("  $(InFile2)")
println("  $(InFile3)")
println("  $(InFile4)")
println("  $(TimeLimitInSeconds)")
println("  $(ScoringMethod)")
println("  $(NetworkModel)")

include("MyJulia2.jl")

compute_solution2(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds, ScoringMethod, NetworkModel, output_dir="test-dataset/scenario_1")


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

using Test

using PyCall

@testset "PyCall" begin
    @info "PyCall uses python command \"$(PyCall.python)\""
end


text:
# Constructs a Differentiable Optimizer model from a MOI Optimizer.
# Supports `forward_differentiate!` and `reverse_differentiate!` methods for solving and differentiating the model respectectively.

# ## Note
# Currently supports differentiating linear and quadratic programs only.

const MOIDD = MOI.Utilities.DoubleDicts

Base.@kwdef mutable struct InputCache
    dx::Dict{VI, Float64} = Dict{VI, Float64}()# dz for QP
    # ds
    # dy #= [d\lambda, d\nu] for QP
    # FIXME Would it be possible to have a DoubleDict where the value depends
    #       on the function type ? Here, we need to create two dicts to have
    #       concrete value types.
    # `scalar_constraints` and `vector_constraints` includes `A` and `b` for CPs
    # or `G` and `h` for QPs
    scalar_constraints::MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.ScalarAffineFunction{Float64}}() # also includes G for QPs
    vector_constraints::MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}} = MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}}() # also includes G for QPs
    objective::Union{Nothing,MOI.AbstractScalarFunction} = nothing
end

function Base.empty!(cache::InputCache)
    empty!(cache.dx)
    empty!(cache.scalar_constraints)
    empty!(cache.vector_constraints)
    cache.objective = nothing
    return
end

"""
    reverse_differentiate!(model::MOI.ModelLike)

Wrapper method for the backward pass / reverse differentiation.
This method will consider as input a currently solved problem and differentials
with respect to the solution set with the [`ReverseVariablePrimal`](@ref) attribute.
The output problem data differentials can be queried with the
attributes [`ReverseObjectiveFunction`](@ref) and [`ReverseConstraintFunction`](@ref).
"""
function reverse_differentiate! end

"""
    forward_differentiate!(model::Optimizer)

Wrapper method for the forward pass.
This method will consider as input a currently solved problem and
differentials with respect to problem data set with
the [`ForwardObjectiveFunction`](@ref) and  [`ForwardConstraintFunction`](@ref) attributes.
The output solution differentials can be queried with the attribute
[`ForwardVariablePrimal`](@ref).
"""
function forward_differentiate! end

"""
    ForwardObjectiveFunction <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to set input data to forward differentiation, that
is, problem input data.
The possible values are any `MOI.AbstractScalarFunction`.
A `MOI.ScalarQuadraticFunction` can only be used in linearly constrained
quadratic models.

For instance, if the objective contains `θ * (x + 2y)`, for the purpose of
computing the derivative with respect to `θ`, the following should be set:
```julia
MOI.set(model, DiffOpt.ForwardObjectiveFunction(), 1.0 * x + 2.0 * y)
```
where `x` and `y` are the relevant `MOI.VariableIndex`.
"""
struct ForwardObjectiveFunction <: MOI.AbstractModelAttribute end

"""
    ForwardConstraintFunction <: MOI.AbstractConstraintAttribute

A `MOI.AbstractConstraintAttribute` to set input data to forward differentiation, that
is, problem input data.

For instance, if the scalar constraint of index `ci` contains `θ * (x + 2y) <= 5θ`,
for the purpose of computing the derivative with respect to `θ`, the following
should be set:
```julia
MOI.set(model, DiffOpt.ForwardConstraintFunction(), ci, 1.0 * x + 2.0 * y - 5.0)
```
Note that we use `-5` as the `ForwardConstraintFunction` sets the tangent of the
ConstraintFunction so we consider the expression `θ * (x + 2y - 5)`.
"""
struct ForwardConstraintFunction <: MOI.AbstractConstraintAttribute end

"""
    ForwardVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to get output data from forward
differentiation, that is, problem solution.

For instance, to get the tangent of the variable of index `vi` corresponding to
the tangents given to `ForwardObjectiveFunction` and `ForwardConstraintFunction`, do the
following:
```julia
MOI.get(model, DiffOpt.ForwardVariablePrimal(), vi)
```
"""
struct ForwardVariablePrimal <: MOI.AbstractVariableAttribute end
MOI.is_set_by_optimize(::ForwardVariablePrimal) = true

"""
    ReverseVariablePrimal <: MOI.AbstractVariableAttribute

A `MOI.AbstractVariableAttribute` to set input data to
reverse differentiation, that is, problem solution.

For instance, to set the tangent of the variable of index `vi`, do the
following:
```julia
MOI.set(model, DiffOpt.ReverseVariablePrimal(), x)
```
"""
struct ReverseVariablePrimal <: MOI.AbstractVariableAttribute end

"""
    ReverseObjectiveFunction <: MOI.AbstractModelAttribute

A `MOI.AbstractModelAttribute` to get output data to reverse differentiation,
that is, problem input data.

For instance, to get the tangent of the objective function corresponding to
the tangent given to `ReverseVariablePrimal`, do the
following:
```julia
func = MOI.get(model, DiffOpt.ReverseObjectiveFunction())
```
Then, to get the sensitivity of the linear term with variable `x`, do
```julia
JuMP.coefficient(func, x)
```
To get the sensitivity with respect to the quadratic term with variables `x`
and `y`, do either
```julia
JuMP.coefficient(func, x, y)
```
or
```julia
DiffOpt.quad_sym_half(func, x, y)
```

!!! warning
    These two lines are **not** equivalent in case `x == y`, see
    [`quad_sym_half`](@ref) for the details on the difference between these two
    functions.
"""
struct ReverseObjectiveFunction <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::ReverseObjectiveFunction) = true

"""
    ReverseConstraintFunction

An `MOI.AbstractConstraintAttribute` to get output data to reverse differentiation, that
is, problem input data.

For instance, if the following returns `x + 2y + 5`, it means that the tangent
has coordinate `1` for the coefficient of `x`, coordinate `2` for the
coefficient of `y` and `5` for the function constant.
If the constraint is of the form `func == constant` or `func <= constant`,
the tangent for the constant on the right-hand side is `-5`.
```julia
MOI.get(model, DiffOpt.ReverseConstraintFunction(), ci)
```
"""
struct ReverseConstraintFunction <: MOI.AbstractConstraintAttribute end
MOI.is_set_by_optimize(::ReverseConstraintFunction) = true

"""
    abstract type AbstractModel <: MOI.ModelLike end

Model supporting [`forward_differentiate!`](@ref) and
[`reverse_differentiate!`](@ref).
"""
abstract type AbstractModel <: MOI.ModelLike end

MOI.supports_incremental_interface(::AbstractModel) = true

MOI.is_valid(model::AbstractModel, idx::MOI.Index) = MOI.is_valid(model.model, idx)

function MOI.add_variable(model::AbstractModel)
    return MOI.add_variable(model.model)
end

function MOI.add_variables(model::AbstractModel, n)
    return MOI.add_variables(model.model, n)
end

function MOI.Utilities.pass_nonvariable_constraints(
    dest::AbstractModel,
    src::MOI.ModelLike,
    idxmap::MOIU.IndexMap,
    constraint_types,
)
    MOI.Utilities.pass_nonvariable_constraints(dest.model, src, idxmap, constraint_types)
end

function MOI.Utilities.final_touch(model::AbstractModel, index_map)
    MOI.Utilities.final_touch(model.model, index_map)
end

function MOI.supports_constraint(model::AbstractModel, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    return MOI.supports_constraint(model.model, F, S)
end

function MOI.add_constraint(model::AbstractModel, func::MOI.AbstractFunction, set::MOI.AbstractSet)
    return MOI.add_constraint(model.model, func, set)
end

function _enlarge_set(vec::Vector, idx, value)
    m = last(idx)
    if length(vec) < m
        n = length(vec)
        resize!(vec, m)
        fill!(view(vec, (n+1):m), NaN)
    end
    vec[idx] = value
    return
end

function MOI.set(model::AbstractModel, ::MOI.VariablePrimalStart, vi::MOI.VariableIndex, value)
    MOI.throw_if_not_valid(model, vi)
    _enlarge_set(model.x, vi.value, value)
end

function MOI.set(model::AbstractModel, ::ForwardObjectiveFunction, objective)
    model.input_cache.objective = objective
    return
end
function MOI.set(model::AbstractModel, ::ReverseVariablePrimal, vi::VI, val)
    model.input_cache.dx[vi] = val
    return
end
function MOI.set(model::AbstractModel,
    ::ForwardConstraintFunction,
    ci::CI{MOI.ScalarAffineFunction{T},S},
    func::MOI.ScalarAffineFunction{T},
) where {T,S}
    model.input_cache.scalar_constraints[ci] = func
    return
end
function MOI.set(model::AbstractModel,
    ::ForwardConstraintFunction,
    ci::CI{MOI.VectorAffineFunction{T},S},
    func::MOI.VectorAffineFunction{T},
) where {T,S}
    model.input_cache.vector_constraints[ci] = func
    return
end

function lazy_combination(op::F, α, a, β, b) where {F<:Function}
    return LazyArrays.ApplyArray(
        op,
        LazyArrays.@~(α .* a),
        LazyArrays.@~(β .* b),
    )
end

function lazy_combination(op::F, α, a, β, b, I::UnitRange) where {F<:Function}
    return lazy_combination(op, α, view(a, I), β, view(b, I))
end
function lazy_combination(op::F, a, b, i::Integer, args::Vararg{Any,N}) where {F<:Function,N}
    return lazy_combination(op, a[i], b, b[i], a, args...)
end
function lazy_combination(op::F, a, b, i::UnitRange, I::UnitRange) where {F<:Function}
    return lazy_combination(op, view(a, i), b, view(b, i), a, I)
end

_lazy_affine(vector, constant::Number) = VectorScalarAffineFunction(vector, constant)
_lazy_affine(matrix, vector) = MatrixVectorAffineFunction(matrix, vector)
function _get_db end
function _get_dA end
function MOI.get(model::AbstractModel, ::ReverseConstraintFunction, ci::CI)
    return _lazy_affine(_get_dA(model, ci), _get_db(model, ci))
end

"""
    π(v::Vector{Float64}, model::MOI.ModelLike, cones::ProductOfSets)

Given a `model`, its `cones`, find the projection of the vectors `v`
of length equal to the number of rows in the conic form onto the cartesian
product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function π(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets) where T
    return map_rows(model, cones, Flattened{T}()) do ci, r
        MOSD.projection_on_set(
            MOSD.DefaultDistance(),
            v[r],
            MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci))
        )
    end
end

"""
    Dπ(v::Vector{Float64}, model, cones::ProductOfSets)

Given a `model`, its `cones`, find the gradient of the projection of
the vectors `v` of length equal to the number of rows in the conic form onto the
cartesian product of the cones corresponding to these rows.
For more info, refer to https://github.com/matbesancon/MathOptSetDistances.jl
"""
function Dπ(v::Vector{T}, model::MOI.ModelLike, cones::ProductOfSets) where T
    return BlockDiagonals.BlockDiagonal(
        map_rows(model, cones, Nested{Matrix{T}}()) do ci, r
            MOSD.projection_gradient_on_set(
                MOSD.DefaultDistance(),
                v[r],
                MOI.dual_set(MOI.get(model, MOI.ConstraintSet(), ci)),
            )
        end
    )
end

# See the docstring of `map_rows`.
struct Nested{T} end
struct Flattened{T} end

# Store in `x` the values `y` corresponding to the rows `r` and the `k`th
# constraint.
function _assign_mapped!(x, y, r, k, ::Nested)
    x[k] = y
end
function _assign_mapped!(x, y, r, k, ::Flattened)
    x[r] = y
end

# Map the rows corresponding to `F`-in-`S` constraints and store it in `x`.
function _map_rows!(f::Function, x::Vector, model, cones::ProductOfSets, ::Type{F}, ::Type{S}, map_mode, k) where {F, S}
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        r = MOI.Utilities.rows(cones, ci)
        k += 1
        _assign_mapped!(x, f(ci, r), r, k, map_mode)
    end
    return k
end

# Allocate a vector for storing the output of `map_rows`.
_allocate_rows(cones, ::Nested{T}) where {T} = Vector{T}(undef, length(cones.dimension))
_allocate_rows(cones, ::Flattened{T}) where {T} = Vector{T}(undef, MOI.dimension(cones))

"""
    map_rows(f::Function, model, cones::ProductOfSets, map_mode::Union{Nested{T}, Flattened{T}})

Given a `model`, its `cones` and `map_mode` of type `Nested` (resp.
`Flattened`), return a `Vector{T}` of length equal to the number of cones (resp.
rows) in the conic form where the value for the index (resp. rows) corresponding
to each cone is equal to `f(ci, r)` where `ci` is the corresponding constraint
index in `model` and `r` is a `UnitRange` of the corresponding rows in the conic
form.
"""
function map_rows(f::Function, model, cones::ProductOfSets, map_mode::Union{Nested, Flattened})
    x = _allocate_rows(cones, map_mode)
    k = 0
    for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
        # Function barrier for type unstability of `F` and `S`
        # `con_map` is a `MOIU.DoubleDicts.MainIndexDoubleDict`, we index it at `F, S`
        # which returns a `MOIU.DoubleDicts.IndexWithType{F, S}` which is type stable.
        # If we have a small number of different constraint types and many
        # constraint of each type, this mostly removes type unstabilities
        # as most the time is in `_map_rows!` which is type stable.
        k = _map_rows!(f, x, model, cones, F, S, map_mode, k)
    end
    return x
end

function _fill(neg::Function, gradient_cache, input_cache, cones, args...)
    _fill(S -> true, neg, gradient_cache, input_cache, cones, args...)
end
function _fill(filter::Function, neg::Function, gradient_cache, input_cache, cones, args...)
    for (F, S) in keys(input_cache.scalar_constraints.dict)
        filter(S) || continue
        _fill(args..., neg(S), cones, input_cache.scalar_constraints[F,S])
    end
    for (F, S) in keys(input_cache.vector_constraints.dict)
        _fill(args..., neg(S), cones, input_cache.vector_constraints[F,S])
    end
    return
end

function _fill(vector::Vector, neg::Bool, cones, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, ci)
        vector[r] = neg ? -MOI.constant(func) : MOI.constant(func)
    end
end
function _fill(I::Vector, J::Vector, V::Vector, neg::Bool, cones, dict)
    for (ci, func) in dict
        r = MOI.Utilities.rows(cones, ci)
        for term in func.terms
            _push_term(I, J, V, neg, r, term)
        end
    end
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::Integer, term::MOI.ScalarAffineTerm)
    push!(I, r)
    push!(J, term.variable.value)
    push!(V, neg ? -term.coefficient : term.coefficient)
end
function _push_term(I::Vector, J::Vector, V::Vector, neg::Bool, r::UnitRange, term::MOI.VectorAffineTerm)
    _push_term(I, J, V, neg, r[term.output_index], term.scalar_term)
end

function MOI.supports(model::AbstractModel, attr::MOI.AbstractModelAttribute)
    return MOI.supports(model.model, attr)
end

function MOI.set(model::AbstractModel, attr::MOI.AbstractModelAttribute, value)
    MOI.set(model.model, attr, value)
end

function MOI.get(model::AbstractModel, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.model, attr)
end


text:
using HDF5, Serialization

file = h5open("imaging_data/7_20X.h5", "r")

raw_data = Array(file["data"]["neurons"])

# Convert to column major and reinterpret values as UInt16 (orginal form)
raw_data = UInt16.(permutedims(raw_data, reverse(1:ndims(raw_data))))

# use native Julia serializer so we can grab the stack from disk faster while developing
serialize("imaging_data/neurondata.blob", raw_data)

max_stars_repo_name: wsphillips/jlcaimaging

max_issues_repo_name: wsphillips/jlcaimaging

max_forks_repo_name: wsphillips/jlcaimaging