text:
using Mdp.VI
@testset "Dict inputs" begin
    # Define simple domains (from Massau Kolobov 2016 book, p.41)
    S = 1:6 # States
    A = 1:2
    T = Dict(   # Transitions, prob = T[(state, action, next_state)]
        (1,1,3)=>1.,
        (1,2,2)=>1.,
        (2,1,3)=>1.,
        (3,1,5)=>1.,
        (3,2,2)=>1.,
        (4,1,5)=>1.,
        (5,1,6)=>1.,
        (5,2,6)=>.6,
        (5,2,4)=>.4,
        (6,1,6)=>1.,    # Manually added goal state transition
    )
    C = Dict( # Costs, cost = C[(state, action, next_state)]
        (1,1,3)=>1,
        (1,2,2)=>1,
        (2,1,3)=>1,
        (3,1,5)=>1,
        (3,2,2)=>1,
        (4,1,5)=>1,
        (5,1,6)=>5,
        (5,2,6)=>2,
        (5,2,4)=>2,
        (6,1,6)=>0,   # Manually added goal state transition cost
    )

    V_init = [3.,3.,2.,2.,1.,0.]
    @testset "Bellman backup" begin
        @test bellman_backup(1, S, A, T, C, V_init)==3.
        @test bellman_backup(2, S, A, T, C, V_init)==3.
        @test bellman_backup(3, S, A, T, C, V_init)==2.
        @test bellman_backup(4, S, A, T, C, V_init)==2.
        @test bellman_backup(5, S, A, T, C, V_init)==2.8
        @test bellman_backup(6, S, A, T, C, V_init)==0.
    end

    @test isapprox(value_iteration(S,A,T,C, ϵ=0.001, V_old=[3.,3.,2.,2.,1.,0.]), [5.99921,5.99921,4.99969, 4.99969,3.99969,0.], atol=0.001)
end

@testset "Function inputs" begin
    # Define simple domains (from Massau Kolobov 2016 book, p.41)
    S = 1:6 # States
    A = 1:2
    T = Dict(   # Transitions, prob = T[(state, action, next_state)]
        (1,1,3)=>1.,
        (1,2,2)=>1.,
        (2,1,3)=>1.,
        (3,1,5)=>1.,
        (3,2,2)=>1.,
        (4,1,5)=>1.,
        (5,1,6)=>1.,
        (5,2,6)=>.6,
        (5,2,4)=>.4,
        (6,1,6)=>1.,    # Manually added goal state transition
    )
    C = Dict( # Costs, cost = C[(state, action, next_state)]
        (1,1,3)=>1,
        (1,2,2)=>1,
        (2,1,3)=>1,
        (3,1,5)=>1,
        (3,2,2)=>1,
        (4,1,5)=>1,
        (5,1,6)=>5,
        (5,2,6)=>2,
        (5,2,4)=>2,
        (6,1,6)=>0,   # Manually added goal state transition cost
    )

    # V_init = [3.,3.,2.,2.,1.,0.]
    # @testset "Bellman backup" begin
    #     @test bellman_backup(1, S, A, T, C, V_init)==3.
    #     @test bellman_backup(2, S, A, T, C, V_init)==3.
    #     @test bellman_backup(3, S, A, T, C, V_init)==2.
    #     @test bellman_backup(4, S, A, T, C, V_init)==2.
    #     @test bellman_backup(5, S, A, T, C, V_init)==2.8
    #     @test bellman_backup(6, S, A, T, C, V_init)==0.
    # end
    #
    # @test isapprox(value_iteration(S,A,T,C, ϵ=0.001, V_old=[3.,3.,2.,2.,1.,0.]), [5.99921,5.99921,4.99969, 4.99969,3.99969,0.], atol=0.001)
end


text:
using WebAPI
using Test

@testset "WebAPI.jl" begin
    # Write your tests here.
end

ext: jl


text:
"""
Add observational noise to a dataset from a normal distribution with zero mean
and standard deviation equivalent to some fraction of the empirical standard
deviation along each variable of the dataset.
"""
function add_observational_noise!(x::AbstractArray{T, 2}, percent_noise) where T

    dims = size(x)
    if percent_noise > 0
        if dims[1] > dims[2]
            nvars = dims[2]
            npts = dims[1]

            for i = 1:nvars
                dist = Normal(0, (percent_noise/100) * std(x[:, i]))
                x[:, i] .+= rand(dist, npts)
            end
        else
            nvars = dims[1]
            npts = dims[2]

            for i = 1:nvars
                dist = Normal(0, (percent_noise/100) * std(x[i, :]))
                x[i, :] .+= rand(dist, npts)
            end
        end
    end

    return x
end

function add_observational_noise!(data::AbstractDataset{dim, T}, percent_noise) where {dim, T}
    D = dim
    L = length(data)
    noisy_data = zeros(T, L, D)
    if percent_noise > 0
        σs = [Normal(0, (percent_noise/100)*std(col)) for col in columns(data)]
        for i = 1:D
            noisy_data[:, i] = data[:, i] .+ rand(σs[i], L)
        end
        return Dataset(noisy_data)
    end
    return data
end


text:
@testset "Georeferencing" begin
  table = DataFrame(x=rand(3), y=[1,2,3], z=["a","b","c"])
  tuple = (x=rand(3), y=[1,2,3], z=["a","b","c"])

  # explicit domain types
  d = georef(table, PointSet(rand(2,3)))
  @test domain(d) isa PointSet
  d = georef(tuple, PointSet(rand(2,3)))
  @test domain(d) isa PointSet
  d = georef(table, RegularGrid(3))
  @test domain(d) isa RegularGrid
  d = georef(tuple, RegularGrid(3))
  @test domain(d) isa RegularGrid

  # coordinates of point set
  d = georef(table, rand(2,3))
  @test domain(d) isa PointSet
  d = georef(tuple, rand(2,3))
  @test domain(d) isa PointSet

  # coordinates names in table
  d = georef(table, (:x,:y))
  @test domain(d) isa PointSet

  # regular grid data
  d = georef(tuple)
  @test domain(d) == RegularGrid(3)
  tuple2D = (x=rand(10,10), y=rand(Int,10,10))
  d = georef(tuple2D)
  @test domain(d) == RegularGrid(10,10)
  d = georef(tuple2D, (1.,2.), (3.,4.))
  @test domain(d) == RegularGrid((10,10), (1.,2.), (3.,4.))
  d = georef(tuple2D, spacing=(3.,4.))
  @test domain(d) == RegularGrid((10,10), (0.,0.), (3.,4.))
  tuple3D = (x=rand(10,10,10), y=rand(10,10,10))
  d = georef(tuple3D)
  @test domain(d) == RegularGrid(10,10,10)
end


text:
@with_kw mutable struct ProgressiveHedgingData{T <: AbstractFloat}
    Q::T = 1e10
    δ::T = 1.0
    δ₁::T = 1.0
    δ₂::T = 1.0
    iterations::Int = 0
end

@with_kw mutable struct ProgressiveHedgingParameters{T <: AbstractFloat}
    τ::T = 1e-5
    log::Bool = true
    keep::Bool = true
    offset::Int = 0
    indent::Int = 0
end

"""
    ProgressiveHedging

Functor object for the progressive-hedging algorithm. Create using the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Parameters
- `τ::AbstractFloat = 1e-6`: Relative tolerance for convergence checks.
- `log::Bool = true`: Specifices if progressive-hedging procedure should be logged on standard output or not.
...
"""
struct ProgressiveHedging{T <: AbstractFloat,
                          A <: AbstractVector,
                          SP <: StochasticProgram,
                          S <: LQSolver,
                          PT <: PenaltyTerm,
                          E <: AbstractExecution,
                          P <: AbstractPenalization} <: AbstractProgressiveHedgingSolver
    stochasticprogram::SP
    data::ProgressiveHedgingData{T}
    parameters::ProgressiveHedgingParameters{T}

    # Estimate
    ξ::A
    Q_history::A
    dual_gaps::A

    # Subproblems
    nscenarios::Int

    # Execution
    execution::E

    penalization::P

    # Params
    progress::ProgressThresh{T}

    function ProgressiveHedging(stochasticprogram::StochasticProgram,
                                x₀::AbstractVector,
                                subsolver::MPB.AbstractMathProgSolver,
                                executer::Execution,
                                penalizer::AbstractPenalizer,
                                penaltyterm::PenaltyTerm; kw...)
        if nworkers() > 1 && executer isa Serial
            @warn "There are worker processes, consider using distributed version of algorithm"
        end
        executer = if nworkers() == 1 && !(executer isa Serial)
            @warn "There are no worker processes, defaulting to serial version of algorithm"
            Serial()
        else
            executer
        end
        first_stage = StochasticPrograms.get_stage_one(stochasticprogram)
        length(x₀) != first_stage.numCols && error("Incorrect length of starting guess, has ", length(x₀), " should be ", first_stage.numCols)

        T = promote_type(eltype(x₀), Float32)
        x₀_ = convert(AbstractVector{T}, copy(x₀))
        A = typeof(x₀_)
        SP = typeof(stochasticprogram)
        S = LQSolver{typeof(MPB.LinearQuadraticModel(subsolver)), typeof(subsolver)}
        PT = typeof(penaltyterm)
        n = StochasticPrograms.nscenarios(stochasticprogram)
        execution = executer(T,A,S,PT)
        E = typeof(execution)
        penalization = penalizer()
        P = typeof(penalization)
        params = ProgressiveHedgingParameters{T}(; kw...)

        ph = new{T,A,SP,S,PT,E,P}(stochasticprogram,
                                  ProgressiveHedgingData{T}(),
                                  params,
                                  x₀_,
                                  A(),
                                  A(),
                                  n,
                                  execution,
                                  penalization,
                                  ProgressThresh(1.0, 0.0, "$(indentstr(params.indent))Progressive Hedging"))
        # Initialize solver
        initialize!(ph, subsolver, penaltyterm)
        return ph
    end
end
ProgressiveHedging(stochasticprogram::StochasticProgram,
                   subsolver::MPB.AbstractMathProgSolver,
                   execution::Execution,
                   penalizer::AbstractPenalizer,
                   penaltyterm::PenaltyTerm; kw...) = ProgressiveHedging(stochasticprogram,
                                                                         rand(decision_length(stochasticprogram)),
                                                                         subsolver,
                                                                         execution,
                                                                         penalizer,
                                                                         penaltyterm; kw...)

function (ph::ProgressiveHedging)()
    # Reset timer
    ph.progress.tfirst = ph.progress.tlast = time()
    # Start workers (if any)
    start_workers!(ph)
    # Start procedure
    while true
        status = iterate!(ph)
        if status != :Valid
            close_workers!(ph)
            return status
        end
    end
end


text:
using DiffEqBase, DiffEqJump

using Plots;
doPlot = true
using BenchmarkTools
doBenchMark = false

Nsims = 8000

# average number of proteins in a simulation
function runSSAs(jump_prob)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[3,end]
    end
    mean(Psamp)
end

# DNA repression model DiffEqBiological
# using DiffEqBiological
# rs = @reaction_network dtype begin
#     k1, DNA --> mRNA + DNA
#     k2, mRNA --> mRNA + P
#     k3, mRNA --> 0
#     k4, P --> 0
#     k5, DNA + P --> DNAR
#     k6, DNAR --> DNA + P
# end k1 k2 k3 k4 k5 k6

# model using mass action jumps
# ids: DNA=1, mRNA = 2, P = 3, DNAR = 4
reactstoch =
[
    [1 => 1],
    [2 => 1],
    [2 => 1],
    [3 => 1],
    [1 => 1, 3 => 1],
    [4 => 1]
]
netstoch =
[
    [2 => 1],
    [3 => 1],
    [2 => -1],
    [3 => -1],
    [1 => -1, 3 => -1, 4 => 1],
    [1 => 1, 3 => 1, 4 => -1]
]
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
majumps = MassActionJumps(rates, reactstoch, netstoch)

tf     = 20000.0
prob   = DiscreteProblem([1,0,0,0], (0.0, tf), rates)

# SSAs to test
SSAalgs = (Direct(), DirectFW(), FRM(), FRMFW())

# plotting one full trajectory
if doPlot
    plothand = plot()
    for alg in SSAalgs
        jump_prob = JumpProblem(prob, alg, majumps)
        sol = solve(jump_prob, SSAStepper())
        plot!(plothand, sol.t, sol[3,:])
    end
    display(plothand)
end

# benchmark performance
if doBenchMark
    # exact methods
    for alg in SSAalgs
        println("Solving with method: ", typeof(method), ", using SSAStepper")
        jump_prob = JumpProblem(prob, method, rs)
        @btime solve($jump_prob, SSAStepper())
    end
    println()
end

# write your own tests here
#@test 1 == 1

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_read.jl - https://github.com/strickek/ReadWriteDlm2.jl

"""

    readdlm2(source; options...)
    readdlm2(source, T::Type; options...)
    readdlm2(source, delim::AbstractChar; options...)
    readdlm2(source, delim::AbstractChar, T::Type; options...)
    readdlm2(source, delim::AbstractChar, eol::AbstractChar; options...)
    readdlm2(source, delim::AbstractChar, T::Type, eol::AbstractChar; options...)

Read a matrix from `source`. The `source` can be a text file, stream or byte
array. Each line (separated by `eol`, this is `'\\n'` by default) gives one row.
The columns are separated by `';'`, another `delim` can be defined.

Pre-processing of `source` with regex substitution changes the decimal marks
from `d,d` to `d.d`. For default `rs` the keyword argument `decimal=','` sets
the decimal Char in the `r`-string of `rs`. When a special regex substitution
tuple `rs=(r.., s..)` is defined, the argument `decimal` is not used.
Pre-processing can be switched off with: `rs=()`.

In addition to stdlib readdlm(), data is also parsed for `Dates` formats,
the `Time` format `\"HH:MM[:SS[.s{1,9}]]\"` and for complex and rational numbers.
To deactivate parsing dates/time set: `dfs=\"\", dtfs=\"\"`.
`locale` defines the language of day (`E`, `e`) and month (`U`, `u`) names.

The result will be a (heterogeneous) array of default element type `Any`. If
`header=true` it will be a tuple containing the data array and a vector for
the columnnames. Other (abstract) types for the data array elements could be
defined. If data is empty, a `0×0 Array{T,2}` is returned.

With `tables=true`[, `header=true`] option[s] a `Tables` interface compatible
`MatrixTable` with individual column types is returned, which for example
can be used as Argument for `DataFrame()`.

# Additional Keyword Arguments

* `decimal=','`: Decimal mark Char used by default `rs`, irrelevant if `rs`-tuple is not the default one
* `rs=(r\"(\\d),(\\d)\", s\"\\1.\\2\")`: Regex (r,s)-tuple, the default change d,d to d.d if `decimal=','`
* `dtfs=\"yyyy-mm-ddTHH:MM:SS.s\"`: Format string for DateTime parsing
* `dfs=\"yyyy-mm-dd\"`: Format string for Date parsing
* `locale=\"english\"`: Language for parsing dates names
* `tables=false`: Return `Tables` interface compatible MatrixTable if `true`
* `dfheader=false`: 'dfheader=true' is shortform for `tables=true, header=true`
* `missingstring=\"na\"`: How missing values are represented

Find more information about `readdlm()` functionality and (keyword) arguments -
which are also supported by `readdlm2()` - in `help` for `readdlm()`.

# Code Example
```jldoctest
julia> using ReadWriteDlm2

julia> A = Any[1 1.2; "text" missing];

julia> writedlm2("test.csv", A)

julia> readdlm2("test.csv")
2×2 Array{Any,2}:
 1        1.2
  "text"   missing
```
"""
readdlm2(input; opts...) =
    readdlm2auto(input, ';', Nothing, '\n', true; opts...)

readdlm2(input, T::Type; opts...) =
    readdlm2auto(input, ';', T, '\n', false; opts...)

readdlm2(input, dlm::AbstractChar; opts...) =
    readdlm2auto(input, dlm, Nothing, '\n', true; opts...)

readdlm2(input, dlm::AbstractChar, T::Type; opts...) =
    readdlm2auto(input, dlm, T, '\n', false; opts...)

readdlm2(input, dlm::AbstractChar, eol::AbstractChar; opts...) =
    readdlm2auto(input, dlm, Nothing, eol, true; opts...)

readdlm2(input, dlm::AbstractChar, T::Type, eol::AbstractChar; opts...) =
    readdlm2auto(input, dlm, T, eol, false; opts...)

function readdlm2auto(input, dlm, T, eol, auto;
    decimal::AbstractChar=',',
    rs::Tuple=(r"(\d),(\d)", s"\1.\2"),
    dtfs::AbstractString="yyyy-mm-ddTHH:MM:SS.s",
    dfs::AbstractString="yyyy-mm-dd",
    locale::AbstractString="english",
    tables::Bool=false,
    dfheader::Bool=false,
    missingstring::AbstractString="na",
    opts...)

    if dfheader == true
        tables = true
        opts = [opts...]
        push!(opts, :header => true)
    end

    ((!isempty(dtfs) && !occursin(Regex("[^YymdHMSs]"), dtfs)) ||
    (!isempty(dfs) && !occursin(Regex("[^YymdHMSs]"), dfs))) && info(
        """
        Format string for DateTime(`$dtfs`) or Date(`$dfs`)
        contains numeric code elements only. At least one non-numeric
        code element or character is needed for parsing dates.
        """)

    # "parsing-logic" - defaults
    doparsedatetime = !isempty(dtfs)
    doparsedate = !isempty(dfs)
    doparsetime = false
    doparsecomplex = true
    doparserational = true
    doparsemissing = false
    doparsenothing = false
    convertarray = true
    T2 = Any

    if T == Nothing            # to know wether T ist Any by default
            T = Any
            anybydefault = true
        else
            anybydefault = false
        end

    # "parsing-logic" for different T::Types
    ((typeintersect(DateTime, T)) == Union{}) && (doparsedatetime = false)
    ((typeintersect(Date, T)) == Union{}) && (doparsedate = false)
    doparsetime =
        ((doparsedatetime && doparsedate) || ((Time <: T) && !(Any <:T)))
    T <: Union{AbstractFloat, AbstractString, Char} &&
        (T2 = T; convertarray = false)
    ((typeintersect(Complex, T)) == Union{}) && (doparsecomplex = false)
    ((typeintersect(Rational, T)) == Union{}) && (doparserational = false)
    Missing <: T && (doparsemissing = true)
    Nothing <: T && (doparsenothing = true)
    (Any <: T) && (convertarray = false)
    if tables == true
        T2 = Any
        convertarray = false
        !anybydefault && (convertarray = true)
        T <: Union{AbstractString, Char} &&
            (T2 = T; convertarray = false)
    end

    s = read(input, String)

    # empty input data - return empty array
    if (isempty(s) || (s == string(eol)))
        return Array{T2}(undef, 0, 0)
    end

    if (!isempty(rs) && (decimal != '.')) # do pre-processing of decimal mark

        # adopt dfs if decimal between two digits of formatstring
        if (rs == (r"(\d),(\d)", s"\1.\2")) && doparsedate
            drs = (Regex("([YymdHM])$decimal([YymdHM])"), s"\1.\2")
            dfs = replace(dfs, drs[1] => drs[2])
        end

        # adopt dtfs if decimal between two digits of formatstring
        if (rs == (r"(\d),(\d)", s"\1.\2")) && doparsedatetime
            drs = (Regex("([YymdHMSs])$decimal([YymdHMSs])"), s"\1.\2")
            dtfs = replace(dtfs, drs[1] => drs[2])
        end

        # Change default regex substitution Tuple if decimal != ','
        if ((rs == (r"(\d),(\d)", s"\1.\2")) && (decimal != ','))
            rs = (Regex("(\\d)$decimal(\\d)"), s"\1.\2")
        end

        # Error if decimal mark to replace is also the delim Char
        "1"*string(dlm)*"1" != replace("1"*string(dlm)*"1", rs[1] => rs[2]) &&
        error(
            """
            Error: Decimal mark to replace is also the delim Char.
            Pre-processing with decimal mark Regex substitution for
            `$(dlm)` (= delim!!) is not allowed - change rs/decimal or delim!
            """)

        # Regex substitution decimal
        s = replace(s, rs[1] => rs[2])

    end

    # Using stdlib DelimitedFiles internal function to read dlm-string
    z = readdlm_string(s, dlm, T2, eol, auto, val_opts(opts))

    # Formats, Regex for Date/DateTime parsing
    doparsedatetime &&
        (dtdf = DateFormat(dtfs, locale); rdt = dfregex(dtfs, locale))
    doparsedate &&
        (ddf = DateFormat(dfs, locale); rd = dfregex(dfs, locale))

    if tables == true  # return MatrixTable for Tables Interface
        if isa(z, Tuple)
             y, h = z
             headerexist = true
             cn = Symbol.(reshape(h, :))
        else
             y = z
             headerexist = false
        end

        rows, cols = size(y)
        coltypes = Array{Type, 1}(undef, cols)
        Tn = typeintersect(T, Number)
        for c in 1:cols # iterate columns
            colcontainmissing = false
            colcontainnothing = false
            coltype = Union{}
        for r in 1:rows # iterate rows in columns

            if isa(y[r,c], AbstractString)
                if doparsedatetime && occursin(rdt, y[r,c]) # parse DateTime
                    try y[r,c] = DateTime(y[r,c], dtdf) catch; end
                elseif doparsedate && occursin(rd, y[r,c]) # parse Date
                    try y[r,c] = Date(y[r,c], ddf) catch; end
                elseif doparsemissing && y[r,c] == missingstring # parse Missing
                    try y[r,c] = missing catch; end
                elseif doparsenothing && y[r,c] == "nothing" # parse Nothing
                    try y[r,c] = nothing catch; end
                else # parse Time, Complex and Rational
                    try y[r,c] = parseothers(y[r,c], doparsetime, doparsecomplex, doparserational) catch; end
                    if isa(y[r,c], AbstractString)
                        #Substring to String
                        try y[r,c] = String.(split(y[r,c]))[1] catch; end
                    end
                end
            elseif convertarray && !(typeof(y[r,c]) <: Tn)
                try y[r,c] = convert(Tn, y[r,c]) catch; end
            end
            if typeof(y[r,c]) == Missing
                colcontainmissing = true
            elseif typeof(y[r,c]) == Nothing
                colcontainnothing = true
            else
                coltype = typejoin(coltype, typeof(y[r,c]))
            end

        end # end for - iterate rows in columns

        if colcontainmissing && colcontainnothing
            coltypes[c] = Union{Missing, Nothing, coltype}
        elseif colcontainmissing
            coltypes[c] = Union{Missing, coltype}
        elseif colcontainnothing
            coltypes[c] = Union{Nothing, coltype}
        else
            coltypes[c] = coltype
        end

        if convertarray && (coltypes[c] <: T)
            coltypes[c] = T
        end

        end # end for - iterate columns

        vdata = vecofvec(y, coltypes)

        if headerexist
            return Tables.table(vdata, header=cn)
        else
            return Tables.table(vdata)
        end

    else     # return standard Format (Arry or Tuple(data, header))

        if isa(z, Tuple)
             y, h = z
        else
             y = z
        end
        for i in eachindex(y)
            if isa(y[i], AbstractString)
                if doparsedatetime && occursin(rdt, y[i]) # parse DateTime
                    try y[i] = DateTime(y[i], dtdf) catch; end
                elseif doparsedate && occursin(rd, y[i]) # parse Date
                    try y[i] = Date(y[i], ddf) catch; end
                elseif doparsemissing && y[i] == missingstring # parse Missing
                    try y[i] = missing catch; end
                elseif doparsenothing && y[i] == "nothing" # parse Nothing
                    try y[i] = nothing catch; end
                else # parse Time, Complex and Rational
                    try y[i] = parseothers(y[i], doparsetime, doparsecomplex, doparserational) catch; end
                    end
                end
            end
        end

        if convertarray
            isa(z, Tuple) ? z = (convert(Array{T}, y), h) : z = convert(Array{T}, z)
        end

        return z

end # End function readdlm2auto()


text:
# Main data structures for the solvers.
include("ds.jl")
include("ds_accessors.jl")
include("none.jl")

# Modelling helpers.
include("_helpers_inflow.jl")
include("_helpers_output.jl")
include("_helpers_bathymetry.jl")
include("_helpers_purpose.jl")

# Evaluation solvers.
include("evaluate_crossvalidate.jl")
include("evaluate_deviation.jl")
include("evaluate_mc_feasibility.jl")
include("evaluate_purposeshortage.jl")
include("evaluate_hydropower.jl")

# Inverse solvers.
include("inverse_functions.jl")

# Rule curve solvers.
include("rulecurve_min_default.jl")
include("rulecurve_min_safe.jl")


text:
using SeaPearl: DefaultStateRepresentation, VariableVertex, ValueVertex

function plottripartite(sr::DefaultStateRepresentation)
    cpmodel = sr.cplayergraph
    n = cpmodel.totalLength
    nodefillc = []
    for id in 1:n
        v = cpmodel.idToNode[id]
        if isa(v, VariableVertex) push!(nodefillc,"red")
        elseif isa(v, ValueVertex) push!(nodefillc,"blue")
        else  push!(nodefillc,"black") end
    end
    gplot(cpmodel;nodefillc=nodefillc)
end

export plottripartite


text:
"""
    DTCSeparable{Tz<:AbstractVector, Tk<:SeparableKernel} <: Kernel

Specifies a low-rank approximation to a kernel `k` through pseudo-inputs `z`. `z` are the
locations of the pseudo-inputs in _space_, since they are automatically replicated
(implicitly) at each point in time.
"""
struct DTCSeparable{Tz<:AbstractVector, Tk<:Separable} <: Kernel
    z::Tz
    k::Tk
end

"""
    dtcify(z::AbstractVector, x)

Recurse through `x` and replace any `Separable` kernels with a `DTCSeparable`. The resulting
kernel can be used to construct the approximate model utilised by the DTC, and is used to
compute the ELBO.
"""
dtcify(z::AbstractVector, k::Separable) = DTCSeparable(z, k)

dtcify(z::AbstractVector, k::ScaledKernel) = ScaledKernel(dtcify(z, k.kernel), k.σ²)

function dtcify(z::AbstractVector, k::TransformedKernel{<:Kernel, <:ScaleTransform})
    return TransformedKernel(dtcify(z, k.kernel), k.transform)
end

function dtcify(z::AbstractVector, k::KernelSum)
    return KernelSum(dtcify(z, k.kernels[1]), dtcify(z, k.kernels[2]))
end

dtcify(z::AbstractVector, fx::FiniteLTISDE) = FiniteGP(dtcify(z, fx.f), fx.x, fx.Σy)

dtcify(z::AbstractVector, fx::LTISDE) = LTISDE(dtcify(z, fx.f), fx.storage)

dtcify(z::AbstractVector, f::GP) = GP(f.mean, dtcify(z, f.kernel))

"""
    dtc(fx::FiniteLTISDE, y::AbstractVector{<:Real}, z_r::AbstractVector)

Compute the DTC (Deterministic Training Conditional) in state-space form [insert reference].

`fx` and `y` are the same as would be provided to `logpdf`, and `z_r` is a specification of
the spatial location of the pseudo-points at each point in time.

Note that this API is slightly different from AbstractGPS.jl's API, in which `z_r` is
replaced by a `FiniteGP`.

WARNING: this API is unstable, and subject to change in future versions of TemporalGPs. It
was thrown together quickly in pursuit of a conference deadline, and has yet to receive the
attention it deserves.
"""
function dtc(fx::FiniteLTISDE, y::AbstractVector, z_r::AbstractVector)
    return logpdf(dtcify(z_r, fx), y)
end

# This stupid pullback saves an absurb amount of compute time.
function Zygote._pullback(::AContext, ::typeof(count), ::typeof(ismissing), yn)
    return count(ismissing, yn), nograd_pullback
end

"""
    elbo(fx::FiniteLTISDE, y::AbstractVector{<:Real}, z_r::AbstractVector)

Compute the ELBO (Evidence Lower BOund) in state-space form [insert reference].
"""
function AbstractGPs.elbo(fx::FiniteLTISDE, y::AbstractVector, z_r::AbstractVector)

    fx_dtc = dtcify(z_r, fx)

    # Compute diagonals over prior marginals.
    lgssm = build_lgssm(fx_dtc)
    Σs = lgssm.emissions.fan_out.Q
    marg_diags = marginals_diag(lgssm)

    k = fx_dtc.f.f.kernel
    Cf_diags = kernel_diagonals(k, fx_dtc.x)

    # Transform a vector into a vector-of-vectors.
    y_vecs = restructure(y, lgssm.emissions)

    tmp = zygote_friendly_map(
        ((Σ, Cf_diag, marg_diag, yn), ) -> begin
            Σ_, _ = fill_in_missings(Σ, yn)
            return sum(diag(Σ_ \ (Cf_diag - marg_diag.P))) -
                count(ismissing, yn) + size(Σ_, 1)
        end,
        zip(Σs, Cf_diags, marg_diags, y_vecs),
    )

    return logpdf(lgssm, y_vecs) - sum(tmp) / 2
end

Zygote.accum(x::NamedTuple{(:diag, )}, y::Diagonal) = Zygote.accum(x, (diag=y.diag, ))

function kernel_diagonals(k::DTCSeparable, x::RectilinearGrid)
    space_kernel = k.k.l
    time_kernel = k.k.r
    Cr_rpred_diag = kernelmatrix_diag(space_kernel, get_space(x))
    time_vars = kernelmatrix_diag(time_kernel, get_times(x))
    return map(s_t -> Diagonal(Cr_rpred_diag * s_t), time_vars)
end

function kernel_diagonals(k::DTCSeparable, x::RegularInTime)
    space_kernel = k.k.l
    time_kernel = k.k.r
    time_vars = kernelmatrix_diag(time_kernel, get_times(x))
    return map(
        (s_t, x_r) -> Diagonal(kernelmatrix_diag(space_kernel, x_r) * s_t),
        time_vars,
        x.vs,
    )
end

function kernel_diagonals(k::ScaledKernel, x::AbstractVector)
    return k.σ²[1] .* kernel_diagonals(k.kernel, x)
end

function kernel_diagonals(k::KernelSum, x::AbstractVector)
    return kernel_diagonals(k.kernels[1], x) .+ kernel_diagonals(k.kernels[2], x)
end

function lgssm_components(k_dtc::DTCSeparable, x::SpaceTimeGrid, storage::StorageType)

    # Construct temporal model.
    k = k_dtc.k
    ts = get_times(x)
    time_kernel = k.r
    As_t, as_t, Qs_t, emission_proj, x0_t = lgssm_components(time_kernel, ts, storage)
    Hs_t, hs_t = _extract_emission_proj(emission_proj)

    # Compute spatial covariance between inducing inputs, and inducing points + obs. points.
    space_kernel = k.l
    x_space = x.xl
    z_space = k_dtc.z
    K_space_z = kernelmatrix(space_kernel, z_space)
    K_space_zx = kernelmatrix(space_kernel, z_space, x_space)

    # Get some size info.
    M = length(z_space)
    N = length(x_space)
    ident_M = my_I(eltype(storage), M)

    # G is the time-invariant component of the H-matrices. It is only time-invariant because
    # we have the same obsevation locations at each point in time.
    Λu_Cuf = cholesky(Symmetric(K_space_z + 1e-12I)) \ K_space_zx

    # Construct approximately low-rank model spatio-temporal LGSSM.
    As = map(A -> kron(ident_M, A), As_t)
    as = map(a -> repeat(a, M), as_t)
    Qs = map(Q -> kron(K_space_z, Q), Qs_t)
    Cs = Fill(Λu_Cuf, length(ts))
    cs = map(h -> Fill(h, N), hs_t) # This should currently be zero.
    Hs = map(H -> kron(ident_M, H), Hs_t)
    hs = Fill(Zeros(M), length(ts))
    x0 = Gaussian(repeat(x0_t.m, M), kron(K_space_z, x0_t.P))
    return As, as, Qs, (Cs, cs, Hs, hs), x0
end

function lgssm_components(k_dtc::DTCSeparable, x::RegularInTime, storage::StorageType)

    # Construct temporal model.
    k = k_dtc.k
    ts = get_times(x)
    time_kernel = k.r
    As_t, as_t, Qs_t, emission_proj, x0_t = lgssm_components(time_kernel, ts, storage)
    Hs_t, hs_t = _extract_emission_proj(emission_proj)

    # Compute spatial covariance between inducing inputs, and inducing points + obs. points.
    space_kernel = k.l
    z_space = k_dtc.z
    K_space_z = kernelmatrix(space_kernel, z_space)
    K_space_z_chol = cholesky(Symmetric(K_space_z + 1e-9I))

    # Get some size info.
    M = length(z_space)
    N = length(ts)
    ident_M = my_I(eltype(storage), M)

    # Construct approximately low-rank model spatio-temporal LGSSM.
    As = zygote_friendly_map(
        ((I, A), ) -> kron(I, A),
        zip(Fill(ident_M, N), As_t),
    )
    as = zygote_friendly_map(a -> repeat(a, M), as_t)
    Qs = zygote_friendly_map(
        ((K_space_z, Q), ) -> kron(K_space_z, Q),
        zip(Fill(K_space_z, N), Qs_t),
    )
    x_big = _reduce(vcat, x.vs)
    C__ = kernelmatrix(space_kernel, z_space, x_big)
    C = \(K_space_z_chol, C__)
    Cs = partition(Zygote.dropgrad(map(length, x.vs)), C)

    cs = map((h, v) -> fill(h, length(v)), hs_t, x.vs) # This should currently be zero.
    Hs = zygote_friendly_map(
        ((I, H_t), ) -> kron(I, H_t),
        zip(Fill(ident_M, N), Hs_t),
    )
    hs = Fill(Zeros(M), N)
    x0 = Gaussian(repeat(x0_t.m, M), kron(K_space_z, x0_t.P))

    return As, as, Qs, (Cs, cs, Hs, hs), x0
end

_extract_emission_proj((Hs, hs)::Tuple{AbstractVector, AbstractVector}) = Hs, hs

_reduce(::typeof(vcat), xs::Vector{<:Vector{<:Real}}) = reduce(vcat, xs)

function _reduce(::typeof(vcat), xs::Vector{<:ColVecs})
    return ColVecs(reduce(hcat, getfield.(xs, :X)))
end

function partition(lengths::AbstractVector{<:Integer}, A::Matrix{<:Real})
    starts = vcat(1, cumsum(lengths) .+ 1)
    starts = starts[1:end-1]
    return map((s, d) -> collect(view(A, :, s:s+d-1)), starts, lengths)
end

function Zygote._pullback(
    ctx::AContext,
    ::typeof(partition),
    lengths::AbstractVector{<:Integer},
    A::Matrix{<:Real},
)
    partition_pullback(::Nothing) = nothing
    partition_pullback(Δ::Vector) = nothing, nothing, reduce(hcat, Δ)
    return partition(lengths, A), partition_pullback
end

function build_emissions(
    (Cs, cs, Hs, hs)::Tuple{AbstractVector, AbstractVector, AbstractVector, AbstractVector},
    Σs::AbstractVector,
)
    Hst = map(adjoint, Hs)
    Cst = map(adjoint, Cs)
    fan_outs = StructArray{LargeOutputLGC{eltype(Cs), eltype(cs), eltype(Σs)}}((Cst, cs, Σs))
    return StructArray{BottleneckLGC{eltype(Hst), eltype(hs), eltype(fan_outs)}}((Hst, hs, fan_outs))
end

"""
    approx_posterior_marginals(
        ::typeof(dtc),
        fx::FiniteLTISDE,
        y::AbstractVector,
        z_r::AbstractVector,
        x_r::AbstractVector,
    )

Compute the DTC (Deterministic Training Conditional) approximation to the posterior
marginals at the times provided by `fx`, but at the new spatial locations `x_r`, given
observations `y` and spatial pseudo-input locations `z_r`.

WARNING: this API is unstable, and subject to change in future versions of TemporalGPs. It
was thrown together quickly in pursuit of a conference deadline, and has yet to receive the
attention it deserves.
"""
function approx_posterior_marginals(
    ::typeof(dtc),
    fx::FiniteLTISDE,
    y::AbstractVector,
    z_r::AbstractVector,
    x_r::AbstractVector,
)
    fx.f.f.mean isa AbstractGPs.ZeroMean || throw(error("Prior mean of GP isn't zero."))

    # Compute approximate posterior LGSSM.
    lgssm = build_lgssm(dtcify(z_r, fx))
    fx_post = posterior(lgssm, restructure(y, lgssm.emissions))

    # Compute the new emission distributions + approx posterior model.
    x_pr = RectilinearGrid(x_r, get_times(fx.x))
    k_dtc = dtcify(z_r, fx.f.f.kernel)
    new_proj, Σs = dtc_post_emissions(k_dtc, x_pr, fx.f.storage)
    new_fx_post = LGSSM(fx_post.transitions, build_emissions(new_proj, Σs))

    # Compute marginals under modified posterior.
    return reduce(vcat, map(marginals, marginals_diag(new_fx_post)))
end

"""
    approx_posterior_marginals(
        ::typeof(dtc),
        fx::FiniteLTISDE,
        y::AbstractVector,
        z_r::AbstractVector,
        x_r::AbstractVector,
        t::Int,
    )

Same as other method, but only returns the predictions at index `t` in `fx`.

As with the other method of this function, it's a bit of a hack. It's correct of course, but
needs to be tidied up at some point.
"""
function approx_posterior_marginals(
    ::typeof(dtc),
    fx::FiniteLTISDE,
    y::AbstractVector,
    z_r::AbstractVector,
    x_r::AbstractVector,
    t::Int,
)
    ts = get_times(fx.x)
    if t < 1 || t > length(ts)
        throw(error("t = $t must be between 1 and length(ts) = $(length(ts))."))
    end

    # Compute approximate posterior LGSSM.
    lgssm = build_lgssm(dtcify(z_r, fx))
    fx_post = posterior(lgssm, restructure(y, lgssm.emissions))

    # Prep prediction locations.
    x_other = x_r[1:1]
    x_rs = fill(x_other, length(ts))
    x_rs[t] = x_r
    x_pr = RegularInTime(ts, x_rs)

    # Compute the new emission distributions + approx posterior model.
    k_dtc = dtcify(z_r, fx.f.f.kernel)
    new_proj, Σs = dtc_post_emissions(k_dtc, x_pr, fx.f.storage)
    new_fx_post = LGSSM(fx_post.transitions, build_emissions(new_proj, Σs))

    # Compute marginals under modified posterior.
    return marginals(marginals_diag(new_fx_post)[t])
end

function approx_posterior_marginals(
    ::typeof(dtc),
    fx::FiniteLTISDE,
    y::AbstractVector,
    z_r::AbstractVector,
    x_pr::RegularInTime,
)
    ts = get_times(fx.x)
    if ts != get_times(x_pr)
        throw(error("Times don't match."))
    end

    # Compute approximate posterior LGSSM.
    lgssm = build_lgssm(dtcify(z_r, fx))
    fx_post = posterior(lgssm, restructure(y, lgssm.emissions))

    # Compute the new emission distributions + approx posterior model.
    k_dtc = dtcify(z_r, fx.f.f.kernel)
    new_proj, Σs = dtc_post_emissions(k_dtc, x_pr, fx.f.storage)
    new_fx_post = LGSSM(fx_post.transitions, build_emissions(new_proj, Σs))

    # Compute marginals under modified posterior.
    return reduce(vcat, map(marginals, marginals_diag(new_fx_post)))
end

function build_emission_covs(k::DTCSeparable, x_new::RectilinearGrid)
    space_kernel = k.k.l
    z_r = k.z
    C_fp_u = kernelmatrix(space_kernel, get_space(x_new), z_r)
    C_u = cholesky(Symmetric(kernelmatrix(space_kernel, z_r) + ident_eps(z_r, 1e-9)))
    Cr_rpred_diag = kernelmatrix_diag(space_kernel, get_space(x_new))
    spatial_Q_diag = Cr_rpred_diag - diag_Xt_invA_X(C_u, C_fp_u')

    time_kernel = k.k.r
    time_vars = kernelmatrix_diag(time_kernel, get_times(x_new))
    return map(s_t -> Diagonal(spatial_Q_diag * s_t), time_vars)
end

function build_emission_covs(k::DTCSeparable, x_new::RegularInTime)
    space_kernel = k.k.l
    z_r = k.z
    C_u = cholesky(Symmetric(kernelmatrix(space_kernel, z_r) + ident_eps(z_r, 1e-9)))

    time_kernel = k.k.r
    time_vars = kernelmatrix_diag(time_kernel, get_times(x_new))
    return map(zip(time_vars, x_new.vs)) do ((time_var, x_r))
        C_fp_u = kernelmatrix(space_kernel, x_r, z_r)
        Cr_rpred_diag = kernelmatrix_diag(space_kernel, x_r)
        spatial_Q_diag = Cr_rpred_diag - diag_Xt_invA_X(C_u, C_fp_u')
        return Diagonal(spatial_Q_diag * time_var)
    end
end

function dtc_post_emissions(k::DTCSeparable, x_new::AbstractVector, storage::StorageType)
    _, _, _, new_proj, _ = lgssm_components(k, x_new, storage)
    new_Σs = build_emission_covs(k, x_new)
    return new_proj, new_Σs
end

function dtc_post_emissions(k::ScaledKernel, x_new::AbstractVector, storage::StorageType)
    (Cs, cs, Hs, hs), Σs = dtc_post_emissions(k.kernel, x_new, storage)
    σ = sqrt(convert(eltype(storage_type), only(k.σ²)))
    return (Cs, cs, map(H->σ * H, Hs), map(h->σ * h, hs)), map(Σ->σ^2 * Σ, Σs)
end

function dtc_post_emissions(k::KernelSum, x_new::AbstractVector, storage::StorageType)
    (Cs_l, cs_l, Hs_l, hs_l), Σs_l = dtc_post_emissions(k.kernels[1], x_new, storage)
    (Cs_r, cs_r, Hs_r, hs_r), Σs_r = dtc_post_emissions(k.kernels[2], x_new, storage)
    Cs = map(vcat, Cs_l, Cs_r)
    cs = cs_l + cs_r
    Hs = map(blk_diag, Hs_l, Hs_r)
    hs = map(vcat, hs_l, hs_r)
    return (Cs, cs, Hs, hs), map(+, Σs_l, Σs_r)
end


text:
@testset "Test Propagation of the Lorenz equation I" begin
    u0 = [10.0; -5.0; 2.0]
    t0 = 0.0
    tf = 5.0
    Δt = 1e-2

    J = ceil(Int64, (tf-t0)/Δt)

    h(u, t) = [sum(u)]
    F = StateSpace(lorenz63!, h)

    prob = ODEProblem(F.f,u0,(t0,tf))
    sol = solve(prob, RK4(), dt = Δt, adaptive = false);

#     Run it for the different ensemble members
    Nx = 3
    Ny = 1
    Ne = 5
    X = zeros(Ny+Nx, Ne)

    viewstate(X, Ny, Nx) .= repeat(u0, 1, Ne)

    statehist = Array{Float64,2}[]
    push!(statehist, deepcopy(X[Ny+1:Ny+Nx,:]))

    for i = 1:J
        # Forecast
        tspan = (t0+(i-1)*Δt, t0+i*Δt)
        prob = remake(prob; tspan=tspan)

        prob_func(prob,i,repeat) = ODEProblem(prob.f,view(X,Ny+1:Ny+Nx,i),prob.tspan)

        ensemble_prob = EnsembleProblem(prob,output_func = (sol,i) -> (sol[end], false),
        prob_func=prob_func)
        sim = solve(ensemble_prob, RK4(), dt = Δt, adaptive = false, EnsembleThreads(),trajectories = Ne,
        dense = false, save_everystep=false);

        @inbounds for i=1:Ne
            view(X, Ny+1:Ny+Nx,i) .= sim[i]
        end

        push!(statehist, deepcopy(viewstate(X, Ny, Nx)))
    end

    observe(F.h, X, 5.0, Ny, Nx)

    for i=1:Ne
        @test norm(view(statehist[end], :,i) - sol(5.0))<1e-11
        @test norm(view(X, 1:Ny,i) - h(sol(5.0), 5.0))<1e-11
        @test norm(view(X, Ny+1:Ny+Nx,i) - sol(5.0))<1e-11

    end
end


text:
using BinDeps

@BinDeps.setup

version = v"0.21.2"
libgit2 = library_dependency("libgit2",
    validate = (name, handle) -> begin
        major, minor, patch = Cint[0], Cint[0], Cint[0]
        ccall(dlsym(handle, :git_libgit2_version), Void,
              (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, patch)
        v = VersionNumber(major[1], minor[1])
        if v.major == version.major && v.minor == version.minor
            return true
        else
            return false
        end
    end
)

provides(Sources,
         URI("https://github.com/libgit2/libgit2/archive/v$version.tar.gz"),
         libgit2,
         unpacked_dir="libgit2-$version")

# Windows binaries built in MSYS2 via:
# /d/code/CMake-2.8/bin/cmake .. -DCMAKE_INSTALL_PREFIX=$PWD/usr -G"MSYS Makefiles"
# /d/code/CMake-2.8/bin/cmake --build . --target install
provides(Binaries, URI("http://sourceforge.net/projects/juliadeps-win/files/libgit2-v$version.7z"),
         libgit2, unpacked_dir="usr$WORD_SIZE/bin", os = :Windows)

prefix = joinpath(BinDeps.depsdir(libgit2),"usr")
srcdir = joinpath(BinDeps.depsdir(libgit2),"src","libgit2-$version")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libgit2)
        @build_steps begin
            ChangeDirectory(srcdir)
            `cmake . -DCMAKE_INSTALL_PREFIX=$prefix
                     -DCMAKE_BUILD_TYPE=Release
                     -DTHREADSAFE=ON
                     -DBUILD_CLAR=OFF`
            `make install`
        end
    end), libgit2, os = :Unix)

@BinDeps.install [:libgit2 => :libgit2]

cd(Pkg.dir("LibGit2")) do
    run(`git submodule init`)
    run(`git submodule update`)
end


text:
###################################################
## Use timings to compute shortest paths
###################################################

"""
Compute shortest path using the given timings
"""
function shortestPaths!(r::RoutingPaths)
    g = r.network.graph
    if !is_strongly_connected(g)
        error("The network needs to be strongly connected")
    end

    pathTimes = zeros(nv(g), nv(g))
    parents  = zeros(Int,nv(g),nv(g))

    for orig in 1:nv(g)
        d = dijkstra_shortest_paths(g, orig, r.times)
        pathTimes[orig,:] = d.dists
        parents[orig,:] = d.parents
    end
    r.pathTimes=pathTimes
    r.pathPrevious=parents
    return r
end

"""
Compute shortest path using the given timings
- parallel version, uses all available processors
- if batch_number is provided, use pmap method
"""
function parallelShortestPaths!(r::RoutingPaths; batch_number::Int=0)
    if batch_number > 0
        return parallelShortestPaths2(r,batch_number)
    end

    if nprocs() <= 1
        println("Only one thread: use shortestPaths! instead")
    end

    g = r.network.graph
    if !is_strongly_connected(g)
        error("The network needs to be strongly connected")
    end

    pathTimes = SharedArray(Float64, (nv(g), nv(g)))
    parents  = SharedArray(Int,(nv(g),nv(g)))

    @sync @parallel for orig in 1:nv(g)
        d = dijkstra_shortest_paths(g, orig, r.times)
        pathTimes[orig,:] = d.dists
        parents[orig,:] = d.parents
    end
    r.pathTimes=pathTimes
    r.pathPrevious=parents
    return r
end

"""
 another version, using pmap after grouping : usualy slower for now
"""
function parallelShortestPaths2!(r::RoutingPaths, batch_number::Int=10)

    if nworkers() <= 1
        println("Only one thread: use shortestPaths! instead")
    end
    g = r.network.graph
    if !is_strongly_connected(g)
        error("The network needs to be strongly connected")
    end

    pathTimes = SharedArray(Float64, (nv(g), nv(g)))
    parents  = SharedArray(Int,(nv(g),nv(g)))

    @sync for w in workers()
        @spawnat(w, eval(Main, Expr(:(=), :r, r)))
    end

    function range_dijkstra(range)
        for orig in range
            d = dijkstra_shortest_paths(r.network.graph, orig, r.times)
            pathTimes[orig,:] = d.dists
            parents[orig,:] = d.parents
        end
    end
    steps = round(Int,linspace(1,nv(g), batch_number+1))
    ranges = [steps[i]:steps[i+1] for i in 1:(length(steps)-1)]

    pmap(range_dijkstra, ranges)

    r.pathTimes=sdata(pathTimes)
    r.pathPrevious=sdata(parents)
    return r
end

"""
    walking time between i and j, doesn't use routing paths,
    applies dijkstra to undirected version of graph
"""
function walkingtime(n::Network, i::Int, j::Int,
                     times::AbstractArray{Float64,2} = walkingDistances(n))
    d = dijkstra_shortest_paths(Graph(n.graph), i, times)
    return d.dists[j]
end

"""
    returns nodes that are within provided walkingTime of node i, along
    with all walking times computed by dijkstra
    doesn't use routing paths and applies dijkstra to undirected version of graph
"""
function walkingNodes(n::Network, i::Int, walkingTime::Float64,
                      times::AbstractArray{Float64,2} = walkingDistances(n))
    d = dijkstra_shortest_paths(Graph(n.graph), i, times)
    nodes = find(d.dists .< walkingTime)
    return nodes, d.dists[nodes]
end

"""
    Using Dijkstra, returns the all pair shortest path for a subset of the nodes as
    origins and destinations
"""
function allPairTimes(n::Network, origins::Vector{Int}, dests::Vector{Int}, times::AbstractArray{Float64, 2})
    result = Matrix{Float64}(length(origins), length(dests))
    for (i,o) in enumerate(origins)
        d = dijkstra_shortest_paths(n.graph, o, times)
        result[i,:] = d.dists[dests]
    end
    return result
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module RelicToolkit

module Config
    macro enabled(id)
        env = split(get(ENV, "RELIC_TOOLKIT_CURVE", id), ",")
        return id in env ? :(true) : :(false)
    end
end

module BN254
    using ..Config: @enabled
    const ENABLED = @enabled("BN254")
    if ENABLED
        const LIB = "librelic_gmp_pbc_bn254"

        include(joinpath(@__DIR__, "Model.jl"))
        include(joinpath(@__DIR__, "LowLevel.jl"))
        include(joinpath(@__DIR__, "HighLevel.jl"))

        export curve_miller, field_final_exp
        curve_miller(a::EP, b::EP2) = pp_map_oatep_k12!(FP12(undef), a, b)
        curve_miller(a::EP2, b::EP) = pp_map_oatep_k12!(FP12(undef), b, a)
        field_final_exp(a::FP12) = pp_exp_k12!(FP12(undef), a)
    end
end

module BLS381
    using ..Config: @enabled
    const ENABLED = @enabled("BLS381")
    if ENABLED
        const LIB = "librelic_gmp_pbc_bls381"

        include(joinpath(@__DIR__, "Model.jl"))
        include(joinpath(@__DIR__, "LowLevel.jl"))
        include(joinpath(@__DIR__, "HighLevel.jl"))

        export curve_miller, field_final_exp
        curve_miller(a::EP, b::EP2) = pp_map_oatep_k12!(FP12(undef), a, b)
        curve_miller(a::EP2, b::EP) = pp_map_oatep_k12!(FP12(undef), b, a)
        field_final_exp(a::FP12) = pp_exp_k12!(FP12(undef), a)
    end
end

end # module


text:
module test_server

using HTTP, HTTP.IOExtras, Sockets, Test, MbedTLS

function testget(url, m=1)
    r = []
    @sync for i in 1:m
        l = rand([0,0,10,1000,10000])
        body = Vector{UInt8}(rand('A':'Z', l))
        # println("sending request...")
        @async push!(r, HTTP.request("GET", "$url/$i", [], body))
    end
    return r
end

@testset "HTTP.listen" begin
    port = 8087 # rand(8000:8999)

    # echo response
    handler = (http) -> begin
        request::HTTP.Request = http.message
        request.body = read(http)
        closeread(http)
        request.response::HTTP.Response = HTTP.Response(request.body)
        request.response.request = request
        startwrite(http)
        write(http, request.response.body)
    end

    server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, "127.0.0.1"), port))
    tsk = @async HTTP.listen(handler, "127.0.0.1", port; server=server)
    sleep(3.0)
    @test !istaskdone(tsk)
    r = testget("http://127.0.0.1:$port")
    @test r[1].status == 200
    close(server)
    sleep(0.5)
    @test istaskdone(tsk)

    server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, "127.0.0.1"), port))
    tsk = @async HTTP.listen(handler, "127.0.0.1", port; server=server)

    handler2 = HTTP.Handlers.RequestHandlerFunction(req->HTTP.Response(200, req.body))

    server2 = Sockets.listen(Sockets.InetAddr(parse(IPAddr, "127.0.0.1"), port+100))
    tsk2 = @async HTTP.serve(handler2, "127.0.0.1", port+100; server=server2)
    sleep(0.5)
    @test !istaskdone(tsk)
    @test !istaskdone(tsk2)

    r = testget("http://127.0.0.1:$port")
    @test r[1].status == 200

    r = testget("http://127.0.0.1:$(port+100)")
    @test r[1].status == 200

    rs = testget("http://127.0.0.1:$port/", 20)
    foreach(rs) do r
        @test r.status == 200
    end

    r = HTTP.get("http://127.0.0.1:$port/"; readtimeout=30)
    @test r.status == 200
    @test String(r.body) == ""

    # large headers
    tcp = Sockets.connect(ip"127.0.0.1", port)
    x = "GET / HTTP/1.1\r\n$(repeat("Foo: Bar\r\n", 10000))\r\n";
    @show length(x)
    write(tcp, "GET / HTTP/1.1\r\n$(repeat("Foo: Bar\r\n", 10000))\r\n")
    sleep(0.1)
    @test occursin(r"HTTP/1.1 413 Request Entity Too Large", String(read(tcp)))

    # invalid HTTP
    tcp = Sockets.connect(ip"127.0.0.1", port)
    write(tcp, "GET / HTP/1.1\r\n\r\n")
    sleep(0.1)
    @test occursin(r"HTTP/1.1 400 Bad Request", String(readavailable(tcp)))

    # no URL
    tcp = Sockets.connect(ip"127.0.0.1", port)
    write(tcp, "SOMEMETHOD HTTP/1.1\r\nContent-Length: 0\r\n\r\n")
    sleep(0.1)
    r = String(readavailable(tcp))
    @test occursin(r"HTTP/1.1 400 Bad Request", r)

    # Expect: 100-continue
    tcp = Sockets.connect(ip"127.0.0.1", port)
    write(tcp, "POST / HTTP/1.1\r\nContent-Length: 15\r\nExpect: 100-continue\r\n\r\n")
    sleep(0.1)
    client = String(readavailable(tcp))
    @test client == "HTTP/1.1 100 Continue\r\n\r\n"

    write(tcp, "Body of Request")
    sleep(0.1)
    client = String(readavailable(tcp))

    println("client:")
    println(client)
    @test occursin("HTTP/1.1 200 OK\r\n", client)
    @test occursin("Transfer-Encoding: chunked\r\n", client)
    @test occursin("Body of Request", client)

    hello = (http) -> begin
        request::HTTP.Request = http.message
        request.body = read(http)
        closeread(http)
        request.response::HTTP.Response = HTTP.Response("Hello")
        request.response.request = request
        startwrite(http)
        write(http, request.response.body)
    end
    close(server)
    close(server2)

    # keep-alive vs. close: issue #81
    port += 1
    server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, "127.0.0.1"), port))
    tsk = @async HTTP.listen(hello, "127.0.0.1", port; server=server, verbose=true)
    sleep(0.5)
    @test !istaskdone(tsk)
    tcp = Sockets.connect(ip"127.0.0.1", port)
    write(tcp, "GET / HTTP/1.0\r\n\r\n")
    sleep(0.5)
    client = String(readavailable(tcp))
    @show client
    @test client == "HTTP/1.1 200 OK\r\n\r\nHello"

    # SO_REUSEPORT
    println("Testing server port reuse")
    t1 = @async HTTP.listen(hello, "127.0.0.1", 8089; reuseaddr=true)
    sleep(0.5)
    @test !istaskdone(t1)

    println("Starting second server listening on same port")
    t2 = @async HTTP.listen(hello, "127.0.0.1", 8089; reuseaddr=true)
    sleep(0.5)
    @test Sys.iswindows() ? istaskdone(t2) : !istaskdone(t2)

    println("Starting server on same port without port reuse (throws error)")
    try
        HTTP.listen(hello, "127.0.0.1", 8089)
    catch e
        @test e isa Base.IOError
        @test startswith(e.msg, "listen")
        @test e.code == Base.UV_EADDRINUSE
    end

    # test automatic forwarding of non-sensitive headers
    # this is a server that will "echo" whatever headers were sent to it
    t1 = @async HTTP.listen("127.0.0.1", 8090) do http
        request::HTTP.Request = http.message
        request.body = read(http)
        closeread(http)
        request.response::HTTP.Response = HTTP.Response(200, request.headers)
        request.response.request = request
        startwrite(http)
        write(http, request.response.body)
    end

    sleep(0.5)
    @test !istaskdone(t1)

    # test that an Authorization header is **not** forwarded to a domain different than initial request
    @test_skip !HTTP.hasheader(HTTP.get("http://httpbin.org/redirect-to?url=http://127.0.0.1:8090", ["Authorization"=>"auth"]), "Authorization")

    # test that an Authorization header **is** forwarded to redirect in same domain
    @test_skip HTTP.hasheader(HTTP.get("http://httpbin.org/redirect-to?url=https://httpbin.org/response-headers?Authorization=auth"), "Authorization")

    # 318
    dir = joinpath(dirname(pathof(HTTP)), "../test")
    sslconfig = MbedTLS.SSLConfig(joinpath(dir, "resources/cert.pem"), joinpath(dir, "resources/key.pem"))
    tsk = @async try
        HTTP.listen("127.0.0.1", 8092; sslconfig = sslconfig, verbose=true) do http::HTTP.Stream
            while !eof(http)
                println("body data: ", String(readavailable(http)))
            end
            HTTP.setstatus(http, 200)
            HTTP.startwrite(http)
            write(http, "response body\n")
            write(http, "more response body")
        end
    catch err
        @error err exception = (err, catch_backtrace())
    end
    clientoptions = (;
        require_ssl_verification = false,
    )
    r = HTTP.request("GET", "https://127.0.0.1:8092"; clientoptions...)
    @test_throws HTTP.IOError HTTP.request("GET", "http://127.0.0.1:8092"; clientoptions...)

end # @testset

@testset "HTTP.listen: rate_limit" begin
    io = IOBuffer()
    logger = Base.CoreLogging.SimpleLogger(io)
    server = listen(IPv4(0), 8080)
    @async Base.CoreLogging.with_logger(logger) do
        HTTP.listen("0.0.0.0", 8080; server=server, rate_limit=2//1) do http
            HTTP.setstatus(http, 200)
            HTTP.setheader(http, "Content-Length" => "0")
            HTTP.startwrite(http)
            HTTP.close(http.stream) # close to force a new connection everytime
        end
    end
    # Test requests from the same IP within the limit
    for _ in 1:5
        sleep(0.6) # rate limit allows 2 per second
        @test HTTP.get("http://127.0.0.1:8080").status == 200
    end
    # Test requests from the same IP over the limit
    try
        for _ in 1:5
            sleep(0.2) # rate limit allows 2 per second
            r = HTTP.get("http://127.0.0.1:8080"; retry=false)
            @test r.status == 200
        end
    catch e
        @test e isa HTTP.IOExtras.IOError
    end

    close(server)
    @test occursin("discarding connection from 127.0.0.1 due to rate limiting", String(take!(io)))

    # # Tests to make sure the correct client IP is used (https://github.com/JuliaWeb/HTTP.jl/pull/701)
    # # This test requires a second machine and thus commented out
    #
    # Machine 1
    # @async HTTP.listen("0.0.0.0", 8080; rate_limit=2//1) do http
    #     HTTP.setstatus(http, 200)
    #     HTTP.setheader(http, "Content-Length" => "0")
    #     HTTP.startwrite(http)
    #     HTTP.close(http.stream) # close to force a new connection everytime
    # end
    # while true
    #     sleep(0.6)
    #     print("#") # to show some progress
    #     HTTP.get("http://$(MACHINE_1_IPV4):8080"; retry=false)
    # end
    #
    # # Machine 2
    # while true
    #     sleep(0.6)
    #     print("#") # to show some progress
    #     HTTP.get("http://$(MACHINE_1_IPV4):8080"; retry=false)
    # end

end

@testset "on_shutdown" begin
    @test HTTP.Servers.shutdown(nothing) === nothing

    IOserver = Sockets.listen(Sockets.InetAddr(parse(IPAddr, "127.0.0.1"), 8052))

    # Shutdown adds 1
    TEST_COUNT = Ref(0)
    shutdown_add() = TEST_COUNT[] += 1
    server = HTTP.Servers.Server(nothing, IOserver, "host", "port", shutdown_add)
    close(server)

    # Shutdown adds 1, performed twice
    @test TEST_COUNT[] == 1
    server = HTTP.Servers.Server(nothing, IOserver, "host", "port", [shutdown_add, shutdown_add])
    close(server)
    @test TEST_COUNT[] == 3

    # First shutdown function errors, second adds 1
    shutdown_throw() = throw(ErrorException("Broken"))
    server = HTTP.Servers.Server(nothing, IOserver, "host", "port", [shutdown_throw, shutdown_add])
    @test_logs (:error, r"shutdown function .* failed") close(server)
    @test TEST_COUNT[] == 4
end # @testset

@testset "access logging" begin
    function handler(http)
        if http.message.target == "/internal-error"
            error("internal error")
        end
        if http.message.target == "/close"
            HTTP.setstatus(http, 444)
            close(http.stream)
            return
        end
        HTTP.setstatus(http, 200)
        HTTP.setheader(http, "Content-Type" => "text/plain")
        msg = "hello, world"
        HTTP.setheader(http, "Content-Length" => string(sizeof(msg)))
        HTTP.startwrite(http)
        if http.message.method == "GET"
            HTTP.write(http, msg)
        end
    end
    function with_testserver(f, fmt)
        l = Sockets.listen(ip"0.0.0.0", 1234)
        logger = Test.TestLogger()
        tsk = @async begin
            Base.CoreLogging.with_logger(logger) do
                HTTP.listen(handler, Sockets.localhost, 1234; server=l, access_log=fmt)
            end
        end
        try
            f()
        finally
            close(l)
        end
        return filter!(x -> x.group == :access, logger.logs)
    end

    # Common Log Format
    logs = with_testserver(common_logfmt) do
        HTTP.get("http://localhost:1234")
        HTTP.get("http://localhost:1234/index.html")
        HTTP.get("http://localhost:1234/index.html?a=b")
        HTTP.head("http://localhost:1234")
        HTTP.get("http://localhost:1234/internal-error"; status_exception=false)
        sleep(1) # necessary to properly forget the closed connection from the previous call
        try HTTP.get("http://localhost:1234/close"; retry=false) catch end
        HTTP.get("http://localhost:1234", ["Connection" => "close"])
    end
    @test length(logs) == 7
    @test all(x -> x.group === :access, logs)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET / HTTP/1.1\" 200 12$", logs[1].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /index.html HTTP/1.1\" 200 12$", logs[2].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /index.html\?a=b HTTP/1.1\" 200 12$", logs[3].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"HEAD / HTTP/1.1\" 200 0$", logs[4].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /internal-error HTTP/1.1\" 500 \d+$", logs[5].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /close HTTP/1.1\" 444 0$", logs[6].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET / HTTP/1.1\" 200 12$", logs[7].message)

    # Combined Log Format
    logs = with_testserver(combined_logfmt) do
        HTTP.get("http://localhost:1234", ["Referer" => "julialang.org"])
        HTTP.get("http://localhost:1234/index.html")
        useragent = HTTP.MessageRequest.USER_AGENT[]
        HTTP.setuseragent!(nothing)
        HTTP.get("http://localhost:1234/index.html?a=b")
        HTTP.setuseragent!(useragent)
        HTTP.head("http://localhost:1234")
    end
    @test length(logs) == 4
    @test all(x -> x.group === :access, logs)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET / HTTP/1.1\" 200 12 \"julialang\.org\" \"HTTP\.jl/.*\"$", logs[1].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /index.html HTTP/1.1\" 200 12 \"-\" \"HTTP\.jl/.*\"$", logs[2].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"GET /index.html\?a=b HTTP/1.1\" 200 12 \"-\" \"-\"$", logs[3].message)
    @test occursin(r"^127.0.0.1 - - \[(\d{2})/.*/(\d{4}):\d{2}:\d{2}:\d{2}.*\] \"HEAD / HTTP/1.1\" 200 0 \"-\" \"HTTP\.jl/.*\"$", logs[4].message)

    # Custom log format
    fmt = logfmt"$http_accept $sent_http_content_type $request $request_method $request_uri $remote_addr $remote_port $remote_user $server_protocol $time_iso8601 $time_local $status $body_bytes_sent"
    logs = with_testserver(fmt) do
        HTTP.get("http://localhost:1234", ["Accept" => "application/json"])
        HTTP.get("http://localhost:1234/index.html")
        HTTP.get("http://localhost:1234/index.html?a=b")
        HTTP.head("http://localhost:1234")
    end
    @test length(logs) == 4
    @test all(x -> x.group === :access, logs)
    @test occursin(r"^application/json text/plain GET / HTTP/1\.1 GET / 127\.0\.0\.1 \d+ - HTTP/1\.1 \d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.* \d+/.*/\d{4}:\d{2}:\d{2}:\d{2}.* 200 12$", logs[1].message)
    @test occursin(r"^\*/\* text/plain GET /index\.html HTTP/1\.1 GET /index\.html 127\.0\.0\.1 \d+ - HTTP/1\.1 \d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.* \d+/.*/\d{4}:\d{2}:\d{2}:\d{2}.* 200 12$", logs[2].message)
    @test occursin(r"^\*/\* text/plain GET /index\.html\?a=b HTTP/1\.1 GET /index\.html\?a=b 127\.0\.0\.1 \d+ - HTTP/1\.1 \d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.* \d+/.*/\d{4}:\d{2}:\d{2}:\d{2}.* 200 12$", logs[3].message)
    @test occursin(r"^\*/\* text/plain HEAD / HTTP/1\.1 HEAD / 127\.0\.0\.1 \d+ - HTTP/1\.1 \d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}.* \d+/.*/\d{4}:\d{2}:\d{2}:\d{2}.* 200 0$", logs[4].message)
end

end # module


text:
using StaticNumbers
using Test
using Test: @inferred

@testset "simple tests" begin
    @test static(1) === StaticInteger{1}()
    @test static(1) == 1
    @test sin(static(1)) == sin(1)
    @test static(1+im) === StaticNumber{1+im}()

    @test static(1) + static(1) == 2

    @test promote_type(StaticInteger{1}, StaticInteger{1}) == Int
    @test promote_type(StaticInteger{1}, StaticInteger{2}) == Int

    @test zero(static(1)) === 0
    @test zero(static(1.0)) === 0.0
    @test zero(typeof(static(1))) === 0
    @test zero(typeof(static(1.0))) === 0.0

    @test static(true) == true
    @test static(false) == false

    @test unstatic(1) === 1
    @test unstatic(static(1)) === 1
end

# @testset "ambiguities" begin
#     ambiguities = detect_ambiguities(Base, StaticNumbers)
#     for a in ambiguities
#         println(a[1], "\n", a[2], "\n")
#     end
#     @test length(detect_ambiguities(StaticNumbers)) == 0
#     @test length(ambiguities) <= 5 + 10*(VERSION < v"1.3")
# end

@testset "static math" begin
    if VERSION <= v"1.5"
        # Don't test with pi on Julia versions that don't support one(pi)
        lst = (-1.0, -1, -1//2, -0.5, 0, 0.0, false, true, 2, 3, 1.5, 2.0, 3.1, 3//2, 3.0+im, Inf)
    else
        lst = (-1.0, -1, -1//2, -0.5, 0, 0.0, false, true, 2, 3, 1.5, 2.0, 3.1, pi, 3//2, 3.0+im, Inf)
    end

    for x in lst
        for f in (:round, :ceil, :floor, :sign, :cos, :sin, :log, :exp, :isfinite, :isnan, :abs, :abs2, :iszero, :isone)
            r = try
                @eval $f($x)
            catch
                nothing
            end
            if r !== nothing
                #println(f, (x,), " == ", r)
                @test @eval $f(static($x)) == $r
            end
        end
        for y in lst
            @test static(x) + y === x + y
            @test x + static(y) === x + y
            @test static(x) + static(y) === x + y
            for f in (:-, :*, :/, :^, :widemul, :rem, :mod, :div, :(<<), :(>>), :(==), :(<), :(<=), :(>), :(>=))
                r = try
                    @eval $f($x,$y)
                catch
                    nothing
                end
                if r !== nothing
                    #println(f, (x,y), " ≈ ", r)
                    if isnan(r)
                        @test @eval isnan($f(static($x), $y))
                        @test @eval isnan($f($x, static($y)))
                        @test @eval isnan($f(static($x), static($y)))
                    else
                        @test @eval $f(static($x), $y) ≈ $r
                        @test @eval $f($x, static($y)) ≈ $r
                        @test @eval $f(static($x), static($y)) ≈ $r
                    end
                end
            end
        end
    end
end

@testset "static types" begin
    @test static(1.5) === StaticReal{1.5}()

    @test Static{1}() === static(1)

    @test Float64(static(1.5)) === 1.5

    @test big(static(1)) == big(1)
    @test typeof(big(static(1))) == BigInt

    @test Complex{Float64}(static(2)) === ComplexF64(2)

    @test 2+static(2) == 4
    @test 1+static(1//2) === 3//2
    @test static(1//2)+1 === 3//2

    @test static(static(1)) == static(1)
    @test_throws ErrorException StaticInteger{StaticInteger{1}()}()
    @test_throws InexactError convert(StaticInteger{1}, 2)

    @test StaticInteger{0}() < 1
    @test static(1) isa Integer
    @test +static(1) === static(1)
    @test -static(1) == -1
    @test ~static(1) == -2
    @test isqrt(static(1)) == 1
    @test zero(static(1)) == 0
    @test one(static(2)) == 1
    @test oneunit(static(2)) == 1

    # Skipping these since it might be a bad idea to add methods to Val
    #@test Val(static(false)) === Val(false)
    #@test Val(static(1)) === Val(1)
    #@test Val(static(3.1)) === Val(3.1)
    #@test Val(static(3+im)) === Val(3+im)

    @test static(1):static(3) isa LengthUnitRange{Int64,<:StaticInteger,<:StaticInteger}
    @test static(1:3) === static(1):static(3)
    @test StaticOneTo(3) === static(1):static(3)
    @test static(1):static(3) === @stat 1:3

    @test StaticReal{3.1}(3.1) === static(3.1)
    @test StaticReal{3.25}(3.25f0) === static(3.25)
    @test_throws InexactError StaticReal{3.1}(3.11)

    @test StaticNumber{3+im}(3+im) === static(3+im)
    @test StaticNumber{3+im}(Int32(3)+im) === static(3+im)
    @test_throws InexactError StaticReal{3+im}(3)

    @test static(static(3.14)) === static(3.14)
    @test static(static(3+im)) === static(3+im)

    @test promote_type(StaticInteger{3}, StaticInteger{3}) === Int
    @test promote_type(typeof(static(3.14)), typeof(pi)) === Float64

    for x in (static(0), static(0.0), static(0//1), static(0 + 0im), static(0.0 + 0.0im))
        @test promote(x,x) === (x+0, x+0)
        @test Base.promote_typeof(x,x) === typeof(x+0)
        @test promote_rule(typeof(x), Bool) === typeof(x+0)
    end
    @test promote(static(1), static(1)) === (1,1)
    @test Base.promote_typeof(static(1), static(1)) === Int

    @test static(1):static(2):5 isa LengthStepRange{Int64,StaticInteger{-1},StaticInteger{2},Int64}

    @test Rational{Int}(static(3//2)) === 3//2

    @test Base.promote_rule(StaticInteger{3}, StaticReal{3.14}) === Float64

    @test widemul(static(1), true) === widemul(1, true)
    @test widemul(static(1.0), false) === widemul(1.0, false)

    @test isnan(static(NaN)) === true
    @test isnan(static(0.0)) === false
end

@testset "show" begin
    @test sprint(show, static(1)) == "static(1)"
    @test sprint(show, static(1//2)) == "static(1//2)"
    @test sprint(show, static(0.5)) == "static(0.5)"
    @test sprint(show, static(1:2:3)) == "static(1:2:3)"
    @test sprint(show, static(2:2:3)) == "static(2:2:2)"
end

@testset "broadcasting" begin
    @test static(1:2) .^ 2 == [1, 4]
    @test static(1:2) .* static(3:4)' == [3 4; 6 8]
    @test static(1:2) .* static(3:4) == [3, 8]
end

@testset "trystatic" begin
    @test trystatic(0, static(1)) === 0
    @test trystatic(1, static(1)) === static(1)
    @test trystatic(0, static(0), static(1)) === static(0)
    @test trystatic(1, static(0), static(1)) === static(1)
    @test trystatic(2, static(0), static(1)) === 2

    @test trystatic(0, 1) === 0
    @test trystatic(static(0), 1) === static(0)
    @test trystatic(1, 1) === static(1)
    @test trystatic(0, 0, 1) === static(0)
    @test trystatic(1, 0, 1) === static(1)
    @test trystatic(2, 0, 1) === 2
    @test trystatic(2, (0, 1)) === 2
    @test trystatic(2, (0, 1, 2)) === static(2)
    @test 2 ⩢ (0, 1) === 2
    @test 2 ⩢ (0, 1, 2) === static(2)
    @test 2 ⩢ 0 ⩢ 1 ⩢ 2 === static(2)
    @test 1 ⩢ 0 ⩢ 1 ⩢ 2 === static(1)
    @test 2 ⩢ 0 ⩢ 1 === 2

    @test trystatic(0, NaN) === 0
    @test trystatic(0, static(NaN)) === 0
    @test trystatic(0.0, static(NaN)) === 0.0
    @test trystatic(NaN, NaN) === static(NaN)
    @test trystatic(NaN, static(NaN)) === static(NaN)
    @test trystatic(0.0/0.0, NaN) === static(NaN)
    @test 0.0 ⩢ NaN === 0.0
    @test NaN ⩢ NaN === static(NaN)

    @test trystatic(2, 1:3) == static(2)
    @test trystatic(4, 1:3) == 4
end

@testset "@generate_static_methods macro" begin
    #println(macroexpand(StaticNumbers, :(@generate_static_methods (0, 1) (Base.Math.sinpi, Base.Math.cospi) (Base.:+, Base.:-) )))
    @generate_static_methods (0, 1) (Base.Math.sinpi, Base.Math.cospi) (Base.:+, Base.:-)
    @test sinpi(StaticInteger{1}()) === StaticInteger{0}()
    @test StaticInteger{1}() - StaticInteger{1}() === StaticInteger{0}()

    # Test that sqrt(-1) doesn't cause problems with @generate_static_methods
    @generate_static_methods (-1, 0, 1) (Base.sqrt,) (Base.rem,)
    @test sqrt(static(1)) === static(1)
    @test rem(static(1), static(1)) === static(0)
end

@testset "various" begin
    @test ntuple(identity, static(5)) === static.(ntuple(identity, Val(5)))
    @inferred ntuple(identity, static(1))
    @inferred ntuple(identity, static(5))
    @inferred ntuple(identity, static(25))
    @inferred ntuple(static, static(1))
    @inferred ntuple(static, static(5))
    @inferred ntuple(static, static(25))
    @inferred ntuple(x->x^2, static(1))
    @inferred ntuple(x->x^2, static(5))
    @inferred ntuple(x->x^2, static(25))

    # Test StaticRanges

    @test range(1, length=static(3)) === staticlength(1:3)
    @test range(1, length=static(3)) == 1:3
    @test range(3, step=static(2), length=static(3)) == 3:2:7
    @test range(3, step=2, length=static(3)) == 3:2:7
    @test range(1, length=3) == 1:3
    @test length(range(1, length=static(3))) === static(3)
    @test length(staticlength(1:3)) === static(3)

    @test range(2, length=static(3)) == 2:4
    @test typeof(range(1, length=static(3))) == typeof(range(2, length=static(3)))

    @test staticlength(1:3) isa LengthRange
    @test staticlength(Base.OneTo(3)) isa LengthRange
    @test length(staticlength(Base.OneTo(3))) === static(3)

    r = LengthStepRange(1, 2, static(3))

    @test r isa LengthRange
    @test r isa LengthRange{Int, Int, Int, <:Static}
    @test r isa LengthRange{Int, Int, Int, <:Static{3}}
    @test r isa LengthRange{Int, Int, Int, StaticInteger{3}}
    @test r isa StaticNumbers.LengthStepRange{Int, Int, Int, StaticInteger{3}}
    @test all(r .== 3:2:7)
    @test r[2] == 5
    @test all(collect(r) .== [3, 5, 7])
    @test all(2*r .== 2*(3:2:7))
    @test all(r*5 .== (3:2:7)*5)
    @test all(2 .* r .== 2 .* (3:2:7))
    @test all(r .* 5 .== (3:2:7) .* 5)
    @test all(7 .+ r .== 7 .+ (3:2:7))
    @test all(r .+ 7 .== (3:2:7) .+ 7)
    @test all(7 .- r .== 7 .- (3:2:7))
    @test all(r .- 7 .== (3:2:7) .- 7)
    @test all(-r .== -(3:2:7))
    @test +r === r
    @test .+r === r
    @test .-(.-r) === r
    @test -(-r) === r
    @test typeof(7 .+ r) == typeof(r)
    @test typeof(r .+ 7) == typeof(r)
    @test typeof(7 .- r) == typeof(r)
    @test typeof(r .- 7) == typeof(r)
    @test typeof(7 .* r) == typeof(r)
    @test typeof(r .* 7) == typeof(r)
    @test typeof(7 .+ 5 .* r) == typeof(r)
    @test typeof(r .* 5 .+ 7) == typeof(r)
    @test all( (7:3:100)[r] .== (7:3:100)[3:2:7] )
    @test LengthRange(r) == r

    @test !(static(false)) === true
    @test !(static(true)) === false
    @test true isa StaticOrBool
    @test static(true) isa StaticOrBool

    @test 1 isa StaticOrInt
    @test static(1) isa StaticOrInt

    @test @stat(oftype(static(1), static(2))) === static(2)
    @test @stat(oftype(1, static(2))) === static(2)
    @test @stat(oftype(1.0, 2)) === 2.0
    @test @stat(oftype(static(1.0), 2)) === 2.0
    @test @stat(oftype(1.0, static(2))) === static(2.0)

    ur = LengthUnitRange(2, static(3))
    @test ur isa LengthRange
    @test ur isa LengthRange{Int, Int, <:Static, <:Static}
    @test ur isa LengthRange{Int, Int, <:Static{1}, <:Static{3}}
    @test ur isa LengthRange{Int, Int, StaticInteger{1}, StaticInteger{3}}
    @test ur isa StaticNumbers.LengthUnitRange{Int, Int, StaticInteger{3}}
    @test all(ur .== 3:5)
    @test ur[2] == 4
    @test all(collect(ur) .== [3, 4, 5])
    @test all(2*ur .== 2*(3:5))
    @test all(ur*5 .== (3:5)*5)
    @test all(2 .* ur .== 2 .* (3:5))
    @test all(ur .* 5 .== (3:5) .* 5)
    @test all(7 .+ ur .== 7 .+ (3:5))
    @test all(ur .+ 7 .== (3:5) .+ 7)
    @test all(-ur .== -(3:5))
    @test +ur === ur
    @test .+ur === ur
    @test .-(.-ur) === ur
    @test -(-ur) === ur
    @test typeof(7 .+ ur) == typeof(ur)
    @test typeof(ur .+ 7) == typeof(ur)
    @test length(7 .- ur) == static(3)
    @test typeof(ur .- 7) == typeof(ur)
    @test length(7 .* ur) === static(3)
    @test length(ur .* 7) === static(3)
    @test length(7 .+ 5 .* ur) === static(3)
    @test length(ur .* 5 .+ 7) === static(3)
    @test all( (7:3:100)[ur] .== (7:3:100)[3:5] )

    # Test types
    @test LengthRange(1:3) isa LengthUnitRange
    @test LengthRange(1:2:4) isa LengthStepRange

    # Test that type inferrence is working
    @inferred 2*r
    @inferred 2*ur

    f3(x) = trystatic(x, 3, 4)
    # @show Base.return_types(f3, (Int,))
    g3() = f3(4)
    @inferred g3()

    f4(x) = trystatic(x, static(3), static(4))
    # @show Base.return_types(f4, (Int,))
    g4() = f4(4)
    @inferred g4()

     # f5(x) = trystatic(mod(x,4), LengthStepRange(static(-1),static(1),static(4)))
     # @show Base.return_types(f5, (Int,))
     # g5() = f5(2)
     # @show Base.return_types(g5, ())

    # Test array handling with static ranges
    A = rand(16,16)
    B = rand(static(16),static(16))
    C = A[staticlength(5:8),StaticOneTo(4)]
    @test all(C .== A[staticlength(5:8),StaticOneTo(4)])
    A[StaticOneTo(4),StaticOneTo(4)] = C
    @test all(A[StaticOneTo(4),StaticOneTo(4)] .== C)
    @test all(staticlength(3:4).^2 == (3:4).^2)

    @test Unsigned(static(2)) === Unsigned(2)
end

# Test with a new numeric type
struct MyType1 <: Real
    x::Float64
end

@testset "various II" begin
    @test MyType1(static(3)) === MyType1(3)
    @test MyType1(static(3.0)) === MyType1(3.0)

    @test maybe_static(+, 2, 2) === 4
    @test maybe_static(+, static(2), 2) === 4
    @test maybe_static(+, 2, static(2)) === 4
    @test maybe_static(+, static(2), static(2)) === static(4)
    @inferred maybe_static(+, 2, 2)
    @inferred maybe_static(+, static(2), 2)
    @inferred maybe_static(+, 2, static(2))
    @inferred maybe_static(+, static(2), static(2))

    # Test @tostatic macro

    x = 0
    @test StaticNumbers.@tostatic(x, 0, 3) === static(0)
    @test StaticNumbers.@tostatic(x, -1, 0) === static(0)
    @test StaticNumbers.@tostatic(x, 0, 0) === static(0)

    x = 3
    @test StaticNumbers.@tostatic(x, 0, 3) === static(3)
    @test StaticNumbers.@tostatic(x, 2, 7) === static(3)
    @test StaticNumbers.@tostatic(x, 3, 3) === static(3)
    @test StaticNumbers.@tostatic(x, 3, 4) === static(3)

    StaticNumbers.@tostatic x 2 5 begin
        @test x === static(3)
    end

    x = static(3)
    @test StaticNumbers.@tostatic(x, 0, 0) === static(3)

    @test eachindex(static(2:3)) === StaticOneTo(2)

    @test static(3) & static(3) === static(3)
    @test static(3) | static(3) === static(3)
    @test static(3) ⊻ static(3) === static(0)
    @test fma(static(3), static(3), static(3)) === 3*3+3
    @test fma(static(3.0), static(3.0), static(3.0)) === 3.0*3.0+3.0
end

@testset "tuple indexing" begin
    for t in ((1,2,3,4), (1, 2.0, 3//1, 4.0f0))
        @test t[static(2)] === t[2]
        @inferred t[static(2)]
        @test t[static(2):static(3)] === t[2:3]
        @inferred t[static(2):static(3)]
        @test t[static(2):static(2):static(4)] === t[2:2:4]
        @inferred t[static(2):static(2):static(4)]
        @test t[static(1):static(2):static(3)] === t[1:2:3]
        @inferred t[static(1):static(2):static(2)]
    end
end

@testset "more tests" begin
    @test StaticOneTo(Base.OneTo(3)) === StaticOneTo(3)
    @test StaticOneTo(1:3) === StaticOneTo(3)
    @test LengthUnitRange(2:4) isa LengthUnitRange
    @test LengthUnitRange(2:4) == 2:4
end

# Expand `@stat` in a clean namespace to avoid assuming anything in
# the expanding namespace.
module CleanNamespace
using StaticNumbers: @stat
add_literal_one(x) = @stat x + 1
add(x, y) = @stat x + y
get_one_to(x, i) = @stat x[1:i]
end

@testset "@stat macro" begin
    @test @stat(2+2) === static(4)
    x = 2
    x2 = @stat x + 2
    @test x2 === 4
    y = static(2)
    y2 = @stat y + 2
    @test y2 === static(4)

    T = (1,2,3,4)
    A = [1,2,3,4]
    @test @stat(length(A)) === 4
    @test @stat(length(T)) === static(4)
    @test @stat(firstindex(T)) === static(1)
    @test @stat(lastindex(A)) === 4
    @test @stat(lastindex(T)) === static(4)

    r = @stat 2:5
    @test r[1] === first(r)
    @test first(r) === 2
    @test static(2) === @stat first(r)
    @test r[end] === last(r)
    @test last(r) === 5
    @test static(5) === @stat last(r)
    @test static(3) === @stat r[2]
    i = 2
    @test 3 === @stat r[i]
    #@test static(3:4) === @stat (1:5)[3:4] # TODO

    f(t) = @stat t[2:end-1]
    @test f(T) === (2,3)
    @inferred f((1,2,3,4))

    @stat g(t,k) = t[k .+ (0:1)]
    @test g(T, 2) === (2,3)
    @inferred g((1,2,3,4), 2)

    @test @stat(first(2:3)) === static(2)
    @test @stat(last(2:3)) === static(3)
    @test static(2:3)[2] === 3
    @test static(2:3)[static(2)] === 3
    @test @stat((2:3)[2]) === static(3)

    # At the moment, we don't want @stat to create static unsigned numbers. (Subject to change.)
    @test UInt(1) === @stat UInt(1)
    @test UInt(1) === @stat unsigned(1)

    T = (1,2,3,4)
    f1(T) = @stat T[2:end-1]
    @test (2,3) === f1(T)
    @inferred f1(T)
    f2(T) = @stat T[range(2, length=2)]
    @test (2,3) === f2(T)
    @inferred f2(T)
    @inferred T[range(2, length=static(2))]

    @test CleanNamespace.add_literal_one(0) === 1
    @test CleanNamespace.add_literal_one(static(0)) === static(1)
    @test CleanNamespace.add(static(0), 1) === 1
    @test CleanNamespace.add(static(0), static(1)) === static(1)
end

tuple_test(n) = Tuple(i for i in static(1):n)
tuple_test2(n) = Tuple(tuple_test(n) for i in static(1):n)
tuple_test3(n) = Tuple(Tuple(i+j for i in static(1):n) for j in static(1:n))

@testset "Tuple generators" begin
    @test Tuple(10i+j for i=static(1:3), j=static(1:2)) === Tuple(10i+j for i=1:3, j=1:2)
    @test Tuple(i+j for i=static(2:3), j=static(10:10:20)) === (12, 13, 22, 23)
    @test tuple_test(static(2)) === static.((1, 2))
    @test tuple_test2(static(2)) === (static.((1, 2)), static.((1, 2)))
    @test tuple_test3(static(2)) === ((2, 3), (3, 4))

    if VERSION >= v"1.3"
        # Some of the inference tests are too hard for Julia 1.0
        @inferred Tuple(i+j for i=static(2:3), j=static(10:10:20))
        @inferred tuple_test(static(2))
        @inferred tuple_test2(static(2))
        @inferred tuple_test3(static(2))
    end
end

@testset "examples for doc" begin
    i = 2
    s = static(2)
    @test 4 === s + s
    @test static(4) === @stat s + s
    @test static(4) === @stat s + 2
    @test 4 === @stat s + i
    @test (1, 4, 9, 16) === Tuple(i^2 for i in static(1):static(4))
    @inferred Tuple(i^2 for i in static(1):static(4))
end

# NOTE: There's still a lot of work to do on unsigned
@testset "unsigned" begin
    @test static(UInt(2)) > -1
    @test UInt(2) > static(-1)
    @test static(UInt(2)) > static(-1)
    @test static(UInt(2)) >= -1
    @test UInt(2) >= static(-1)
    @test static(UInt(2)) >= static(-1)
    @test static(1) << UInt(2) === 4
    @test static(4) >> UInt(2) === 1
end

@testset "setindex(::Tuple, _, ::StaticInteger)" begin
    @test @inferred(Base.setindex((1, 2, 3), nothing, static(1))) == (nothing, 2, 3)
    @test @inferred(Base.setindex((1, 2, 3), nothing, static(2))) == (1, nothing, 3)
    @test @inferred(Base.setindex((1, 2, 3), nothing, static(3))) == (1, 2, nothing)
    @test_throws BoundsError Base.setindex((1, 2, 3), nothing, static(0))
    @test_throws BoundsError Base.setindex((1, 2, 3), nothing, static(4))
end

@testset "StaticNumbers.map" begin
    long_tuple = (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35)
    @test StaticNumbers.map(sin, (1,2,3)) === Base.map(sin, (1,2,3))
    @test StaticNumbers.map(sin, long_tuple) === Base.map(sin, long_tuple)
    @test StaticNumbers.map(+, (1,2,3), (1,2,3)) === Base.map(+, (1,2,3), (1,2,3))
    @test StaticNumbers.map(+, long_tuple, long_tuple) === Base.map(+, long_tuple, long_tuple)
    @test StaticNumbers.map(+, (1,2,3), (1,2,3), (1,2,3)) === Base.map(+, (1,2,3), (1,2,3), (1,2,3))
    @test StaticNumbers.map(+, long_tuple, long_tuple, long_tuple) === Base.map(+, long_tuple, long_tuple, long_tuple)
    if VERSION >= v"1.3"
        @inferred StaticNumbers.map(sin, (1,2,3))
        @inferred StaticNumbers.map(sin, long_tuple)
        @inferred StaticNumbers.map(+, (1,2,3), (1,2,3))
        @inferred StaticNumbers.map(+, long_tuple, long_tuple)
        @inferred StaticNumbers.map(+, (1,2,3), (1,2,3), (1,2,3))
        @inferred StaticNumbers.map(+, long_tuple, long_tuple, long_tuple)
    end
end

@testset "Tuple from view" begin
    A = rand(5,5)
    v = view(A, range(2, length=static(3)), range(3, length=static(2)))
    S = @inferred Tuple(v)
    @test S === Tuple(A[2:4,3:4])
end

include("StaticArrays_test.jl")

include("SIMD_test.jl")

# @testset "ambiguities (extended)" begin
#     ambiguities = detect_ambiguities(Base, StaticNumbers, StaticArrays, SIMD)
#     for a in ambiguities
#         println(a[1], "\n", a[2], "\n")
#     end
#     @test length(ambiguities) <= 6 + 10*(VERSION < v"1.3")
# end

# Run all tests in optional packages, together with StaticNumbers,
# to make sure we didn't break anything.
# if false
#     # This currenlty breaks because both SIMD and StaticArrays export "setindex"
#     for m in (StaticArrays, SIMD)
#         println("Running tests from ", m)
#         include(joinpath(dirname(dirname(pathof(m))), "test", "runtests.jl"))
#     end
# end


text:
function part1(lines)
    departure = parse(Int, first(lines))
    ids = split(last(lines), ',')
    inservice = parse.(Int, filter(x -> x != "x", ids))

    wait = inservice .* ceil.(Int, departure ./ inservice) .- departure
    _, i = findmin(wait)

    wait[i] * inservice[i]
end

function part2(lines)
    fields = split(last(lines), ",")

    ranges = map(filter(!=("x") ∘ last, collect(enumerate(fields)))) do elem
        δ = first(elem) - 1
        c = parse(Int, last(elem))

        # BigInt is necessary to prevent integer overflow inside the
        # `intersect` to come
        c .* (one(BigInt):typemax(Int)) .- δ
    end

    first(intersect(ranges...))
end

function main()
    lines = readlines("input.txt")
    @info "Solution" part1(lines) part2(lines)
end

main()

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
using GLVisualize, GLAbstraction, Colors, Reactive, GeometryTypes

if !isdefined(:runtests)
	window = glscreen()
	timesignal = loop(linspace(0f0,1f0,360))
end

# generate some pretty data
function xy_data(x,y,i, N)
	x = ((x/N)-0.5f0)*i
	y = ((y/N)-0.5f0)*i
	r = sqrt(x*x + y*y)
	Float32(sin(r)/r)
end

surf(i, N) = Float32[xy_data(x, y, i, N) for x=1:N, y=1:N]

t = map(t->(t*30f0)+20f0, timesignal)

bb = Signal(AABB{Float32}(Vec3f0(0), Vec3f0(1)))

_view(visualize(const_lift(surf, t, 400), :surface, boundingbox=bb))

if !isdefined(:runtests)
	renderloop(window)
end


text:
__precompile__()

module MultiScaleArrays

import Base: length, push!, deleteat!, getindex, setindex!, eachindex, ndims,
       size, print_matrix, similar, hcat, vcat,
       ==, *, +, /, -, show, vec, reshape

import RecursiveArrayTools: recursivecopy,recursivecopy!

import RecursiveArrayTools: chain

"""
# Defining A MultiScaleModel: The Interface

The required interface is as follows. Leaf types must extend AbstractMultiScaleArrayLeaf, the
highest level of the model or the head extends MultiScaleModelHead, and all
intermediate types extend AbstractMultiScaleArray. The leaf has an array `values::Vector{B}`.
Each type above then contains three fields:

- `nodes::Vector{T}`
- `values::Vector{B}`
- `end_idxs::Vector{Int}`

Note that the ordering of the fields matters.
`B` is the `BottomType`, which has to be the same as the eltype for the array
in the leaf types. `T` is another `AbstractMultiScaleArray`. Thus at each level,
an` AbstractMultiScaleArray` contains some information of its own (`values`), the
next level down in the heirarchy (`nodes`), and caching for indices (`end_idxs`).
You can add and use extra fields as you please, and even make the types immutable.

## The MultiScaleModel API

The resulting type acts as an array. A leaf type `l` acts exactly as an array
with `l[i] == l.values[i]`. Higher nodes also act as a linear array. If `ln` is level
`n` in the heirarchy, then `ln.nodes` is the vector of level `n-1` objects, and `ln.values`
are its "intrinsic values". There is an indexing scheme on `ln`, where:

- `ln[i,j,k]` gets the `k`th `n-3` object in the `j`th `n-2` object in the `i`th level `n-1`
  object. Of course, this recurses for the whole hierarchy.
- `ln[i]` provides a linear index through all `.nodes` and `.values` values in every lower
  level and `ln.values` itself.

Thus `typeof(ln) <: AbstractVector{B}` where `B` is the eltype of its leaves and
all `.values`'s.

In addition, iterators are provided to make it easy to iterate through levels.
For `h` being the head node, `level_iter(h,n)` iterates through all level objects
`n` levels down from the top, while `level_iter_idx(h,n)` is an enumeration
`(node,y,z)` where `node` are the `n`th from the head objects, with `h[y:z]` being
the values it holds in the linear indexing.

## Indexing and Iteration

The head node then acts as the king. It is designed to have functionality which
mimics a vector in order for usage in DifferentialEquations or Optim. So for example

```julia
embryo[12]
```

returns the "12th protein", counting by Embryo > Tissue > Population > Cell in order
of the vectors. The linear indexing exists for every `AbstractMultiScaleArray`.
These types act as full linear vectors, so standard operations do the sensical
operations:

```julia
embryo[10] = 4.0 # changes protein concentration 10
embryo[2,3,1] # Gives the 1st cell in the 3rd population of the second tissue
embryo[:] # generates a vector of all of the protein concentrations
eachindex(embryo) # generates an iterator for the indices
```

Continuous models can thus be written at the protein level and will work seamlessly
with DifferentialEquations or Optim which will treat it like a vector of protein concentrations.
Using the iterators, note that we can get each cell population by looping through
2 levels below the top, so

```julia
for cell in level_iter(embryo,3)
  # Do something with the cells!
end
```

or the multiple level iter, which is the one generally used in
DifferentialEquations.jl functions:

```julia
for (cell, dcell) in LevelIter(3,embryo, dembryo)
    # If these are similar structures, `cell` and `dcell` are the similar parts
    cell_ode(dcell,cell,p,t)
end
```

`LevelIterIdx` can give the indices along with iteration:

```julia
for (cell, y, z) in LevelIterIdx(embryo, 3)
    # cell = embryo[y:z]
end
```

However, the interesting behavior comes from event handling. Since `embryo` will be the
"vector" for the differential equation or optimization problem, it will be the value
passed to the event handling. MultiScaleArrays includes behavior for changing the
structure. For example:

```julia
tissue3 = construct(Tissue, deepcopy([population, population2]))
add_node!(embryo, tissue3) # Adds a new tissue to the embryo
remove_node!(embryo, 2, 1) # Removes population 1 from tissue 2 of the embryo
```

Combined with event handling, this allows for dynamic structures to be derived from
low level behaviors.

## Heterogeneous Nodes via Tuples

Note that tuples can be used as well. This allows for type-stable broadcasting with
heterogeneous nodes. This could be useful for mixing types
inside of the nodes. For example:

```julia
struct PlantSettings{T} x::T end
struct OrganParams{T} y::T end

struct Organ{B<:Number,P} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
    name::Symbol
    params::P
end

struct Plant{B,S,N<:Tuple{Vararg{<:Organ{<:Number}}}} <: AbstractMultiScaleArray{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
    settings::S
end

struct Community{B,N<:Tuple{Vararg{<:Plant{<:Number}}}} <: AbstractMultiScaleArray{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
end

mutable struct Scenario{B,N<:Tuple{Vararg{<:Community{<:Number}}}} <: AbstractMultiScaleArrayHead{B}
    nodes::N
    values::Vector{B}
    end_idxs::Vector{Int}
end

organ1 = Organ([1.1,2.1,3.1], :Shoot, OrganParams(:grows_up))
organ2 = Organ([4.1,5.1,6.1], :Root, OrganParams("grows down"))
organ3 = Organ([1.2,2.2,3.2], :Shoot, OrganParams(true))
organ4 = Organ([4.2,5.2,6.2], :Root, OrganParams(1//3))
plant1 = construct(Plant, (deepcopy(organ1), deepcopy(organ2)), Float64[], PlantSettings(1))
plant2 = construct(Plant, (deepcopy(organ3), deepcopy(organ4)), Float64[], PlantSettings(1.0))
community = construct(Community, (deepcopy(plant1), deepcopy(plant2), ))
scenario = construct(Scenario, (deepcopy(community),))
```

(of course at the cost of mutability).
"""
abstract type AbstractMultiScaleArray{B}     <: AbstractVector{B} end
abstract type AbstractMultiScaleArrayLeaf{B} <: AbstractMultiScaleArray{B} end
abstract type AbstractMultiScaleArrayHead{B} <: AbstractMultiScaleArray{B} end

using DiffEqBase, Statistics, LinearAlgebra, FiniteDiff
import OrdinaryDiffEq, StochasticDiffEq, ForwardDiff, SparseDiffTools

Base.show(io::IO, x::AbstractMultiScaleArray) = invoke(show, Tuple{IO, Any}, io, x)
Base.show(io::IO, ::MIME"text/plain", x::AbstractMultiScaleArray) = show(io, x)

include("shape_construction.jl")
include("addition_deletion.jl")
include("indexing.jl")
include("math.jl")
include("level_iterations.jl")
include("diffeq.jl")
include("print_human_readable.jl")

using TreeViews
TreeViews.hastreeview(x::AbstractMultiScaleArray) = true
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::AbstractMultiScaleArray) = x

# Types
export AbstractMultiScaleArray, AbstractMultiScaleArrayLeaf,
       AbstractMultiScaleArrayHead

# Constructors
export construct, recursivecopy!

# Addition Deletion
export add_node!, remove_node!

# Indexing
export num_nodes, getindices

# Misc
export LevelIterIdx, LevelIter, level_iter

# print_human_readable
export print_human_readable

end # module


text:
using ModelGraphSolvers
using Test

@testset "ModelGraphSolvers.jl" begin
    # Write your own tests here.
end