text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

# Dictionary of medium data

const excludeEnthalpyOfFormation = true                         # If true, enthalpy of formation Hf is not included in specific enthalpy hc
const referenceEnthalpy          = ReferenceEnthalpy_ZeroAt0K   # Choice of reference enthalpy
const h_offset                   = 0.0                          # User defined offset for reference enthalpy, if ReferenceEnthalpy = ReferenceEnthalpy_UserDefined

### Data structures

"""
    data = SingleGasNasaData(;name=Missing, MM=nothing, Hf=nothing, H0=nothing, Tlimit=nothing,
                              alow=Missing, blow=Missing, ahigh=Missing,
                              bhigh=Missing, R=nothing)

Generate a `SingleGasNasaData` object containing the data
for an ideal Gas based on the NASA Glenn coefficients.
"""
mutable struct SingleGasNasaData
    name::AbstractString       # Name of ideal gas
    MM::Float64                # SI.MolarMass "Molar mass"
    Hf::Float64                # SI.SpecificEnthalpy "Enthalpy of formation at 298.15K"
    H0::Float64                # SI.SpecificEnthalpy "H0(298.15K) - H0(0K)"
    Tlimit::Float64            # SI.Temperature Tlimit "Temperature limit between low and high data sets"
    alow::SVector{7,Float64}   # "Low temperature coefficients a"
    blow::SVector{2,Float64}   # "Low temperature constants b"
    ahigh::SVector{7,Float64}  # "High temperature coefficients a"
    bhigh::SVector{2,Float64}  # "High temperature constants b"
    R::Float64                 # SI.SpecificHeatCapacity R "Gas constant"

    SingleGasNasaData(;name=nothing,
                       MM=nothing,
                       Hf=nothing,
                       H0=nothing,
                       Tlimit=nothing,
                       alow=nothing,
                       blow=nothing,
                       ahigh=nothing,
                       bhigh=nothing,
                       R=nothing) =
                       new(name, MM,Hf,H0,Tlimit,alow,blow,ahigh,bhigh,R)
end

"""
    medium = SingleGasNasa(; mediumName     = Missing,
                             reference_p    = 101325,
                             reference_T    = 298.15,
                             p_default      = 101325,
                             T_default      = 293.15,
                             fluidConstants = nothing,
                             fluidLimits    = FluidLimits(TMIN=200.0, TMAX=6000.0),
                             data           = nothing)

Generate a `SingleGasNasa <: PureSubstance` medium object.
"""
struct SingleGasNasa <: PureSubstance
    infos::FluidInfos
    fluidConstants::SVector{1,IdealGasFluidConstants}
    fluidLimits::FluidLimits
    data::SingleGasNasaData

    function SingleGasNasa(; mediumName=nothing,
                             reference_p=101325,
                             reference_T=298.15,
                             p_default=101325,
                             T_default=293.15,
                             fluidConstants=nothing,
                             fluidLimits=FluidLimits(TMIN=200.0, TMAX=6000.0),
                             data=nothing)

        infos = FluidInfos(mediumName           = mediumName,
                           substanceNames       = [mediumName],
                           extraPropertiesNames = fill("",0),
                           ThermoStates         = IndependentVariables_pT,
                           baseProperties       = :BaseProperties_SingleGasNasa,
                           singleState          = false,
                           reducedX             = true,
                           fixedX               = false,
                           reference_p          = reference_p,
                           reference_T          = reference_T,
                           reference_X          = fill(1.0,1),
                           p_default            = p_default,
                           T_default            = T_default,
                           h_default            = h_T(data, T_default))

        fluidConstants.molarMass = data.MM
        fluidLimits.HMIN         = h_T(data, fluidLimits.TMIN)
        fluidLimits.HMAX         = h_T(data, fluidLimits.TMAX)

        new(infos, SVector{1,IdealGasFluidConstants}(fluidConstants), fluidLimits, data)
    end
end

### Functions specific for SingleGasNasa

"cp = cp_T(data::SingleGasNasaData, T) - Compute specific heat capacity at constant pressure from temperature and gas data"
cp_T(data::SingleGasNasaData, T::Float64)::Float64 = if T<data_Tlimit; data_R*
                                                (1/(T*T)*
                                                (data.alow[1]+T*
                                                (data.alow[2]+T*
                                                (1.0*data.alow[3]+T*
                                                (data.alow[4]+T*
                                                (data.alow[5]+T*(data.alow[6]+data.alow[7]*T))))))) else data.R*
                                                (1/(T*T)*
                                                (data.ahigh[1]+T*
                                                (data.ahigh[2]+T*
                                                (1.0*data.ahigh[3]+T*
                                                (data.ahigh[4]+T*
                                                (data.ahigh[5]+T*(data.ahigh[6]+data.ahigh[7]*T))))))) end

"""
    h = h_T(data, T, exclEnthForm=true, refChoice=ReferenceEnthalpy_ZeroAt0K, h_off=0.0)

Return specific enthalpy from temperature and gas data.

# Arguments

- `data::SingleGasNasaData`: Data of the SingleGasNasa medium.
- `T::Float64`: Temperature in [K].
- `exclEnthForm::Bool`: If true, enthalpy of formation Hf is not included in specific enthalpy h.
- `refChoice::ReferenceEnthalpy`: Choice of reference enthalpy.
- `h_off::Float64`: User defined offset for reference enthalpy, if referenceEnthalpy = ReferenceEnthalpy_UserDefined
"""
h_T(data::SingleGasNasaData, T::Float64, exclEnthForm::Bool=excludeEnthalpyOfFormation,
    refChoice::ReferenceEnthalpy=referenceEnthalpy, h_off::Float64=h_offset)::Float64 =
                 (if T<data.Tlimit; data.R*
                 (
                 (-data.alow[1]+T*
                 (data.blow[1]+data.alow[2]*log(T)+T*
                 (1.0*data.alow[3]+T*
                 (0.5*data.alow[4]+T*
                 (1/3*data.alow[5]+T*(0.25*data.alow[6]+0.2*data.alow[7]*T))))))/T) else data.R*
                 (
                 (-data.ahigh[1]+T*
                 (data.bhigh[1]+data.ahigh[2]*log(T)+T*
                 (1.0*data.ahigh[3]+T*
                 (0.5*data.ahigh[4]+T*
                 (1/3*data.ahigh[5]+T*(0.25*data.ahigh[6]+0.2*data.ahigh[7]*T))))))/T) end)+(if exclEnthForm; -data.Hf else 0.0 end)+
                 (if refChoice==ReferenceEnthalpy_ZeroAt0K   ; data.H0 else 0.0 end)+
                 (if refChoice==ReferenceEnthalpy_UserDefined; h_off   else 0.0 end)

"Compute specific entropy from temperature and gas data"
s0_T(data::SingleGasNasaData,   # Ideal gas data
     T::Float64  # Temperature
    )::Float64 = (if T<data.Tlimit; data.R*
                 (data.blow[2]-0.5*data.alow[1]/(T*T)-data.alow[2]/T+data.alow[3]*log(T)+T*
                 (data.alow[4]+T*
                 (0.5*data.alow[5]+T*(1/3*data.alow[6]+0.25*data.alow[7]*T)))) else data.R*
                 (data.bhigh[2]-0.5*data.ahigh[1]/(T*T)-data.ahigh[2]/T+data.ahigh[3]*log(T)+T*
                 (data.ahigh[4]+T*
                 (0.5*data.ahigh[5]+T*(1/3*data.ahigh[6]+0.25*data.ahigh[7]*T)))) end)

### Set states
setState_pTX(m::SingleGasNasa,p,T,X) = ThermodynamicState_pT(p,T)
#setState_phX(m::SingleGasNasa,p,h,X) = ???     solve nonlinear equation with Brent algorithm (needs to be transferred from Modelica to Modia and stord in ModiaMath)
#setState_psX(m::SingleGasNasa,p,s,X) = ???     solve nonlinear equation with Brent algorithm
setState_dTX(m::SingleGasNasa,d,T,X) = ThermodynamicState_pT(d*m.data.R*T, T)

pressure(              m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = state.p
temperature(           m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = state.T
density(               m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = state.p/(m.data.R*state.T)
specificEnthalpy(      m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = h_T(m.data,state.T)
specificInternalEnergy(m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = h_T(m.data,state.T) - m.data.R*state.T
specificHeatCapacityCp(m::SingleGasNasa, state::ThermodynamicState_pT)::Float64 = s0_T(m.data,state.T) - m.data.R*log(state.p/m.infos.reference_p)

function standardCharacteristics(m::SingleGasNasa)::Dict{AbstractString,Any}
    p_ref = m.infos.reference_p
    T     = collect( range(m.fluidLimits.TMIN, stop=min(1600.0, m.fluidLimits.TMAX), length=501) )
    p     = [0.5e5, 1.0e5, 2.0e5]
    nT    = length(T)
    np    = length(p)
    h     = zeros(nT)
    u     = zeros(nT)
    cp    = zeros(nT)
    d     = zeros(nT,np)

    for i in 1:nT
        state = setState_pT(m,p_ref,T[i])
        h[i]  = specificEnthalpy(      m, state)
        u[i]  = specificInternalEnergy(m, state)
        cp[i] = specificHeatCapacityCp(m, state)
    end

    for j in 1:np
        for i in 1:nT
            d[i,j] = to_DensityDisplayUnit( density(m, ThermodynamicState_pT(p[j],T[i]) ) )
        end
    end

    mediumDict = Dict{AbstractString,Any}()
    mediumDict["T"]  = uconvert.(u"°C", T*1U"K")
    mediumDict["h"]  = h*1U"J/kg"
    mediumDict["u"]  = u*1U"J/kg"
    mediumDict["cp"] = cp*1U"J/(kg*K)"
    mediumDict["d(p=0.5 bar)"] = d[:,1]*1U"g/cm^3"
    mediumDict["d(p=1.0 bar)"] = d[:,2]*1U"g/cm^3"
    mediumDict["d(p=2.0 bar)"] = d[:,3]*1U"g/cm^3"
    return mediumDict
end

function standardPlot(m::SingleGasNasa; figure=1)
    mediumDict = standardCharacteristics(m)
    ModiaMath.plot(mediumDict, [("h", "u"), "cp", ("d(p=0.5 bar)" ,
                                                   "d(p=1.0 bar)" ,
                                                   "d(p=2.0 bar)")], xAxis="T", heading=m.infos.mediumName, figure=figure)
end


text:
"""
    MotifSequenceGenerator
This module generates random sequences of motifs, under the constrain that the
sequence has some total length ℓ so that `q - δq ≤ ℓ ≤ q + δq`.
All main functionality is given by the function [`random_sequence`](@ref).
"""
module MotifSequenceGenerator

using Combinatorics, Random, StatsBase

export random_sequence, all_possible_sums

struct DeadEndMotifs <: Exception
  tries::Int
  summands::Int
  tailcut::Int
end
Base.showerror(io::IO, e::DeadEndMotifs) = print(io,
"DeadEndMotifs: Couldn't find a proper sequence with $(e.tries) random tries, "*
"each with summands up to $(e.summands) (total tailcuts: $(e.tailcut)).")

#= Algorithm description
The algorithm works as follows: First a random sequence of motifs is created,
so that it has length of `q - δq ≤ ℓ ≤ q - δq + maximum(motiflengths)`.
The possible tries of random sequences is set by the `tries` keyword (default `5`).
The sequence is optionally sampled given a probability vector.

For each random try, it is first checked whether the sequence is already correct.
If not, the last entry of the sequence is dropped. Then, since the sequence is now
already smaller than `q`, all possible sums of `summands` out of the motif pool
are checked. If some combination of `summands` sums to the difference,
they are added to the sequence.
For multiple satisfactory combinations, a random one is picked.

If the random combination of `summands` does not fit, one more entry is dropped
from the sequence and the process is repeated.

The keyword `tailcut` limits how many times will an element
be dropped (default is `2`).
The keyword `summands` denotes how many possible combinations of sums to check.
Default is `3` which means that we check from 2 up to `3` summands (*all*
possible combinations of summands are checked, which means that the computation
skyrockets for large `summands`!).

If after going though all these combinations of possible sequences we do not find
a proper one, an error is thrown.
=#

"""
    random_sequence(motifs::Vector{M}, q, limits, translate, δq = 0; kwargs...)
Create a random sequence of motifs of type `M`, under the constraint that the
sequence has "length" `ℓ` **exactly** within `q - δq ≤ ℓ ≤ q + δq`.
Return the sequence itself as well as the
sequence of indices of `motifs` used to create it. A vector of probabilities `weights`
can be given as a keyword argument, which then dictates the sampling probability
for each entry of `motifs` for the initial sequence created.

"length" here means an abstracted length defined by the struct `M`,
based on the `limits` and `translate` functions.
It does **not** refer to the amount of elements!

`M` can be anything, given the two functions
* `limits(motif)` : Some function that given the `motif` it returns the
  `(start, fine)` of the the motif in the same units as `q`.
  This function establishes a measure of length, which simply is `fine - start`.
* `translate(motif, t)` : Some function that given the `motif` it returns a *new*
  motif which is translated by `t` (either negative or positive), with
  respect to the same units as `q`.

## Other Keywords
Please see the source code (use `@which`) for a full description of the algorithm.

* `tries = 5` : Up to how many initial random sequences are accepted.
* `taulcut = 2` : Up to how times an element is dropped from the initial guess.
* `summands = 3` : Up to how many motifs may be combined as a sum to
  complete a sequence.
"""
function random_sequence(motifs::Vector{M}, q,
    limits, translate, δq = 0;
    tries = 5, summands = 3, tailcut = 2,
    weights = ones(length(motifs))) where {M}

    ws = _toweight(weights)

    idxs = 1:length(motifs)
    motifs0, motiflens = _motifs_at_origin(motifs, limits, translate)

    q - δq < minimum(motiflens) && throw(ArgumentError(
    "Minimum length of motifs is greater than `q - δq`! "*
    "Impossible to make a sequence."
    ))

    eltype(motiflens) <: AbstractFloat && δq == 0 && throw(ArgumentError(
    "Due to finite precision of floating values, it is almost always impossible "*
    "to create a sequence of *exact* amount of floating length. Please provide "*
    "a δq > 0."
    ))

    worked = false; count = 0; seq = Int[]
    while worked == false
        count > tries && throw(DeadEndMotifs(tries, summands, tailcut))

        seq, seq_length = _random_sequence_try(motiflens, q, δq, ws)

        worked = _complete_sequence!(seq, motiflens, q, δq, summands, tailcut)

        count += 1
    end

    return _instantiate_sequence(motifs0, motiflens, seq, translate), seq
end

_toweight(a) = (s = sum(a); ProbabilityWeights(a./s, 1))

"""
    _motifs_at_origin(motifs, limits, translate) -> (motifs0, motiflens)
Bring all motifs to the origin and compute the motif lengths.
"""
function _motifs_at_origin(motifs::Vector{M}, limits, translate) where M
    motifs0 = similar(motifs)
    a, b = limits(motifs[1])
    motiflens = zeros(typeof(b-a), length(motifs))
    for i in 1:length(motifs)
        start, fine = limits(motifs[i])
        motifs0[i] = start == 0 ? motifs[i] : translate(motifs[i], -start)
        motiflens[i] = fine - start
    end
    return motifs0, motiflens
end

"""
    _random_sequence_try(motiflens, q, δq [, ws]) -> seq, seq_length
Return a random sequence of motif indices
so that the total sequence is *guaranteed* to have total length of
`q - δq ≤ ℓ ≤ q - δq + maximum(motiflens)`.
"""
function _random_sequence_try(motiflens, q, δq, ws = defaultweights(motiflens))
    seq = Int[]; seq_length = 0; idxs = 1:length(motiflens)
    while seq_length < q - δq
        i = sample(idxs, ws)
        push!(seq, i)
        seq_length += motiflens[i]
    end
    return seq, seq_length
end

function _complete_sequence!(seq, motiflens, q, δq, summands, tailcut)
    # notice the ||
    _complete_sequence_extra!(seq, motiflens, q, δq) ||
    _complete_sequence_remainder!(seq, motiflens, q, δq, summands, tailcut) ||
    false
end
function _complete_sequence_extra!(seq, motiflens, q, δq)

    ℓ = sum(motiflens[k] for k in seq)

    if q - δq ≤ ℓ ≤ q + δq
        # Case 0: The sequence is already within δq limits
        return true
    elseif (extra = ℓ - q - δq) > 0
        # Case 1: There is an extra difference. We check if it can be
        # accounted for by deleted some motif with length up to `extra + 2δq`.
        # We find the possible motifs, pick a random one, and pick
        # a random position in the sequence that it exists.
        # Delete that entry of the sequence.
        mi = findall(x -> extra ≤ x ≤ extra + 2δq, motiflens)
        possible = findall(in(mi), seq)
        if !isempty(possible)
            deleteat!(seq, rand(possible))
            return true
        end
    end
    return false
end
function _complete_sequence_remainder!(seq, motiflens, q, δq, summands, tailcut)

    # Case 2: Recursive deletion of last entry of the sequence, and trying to
    # see if it can be completed with some combination of existing motifs
    tcut = 0
    while tcut < tailcut
        tcut += 1
        pop!(seq)
        isempty(seq) && return false

        # At this point ℓ is guaranteed less than q - δq
        remainder = q - δq - sum(motiflens[k] for k in seq)
        @assert remainder > 0

        mi = findall(x -> remainder ≤ x ≤ remainder + 2δq, motiflens)
        if !isempty(mi)
            push!(seq, rand(mi))
            return true
        end

        for n in 2:summands
            everything = all_possible_sums(motiflens, n)
            sums = [e[1] for e in everything]
            cases = findall(x -> remainder ≤ x ≤ remainder + 2δq, sums)
            if !isempty(cases)
                idxs_of_vals = shuffle!(everything[rand(cases)][2])
                push!(seq, idxs_of_vals...)
                return true
            end
        end
    end
    return false
end

# Function provided by Mark Birtwistle in stackoverflow
"""
    all_possible_sums(summands, n)
Compute all possible sums from combining `n` elements from `summands`
(with repetition), only using unique combinations.

Return a vector of tuples: the first
entry of each tuple is the sum, while the second is the indices of summands
used to compute the sum.
"""
function all_possible_sums(summands::Vector{T}, n) where {T}
    m = length(summands)
    r = Vector{Tuple{T, Vector{Int}}}(undef, binomial(m + n - 1, n))
    s = with_replacement_combinations(eachindex(summands), n)
    for (j, i) in enumerate(s)
        r[j] = (sum(summands[j] for j in i), i)
    end
    return r
end

function _instantiate_sequence(motifs0::Vector{M}, motiflens, seq, translate) where M
    ret = M[]
    prev = 0
    for s in seq
        push!(ret, translate(motifs0[s], prev))
        prev += motiflens[s]
    end
    return ret
end

end#Module

ext: jl


text:
import Pkg

modelproto(;kwargs...) = ModelProto(;
    ir_version=6,
    opset_import=[OperatorSetIdProto(version=11)],
    producer_name="ONNX.jl",
    producer_version=string(Pkg.Types.Context().env.project.version), # TODO: Ugh....
    kwargs...)

"""
    graphproto()
Return an [`ONNX.GraphProto`](@ref) with all fields initialized to empty arrays.
"""
graphproto(;kwargs...) = GraphProto(;
    node = NodeProto[],
    initializer = TensorProto[],
    input = ValueInfoProto[],
    output = ValueInfoProto[],
    value_info = ValueInfoProto[],
    kwargs...
)

add!(gp::GraphProto, np::NodeProto) = push!(gp.node, np)

add!(gp::GraphProto, tp::TensorProto) = push!(gp.initializer, tp)

##############################################################################
#                                 Methods                                    #
##############################################################################

onnx_name(v::Variable) = "x$(v.id)"
onnx_name(op::Ghost.AbstractOp) = "x$(op.id)"

"""
    save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g::GraphProto, ::OpConfig{:Backend, Fn}, op::Ghost.Call)

Serialize a single operation from a tape to graph.
"""
function save_node!(g::GraphProto, op::Ghost.Call)
    save_node!(g, OpConfig{:ONNX, typeof(op.fn)}(), op)
end

# function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(conv)}, op::Ghost.Call)

#     nd = NodeProto(
#         input=[onnx_name(v) for v in op.args],
#         output=[onnx_name(op)],
#         name=onnx_name(op),
#         attribute=AttributeProto[],  # TODO
#         op_type="Conv"
#     )
#     push!(g.node, nd)
# end

function save_node!(g::GraphProto, ::OpConfig{:ONNX, typeof(add)}, op::Ghost.Call)
    nd = NodeProto(
        input=[onnx_name(v) for v in op.args],
        output=[onnx_name(op)],
        name=onnx_name(op),
        attribute=AttributeProto[],
        op_type="Add"
    )
    push!(g.node, nd)
end

ValueInfoProto(op::Ghost.AbstractOp) = ValueInfoProto(
    onnx_name(op),
    mrev(size(op.val)),
    eltype(op.val)
)

##############################################################################
#                                    API                                     #
##############################################################################

"""
    save(io::IO, tape::Ghost.Tape{ONNXCtx})
    save(filename::String, tape::Ghost.Tape{ONNXCtx})

Save tape as an ONNX model. The way a particular operation is serialized is
controlled by methods of [save_node!](@ref).

See also: [load!](@ref)
"""
function save(io::IO, tape::Tape{ONNXCtx})
    g = graphproto()
    g.name = "generated_model"
    for op in tape
        if op isa Ghost.Input
            # Ghost.Tape represents both - arguments and model parameters
            # as Input ops. In ONNX, parameters (and constants) must be added
            # to .initializer, so will we have to do when we understand how to
            # handle parameters (and complex structs) in general.
            # For reference, the following commented line shows how to add Input
            # to the .initializer:
            # add!(g, TensorProto(op.val, onnx_name(op)))
            push!(g.input, ValueInfoProto(op))
        elseif op isa Ghost.Call
            save_node!(g, op)
        else
            error("$(typeof(op)) is not yet supported in model export")
        end
    end
    push!(g.output, ValueInfoProto(tape[tape.result]))
    m = modelproto();
    m.graph = g;
    writeproto(io, m)
end

function save(filename::String, tape::Tape{ONNXCtx})
    open(filename, "w") do io
        save(io, tape)
    end
end


text:
using DelimitedFiles
include("TrafficAssign.jl")
link_flow, link_travel_time = solve_TAP("tmp", "../../Joint/tmp_jl/net.txt", "../../Joint/tmp_jl/trips.txt", [1,0,0,0,.15,0])

writedlm("../../Joint/tmp_jl/link_flow.txt", link_flow)
writedlm("../../Joint/tmp_jl/link_travel_time.txt", link_travel_time)

max_stars_repo_name: salomonw/tnet

max_issues_repo_name: salomonw/tnet

max_forks_repo_name: salomonw/tnet


text:
# https://storopoli.io/Bayesian-Julia/pages/3_prob_dist/
# https://nbviewer.org/github/takuizum/julia/blob/master/jmd/summaryPkg.html
using Plots, StatsPlots, LaTeXStrings
using Distributions, Random

# Random draw
Random.seed!(123)
d = DiscreteUniform(1, 6)
rand(d, 1)

# DiscreteUniform
d = DiscreteUniform(1, 6)
params(d)
params(d)[1]
params(d)[2]
span(d)
probval(d)
probval(d) * span(d)
minimum(d)
maximum(d)

plot(
    d,
    label="6-sided Dice",
    markershape=:circle,
    xlabel=L"\theta",
    ylabel="Probability Mass",
    ylims=(0, 0.3)
)

# Bernoulli
d = Bernoulli(0.1)
params(d)
succprob(d)
failprob(d)

plot(
    Bernoulli(0.5),
    markershape=:circle,
    label=L"p=0.5",
    alpha=0.5,
    xlabel=L"\theta",
    ylabel="Probability Mass",
    ylim=(0, 1)
)
plot!(
    Bernoulli(0.2),
    markershape=:circle,
    label=L"p=0.2",
    alpha=0.5
)

# Binomial
d = Binomial(5, 0.5)

plot(Binomial(5, 0.5),
        markershape=:circle,
        label=L"p=0.5",
        alpha=0.5,
        xlabel=L"\theta",
        ylabel="Probability Mass"
    )
plot!(Binomial(5, 0.2),
        markershape=:circle,
        label=L"p=0.2",
        alpha=0.5)

# Poisson
d = Distributions.Poisson(3)
params(d)
mean(d)

plot(Poisson(1),
        markershape=:circle,
        label=L"\lambda=1",
        alpha=0.5,
        xlabel=L"\theta",
        ylabel="Mass"
    )
plot!(Poisson(4),
    markershape=:circle,
    label=L"\lambda=4",
    alpha=0.5)

# NegativeBinomial
plot(NegativeBinomial(1, 0.5),
        markershape=:circle,
        label=L"k=1",
        alpha=0.5,
        xlabel=L"\theta",
        ylabel="Mass"
    )
plot!(NegativeBinomial(2, 0.5),
        markershape=:circle,
        label=L"k=2",
        alpha=0.5)

# Normal
d = Normal(0, 2)
params(d)
mean(d)
std(d)
var(d)

plot(
    Normal(0, 1),
    label=L"\sigma=1",
    lw=5,
    xlabel=L"\theta",
    ylabel="Probability Density",
    xlims=(-4, 4)
)
plot!(Normal(0, 0.5), label=L"\sigma=0.5", lw=5)
plot!(Normal(0, 2), label=L"\sigma=2", lw=5)

# LogNormal
d = LogNormal(0, 2)
params(d)
meanlogx(d)
varlogx(d)
stdlogx(d)

plot(LogNormal(0, 1),
        label=L"\sigma=1",
        lw=5,
        xlabel=L"\theta",
        ylabel="Density",
        xlims=(0, 3)
    )
plot!(LogNormal(0, 0.25), label=L"\sigma=0.25", lw=5)
plot!(LogNormal(0, 0.5), label=L"\sigma=0.5", lw=5)

# Exponential
plot(Exponential(1),
        label=L"\lambda=1",
        lw=5,
        xlabel=L"\theta",
        ylabel="Density",
        xlims=(0, 4.5)
    )
plot!(Exponential(0.5), label=L"\lambda=0.5", lw=5)
plot!(Exponential(1.5), label=L"\lambda=2", lw=5)

# Student-t
plot(TDist(2),
        label=L"\nu=2",
        lw=5,
        xlabel=L"\theta",
        ylabel="Density",
        xlims=(-4, 4)
    )
plot!(TDist(8), label=L"\nu=8", lw=5)
plot!(TDist(30), label=L"\nu=30", lw=5)

# Beta
plot(Beta(1, 1),
        label=L"a=b=1",
        lw=5,
        xlabel=L"\theta",
        ylabel="Density",
        xlims=(0, 1)
    )
plot!(Beta(3, 2), label=L"a=3, b=2", lw=5)
plot!(Beta(2, 3), label=L"a=2, b=3", lw=5)


text:
module MIDIProtocol

using MIDI
using MusicManipulations
using MPIMeasurements

include("Protocol.jl")

end # module


text:
using .Curve: BN, FP, EP, EP2, LIB, Limb
using .Util: Point
using .Config: PRIME, ORDER, @ID

evalpoly(poly, x) = evalpoly(poly, BN(x))
function evalpoly(poly, x::BN)
    @assert length(poly) > 0
    xpow, y = x, first(poly)
    for coeff in Iterators.drop(poly, 1)
        y += xpow * coeff
        xpow = mod(x * xpow, ORDER)
    end
    return isa(y, Point) ? y : mod(y, ORDER)
end

function lagrange_interpolate_c0(coeffs::LagrangeCoeffGenerator, shares::AbstractDict{@ID,V})  where {V<:Union{BN,EP,EP2}}
    length(shares) >= 2 || error("need at least two shares")
    c0 = zero(V)
    for (x, y) in shares
        c0 += coeffs[x] * y
    end
    return isa(c0, Point) ? c0 : mod(c0, ORDER)
end

function lagrange_interpolate_c0(weights::BarycentricWeightGenerator, shares::AbstractDict{@ID,V}) where {V<:Union{BN,EP,EP2}}
    length(shares) >= 2 || error("need at least two shares")
    num, denom = zero(V), zero(BN)
    for (x, y) in shares
        # optimisation: normally inv is done when calculating weights
        #w = invmod(weights[j], ORDER)
        w = invmod(weights[x], ORDER)
        #@info weights[j], w
        num += w * y
        # just addition in practice no risk of overflow
        denom += w
    end
    c0 = invmod(mod(denom, ORDER), ORDER) * num
    return isa(c0, Point) ? c0 : mod(c0, ORDER)
end

lagrange_interpolate_c0(shares::AbstractDict) =
    lagrange_interpolate_c0(LagrangeCoeffGenerator(keys(shares)), shares)


text:
using GaussianDistributions
using Test

include("gaussian.jl")
include("bivariate.jl")


text:
module SpherePacking

using ClusteredLowRankSolver, AbstractAlgebra, BasesAndSamples, SpecialFunctions

export Nsphere_packing, cohnelkies

function spherevolume(n, r)
    sqrt(BigFloat(pi))^n / gamma(BigFloat(n) / 2 + 1) * r^n
end

laguerre(k, alpha, x) = basis_laguerre(k, alpha, x)[end]

function Nsphere_packing(n,d,r,N=length(r),prec=512; ret_sdp = false, scaling_fun=x->1, kwargs...)
    #The problem:
    # min M
    # M - f_ii(0)>=0 for i=1:N
    # F(f)(0) - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N >=0
    # F(f)(t) ⪰ 0 for t>=0
    # -f(w)_ij >= 0 for w>=r_i+r_j for i=1:N j=1:i (symmetric)
    # With f(x) = sum_k a_k k!/π^k L_k^{n/2-1}(π ||x||^2)
    # and thus F(f)(x) = sum_k a_k x^{k}
    # where a_k are symmetric NxN matrices
    # min M
    # s.t.  - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N + \sum_i \sum_j a_{ij,0} E_ij >= 0 (G={1}) (NxN)
    #       0 + sum_k sum_i sum_j a_{ij,k} E_{ij} x^k >= 0        (G = {1,x}) (NxN)
    #       0 - sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) >=0          (G = {1,x - (r_i + r_j)^2}) (1x1)
    #       M - sum_k a_iik k!/pi^k L_k^{n/2-1}(0) >= 0 for i=1:N   G = {1}, 1x1
    setprecision(BigFloat,prec)
    R, (x,) = PolynomialRing(RealField,["x"])

    #Constraint 1: - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N + \sum_i \sum_j a_{ij,0} E_ij is PSD. So -a_ijk + PSD_ij = -(vol B(r_i) vol B(r_j)_^1/2
    con1 = []
    for i=1:N
        for j=1:i
            if i != j
                push!(con1,Constraint(-R(sqrt(spherevolume(n,r[i])*spherevolume(n,r[j]))),
                    Dict(Block(:PSD1,i,j)=>LowRankMatPol([R(1//2)],[[R(1)]]),
                        Block(:PSD1,j,i)=>LowRankMatPol([R(1//2)],[[R(1)]])),
                    Dict(Block(0,i,j)=>R(-1)), #we only use a_ijk for i>=j
                    [[0]]))
            else
                push!(con1,Constraint(-R(sqrt(spherevolume(n,r[i])*spherevolume(n,r[j]))),
                    Dict(Block(:PSD1,i,j)=>LowRankMatPol([R(1)],[[R(1)]])),
                    Dict(Block(0,i,j)=>R(-1)), #we only use a_ijk for i>=j
                    [[0]]))
            end
        end
    end

    # We orthogonalize the laguerre basis with respect to the sample points
    q = basis_laguerre(2d+1, BigFloat(n) / 2 - 1, BigFloat(2 * pi) * x)
    max_coef = [max(coefficients(q[i])...) for i = 1:length(q)]
    basis = [max_coef[i]^(-1) * q[i] for i = 1:length(q)]
    samples = sample_points_rescaled_laguerre(2d+1)
    basis, samples = approximatefekete(basis,samples)

    #constraint 2: sum_k a_{ij,k} x^k = <SOS21_ij, bb^T> + <SOS22_ij, xbb^T>
    #(Polynomial matrix constraint gives a polynomial constraint for each entry)
    con2 = []
    for i=1:N
        for j=1:i
            PSD_part = Dict()
            free_part = Dict()
            if i != j # We have to make the stuff symmetric. The factor still matters because we use the same SOS (l) block
                for k=0:2d+1
                    free_part[Block(k,i,j)] = -2*x^k
                end
                PSD_part[Block(:SOS21,i,j)] = LowRankMatPol([R(1)],[basis[1:d+1]])
                PSD_part[Block(:SOS22,i,j)] = LowRankMatPol([x],[basis[1:d+1]])
                PSD_part[Block(:SOS21,j,i)] = LowRankMatPol([R(1)],[basis[1:d+1]])
                PSD_part[Block(:SOS22,j,i)] = LowRankMatPol([x],[basis[1:d+1]])
            else
                for k=0:2d+1
                    free_part[Block(k,i,j)] = -x^k
                end
                PSD_part[Block(:SOS21,i,j)] = LowRankMatPol([R(1)],[basis[1:d+1]])
                PSD_part[Block(:SOS22,i,j)] = LowRankMatPol([x],[basis[1:d+1]])
            end
            push!(con2,Constraint(R(0),PSD_part,free_part,samples,[scaling_fun(s...) for s in samples]))
        end
    end

    #constraint 3: SOS + (x - (r_i + r_j)^2)*SOS + sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) = 0
    con3 = []
    for i=1:N
        for j=1:i
            #We take a new basis since the domain depends on the radii
            basis = [max_coef[i]^(-1) * q[i] for i = 1:length(q)]
            samples = [x .+ (r[i]+r[j])^2 for x in sample_points_rescaled_laguerre(2d+1)]
            basis,samples = approximatefekete(basis,samples)

            PSD_part = Dict()
            free_part = Dict()
            for k=0:2d+1
                free_part[Block(k,i,j)] = factorial(big(k))/BigFloat(pi)^k * laguerre(k, BigFloat(n) / 2 - 1, BigFloat(pi)*x)
            end
            #since we use 2d+1, we only need the constant of this SOS matrix; in practice it is always the 0 matrix
            PSD_part[Block((:SOS31,i,j))] = LowRankMatPol([R(1)],[basis[1:1]])
            PSD_part[Block((:SOS32,i,j))] = LowRankMatPol([x-(r[i]+r[j])^2],[basis[1:d+1]])
            push!(con3,Constraint(R(0),PSD_part,free_part,samples,[scaling_fun(s...) for s in samples]))
        end
    end

    #constraint 4: M - sum_k a_iik k!/pi^k L_k^{n/2-1}(0) = pos.slack
    con4 = []
    for i=1:N
        PSD_part = Dict()
        free_part = Dict()
        for k=0:2d+1
            free_part[Block(k,i,i)] = R(factorial(big(k))/BigFloat(pi)^k * laguerre(k, BigFloat(n) / 2 - 1,0))
        end
        free_part[:M] = -R(1)
        PSD_part[Block((:slack4,i))] = LowRankMatPol([R(1)],[[R(1)]])
        push!(con4,Constraint(R(0),PSD_part,free_part,[[0]]))
    end

    #objective: M
    obj = Objective(0,Dict(),Dict(:M=>1))

    sos = LowRankPolProblem(false,obj,[con1...,con2...,con3...,con4...])

    sdp = ClusteredLowRankSDP(sos)
    if ret_sdp
        return sdp
    else
        return solvesdp(sdp; kwargs...)
    end
end

function cohnelkies(n,d,r=1;scale_var_fun = k->1, prec=512, model_prec=prec,ret_sdp = false, scaling_fun=x->1,  kwargs...)
    #The problem:
    # min M
    # M - f_ii(0)>=0 for i=1:N
    # F(f)(0) - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N >=0
    # F(f)(t) >=0 for t>=0
    # -f(w)_ij >= 0 for w>=r_i+r_j for i=1:N j=1:i (symmetric)
    # With f(x) = sum_k a_k k!/pi^k L_k^{n/2-1}(pi ||x||^2)
    # and thus F(f)(x) = sum_k a_k x^{k}
    # min M
    # s.t.  - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N + \sum_i \sum_j a_{ij,0} E_ij >= 0 (G={1}) (NxN)
    #       0 + sum_k sum_i sum_j a_{ij,k} E_{ij} x^k >= 0        (G = {1,x}) (NxN)
    #       0 - sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) >=0          (G = {1,x - (r_i + r_j)^2}) (1x1)
    #       M - sum_k a_iik k!/pi^k L_k^{n/2-1}(0) >= 0 for i=1:N   G = {1}, 1x1
    # CohnElkies.jl models the problem for 2d+1, and then one of the PSD matrices is zero. We do that now too.
    setprecision(BigFloat,model_prec)

    R, (x,) = PolynomialRing(RealField,["x"])

    q = basis_laguerre(2d+1, BigFloat(n) / 2 - 1, BigFloat(2 * pi) * x)
    max_coef = [max(coefficients(q[i])...) for i = 1:length(q)]
    basis = [max_coef[i]^(-1) * q[i] for i = 1:length(q)]
    samples = vcat(sample_points_rescaled_laguerre(2d+1),sample_points_chebyshev(4d,0,5)) #TODO test with different sample points. Now they're fixed already and we only improve the basis
    basis, samples = approximatefekete(basis,samples)

    #constraint 1: sum_k a_k x^k = <SOS21_ij, bb^T> + <SOS22_ij, xbb^T>
    PSD_part = Dict()
    free_part = Dict()
    for k=1:2d+1
        free_part[Block(k)] = -scale_var_fun(k)*x^k
    end
    PSD_part[Block(:SOS21)] = LowRankMatPol([R(1)],[basis[1:d+1]])
    PSD_part[Block(:SOS22)] = LowRankMatPol([x],[basis[1:d+1]])

    con1 = Constraint(R(1),PSD_part,free_part,samples, [scaling_fun(s...) for s in samples]) #y_0 = 1 = constant

    #constraint 2: SOS + (x - (r_i + r_j)^2)*SOS + sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) = 0
    basis = [max_coef[i]^(-1) * q[i] for i = 1:length(q)] #2d+1 sample points, because we have polynomials of that degree
    samples = vcat([x .+ r^2 for x in sample_points_rescaled_laguerre(2d+1)],sample_points_chebyshev(4d,r^2,5+r^2))
    basis, samples = approximatefekete(basis,samples)

    PSD_part = Dict()
    free_part = Dict()
    for k=1:2d+1
        free_part[Block(k)] = scale_var_fun(k)*factorial(big(k))/BigFloat(pi)^k * laguerre(k, BigFloat(n) / 2 - 1, BigFloat(pi)*x)
    end
    # because we use 2d+1 points, we need the matrix elements to form a basis of polynomials up to degree 2d+1.
    # only taking the second matrix gives a basis of polynomials of the form (x-r^2)*p(x) with p of degree <= 2d, so we need a constant to form the complete basis.
    # one would think that taking 1 less sample point should make it possible to not use the 1x1 PSD block (2d basis polynomials and 2d points). However, then it gives the wrong answer (~0.133 instead of 0.779)
    #I am not sure why that is the case. the extra constant is nearly 0 as expected
    #for three blocks:
    PSD_part[Block((:SOS31))] = LowRankMatPol([R(1)],[basis[1:1]])
    #for 4 blocks:
    # PSD_part[Block((:SOS31))] = LowRankMatPol([R(1)],[basis[1:d+1]])
    PSD_part[Block((:SOS32))] = LowRankMatPol([x-r^2],[basis[1:d+1]])
    con2 = Constraint(-laguerre(0, BigFloat(n) / 2 - 1, BigFloat(pi)*x),PSD_part,free_part,samples,[scaling_fun(s...) for s in samples])

    #objective: M
    obj = Objective(spherevolume(n, BigFloat(r) / 2)*laguerre(0, BigFloat(n) / 2 - 1, BigFloat(0)),Dict(),Dict(Block(k)=> scale_var_fun(k)*spherevolume(n, BigFloat(r) / 2)*factorial(big(k))/BigFloat(pi)^k * laguerre(k, BigFloat(n) / 2 - 1, BigFloat(0)) for k=1:2d+1))
    #NOTE: these numbers become extremely large for large k. So no wonder that solvers have issues with that

    sos = LowRankPolProblem(false, obj, [con1,con2])
    sdp = ClusteredLowRankSDP(sos)
    if ret_sdp
        return sdp
    else
        return status, sol, time, errorcode = solvesdp(sdp;prec=prec, kwargs...) #status, sol, time
    end
end

end # of module


text:
include("_defaultcase.jl")
include("_delimcase.jl")
include("_purgecase.jl")

include("_decamelize.jl")


text:
using ProgressBars, Plots, StatsBase, Statistics, LaTeXStrings, JLD

@inline function RandomWalker2D(TotalSteps, R₀)
    choices = ((1,0),(0,1),(-1,0),(0,-1))
    TotalWalks = rand(choices,TotalSteps)
    return reduce(.+,TotalWalks) .+ R₀
end

#=
function PlotHeatmap(Data, RWNumbers) #useless piece of shit
    Tuples = collect(Tuple(Data[i,:]) for i ∈ 1:RWNumbers)
    maxindx = max(getindex.(Tuples, 1)...,getindex.(Tuples, 2)...)
    minindx = min(getindex.(Tuples, 1)...,getindex.(Tuples, 2)...)
    delta = maxindx - minindx + 1
    Surf = zeros(Int,delta,delta)

    for tup in Tuples
        Surf[(tup .- (minindx - 1))...] += 1
    end
    return heatmap(Surf)
end
=#

TStepsList = [10,100,1000]
RWNumbersList = [10000,100000,1000000]
R₀ = [0,0]
Data = []
for TotalSteps in TStepsList
    RWNumbers = TotalSteps * 1000
    TotalWalkers = zeros(Int,RWNumbers,2)
    for i ∈ ProgressBar(1:RWNumbers)
        TotalWalkers[i,:] = RandomWalker2D(TotalSteps,R₀)
    end
    push!(Data, copy(TotalWalkers))
end

cd(dirname(@__FILE__))
save("../../Data/Q1/Q1-2DHist.jld", "Data", Data)
# Data = load("../../Data/Q1/Q1-2DHist.jld")["Data"]
#Plot section

scatter(RgList.^2)

# PlotHeatmap(Data, RWNumbers)
dimnum = 3
P1 = histogram2d(collect(Tuple(Data[dimnum][i,:]) for i ∈ 1:RWNumbersList[dimnum]),
    nbin = 200,
    xaxis = nothing,
    yaxis = nothing,
    framestyle = :box
    )
P2 = histogram(Data[dimnum][:,2], orientation = :horizontal, xticks = [0,3000], bins = 150,
    ylabel = L"Y", flip = true,framestyle = :box, c = :purple)
P3 = histogram(Data[dimnum][:,1], flip = true, yticks = [0,3000], bins = 150,
    xlabel = L"X", framestyle = :box, c = :purple)
P4 = plot(legend=false,grid=false,
    foreground_color_subplot=:white,background=:white,framestyle = :box)

P21 = plot(P2,P1,layout = @layout[a{0.1w} grid(1,1, heights = [0.95])], legend = nothing)

P22 = plot(P4,P3,P4,layout = @layout[a{0.08w} grid(1,1) a{0.02w}], legend = nothing)

plot(P21,P22,layout = grid(2,1, heights = [0.85,0.15]),
    plot_title = L"Distribution\ of\ %$(RWNumbersList[dimnum])\ RW\ in\ %$(TStepsList[dimnum])\ steps")

savefig("../../Figs/Q1/Q1-2DHist$dimnum.pdf")


text:
if VERSION >= v"1.3"
    const safe_rand = rand
else
    # rand is thread safe after julia 1.3
    # TODO: delete this when we decide to drop 1.0 compatibility

    # --------------------------------------------------------------------
    # rand() is not threadsafe (https://discourse.julialang.org/t/4683)

    # Because we only require random numbers to sample parameters
    # and not the actual expensive computation, this seems like a better
    # approach than using separate RNG per thread.
    const rand_mutex = Ref{Threads.Mutex}()

    # constant overhead of about 80 ns compared to unsafe rand
    function safe_rand(args...)
        lock(rand_mutex[])
        result = rand(args...)
        unlock(rand_mutex[])
        result
    end
end

# issue #39
if all(map((FixedPoint, Normed)) do T
        isempty(methodswith(AbstractArray{<:T}, Interpolations.tweight))
    end)
    # Reached when upstream package (e.g., FixedPoint) doesn't provide a specialization
    Interpolations.tweight(A::AbstractArray{T}) where T<:FixedPoint = T
end

# backward support for OffsetArrays v0
if isdefined(OffsetArrays, :IdOffsetRange)
    # introduced by OffsetArrays v1.0.0
    OffsetRange = Union{UnitRange, IdentityUnitRange, OffsetArrays.IdOffsetRange}
else
    OffsetRange = Union{UnitRange, IdentityUnitRange}
end


text:
type SolverParams
nvar::Int64
dtsav::Float64
h0::Float64
h0nom::Float64
hmin::Float64
eps::Float64

    function SolverParams(nvar=0)
        this = new()
        this.nvar = nvar;
        this.dtsav = 1e-3;
        this.h0 = 1e-3;
        this.h0nom = 1e-3;
        this.hmin = 1e-16;
        this.eps = 1e-8;
        return this
    end
end


text:
module ProcessInputFile

using TOML
using Random
using Statistics
using LinearAlgebra
using Logging
using LibGit2
using Serialization
using Parameters
using Pkg

using ..UnitCells: UnitCell
using ..Lattices: Lattice
using ..Models: HolsteinModel, SSHModel, AbstractModel
using ..Models: assign_μ!, assign_ω!, assign_λ!, assign_λ₂!, assign_ω₄!
using ..Models: assign_t!, assign_ωᵢⱼ!, assign_hopping!
using ..Models: assign_datafolder!
using ..Models: initialize_model!, update_model!, read_phonons!, mulM!, mulMᵀ!
using ..MuFinder: MuTuner, update_μ!
using ..GreensFunctions: EstimateGreensFunction, update!
using ..InitializePhonons: init_phonons_half_filled!
using ..LangevinDynamics: EulerDynamics, RungeKuttaDynamics, HeunsDynamics
using ..HMC: HybridMonteCarlo
using ..SpecialUpdates: SpecialUpdate, NullUpdate, ReflectionUpdate, SwapUpdate
using ..FourierAcceleration: FourierAccelerator, update_Q!, update_M!
using ..SimulationParams: SimulationParameters
using ..Measurements: initialize_measurements_container, initialize_measurement_folders!
using ..KPMPreconditioners: LeftRightKPMPreconditioner, SymmetricKPMPreconditioner

export process_input_file, initialize_holstein_model

function process_input_file(filename::String,input::Dict)

    # Input file must describe either a Langevin or HMC Simulation but not both.
    @assert haskey(input,"hmc") ⊻ haskey(input,"langevin")

    # Input file must describe either a Holstein or SSH model but not both.
    @assert haskey(input,"holstein") ⊻ haskey(input,"ssh")

    ##################################
    ## DEFINE SIMULATION PARAMETERS ##
    ##################################

    # initialize simulation parameters
    sim_params = initialize_simulation_params(input)

    # copy input file into data folder
    cp(filename, joinpath(sim_params.datafolder,filename))

    ######################
    ## INITIALIZE MODEL ##
    ######################

    model = initialize_model(input)

    ##############################
    ## INITIALIZE PHONON FIELDS ##
    ##############################

    initialize_phonon_fields!(input,model)

    #####################################
    ## TUNE DENSITY/CHEMICAL POTENTIAL ##
    #####################################

    # initialize MuTuner
    μ_tuner = initialize_mutuner(input,model)

    ###########################
    ## DEFINE PRECONDITIONER ##
    ###########################

    preconditioner = initialize_preconditioner(input,model)

    #################################
    ## DEFINE FOURIER ACCELERATION ##
    #################################

    # initialize FourierAccelerator
    fa = initialize_fourieraccelerator(input, model)

    #####################
    ## DEFINE DYNAMICS ##
    #####################

    burnin_dynamics, simulation_dynamics = initialize_dynamics(input,model)

    ############################
    ## DEFINE SPECIAL UDPATES ##
    ############################

    burnin_reflect_update, sim_reflect_update = initialize_reflect_update(input,model,burnin_dynamics,simulation_dynamics)
    burnin_swap_update,    sim_swap_update    = initialize_swap_update(input,model,burnin_dynamics,simulation_dynamics)

    #########################
    ## DEFINE MEASUREMENTS ##
    #########################

    # construct object of estimating Green's function
    Gr = EstimateGreensFunction(model,input["measurements"]["num_random_vectors"])

    # construct measurements container
    container = initialize_measurements_container(model,input["measurements"],sim_params.datafolder)

    # initialize measurement files
    initialize_measurement_folders!(container)

    #################################
    ## INITIALIZE SIMULATION STATS ##
    #################################

    sim_stats = initialize_sim_stats()

    return (model, Gr, μ_tuner, sim_params, simulation_dynamics, burnin_dynamics,
            burnin_reflect_update, sim_reflect_update,
            burnin_swap_update, sim_swap_update,
            fa, preconditioner, container, 1, 1, sim_stats)
end

function process_checkpoint(input::Dict)

    ############################
    ## DESERIALIZE CHECKPOINT ##
    ############################

    chkpnt = deserialize(joinpath(input["simulation"]["datafolder"],"checkpoint.jls"))

    #######################
    ## UNPACK CHECKPOINT ##
    #######################

    @unpack model, μ_tuner, container, burnin_start, sim_start, sim_stats = chkpnt

    ######################################
    ## INITIALIZE SIMULATION PARAMETERS ##
    ######################################

    sim_params = initialize_simulation_params(input)

    #########################
    ## INITIALIZE DYNAMICS ##
    #########################

    burnin_dynamics, simulation_dynamics = initialize_dynamics(input,model,burnin_start,sim_start)

    ###########################
    ## DEFINE SPECIAL UDPATE ##
    ###########################

    burnin_reflect_update, sim_reflect_update = initialize_reflect_update(input,model,burnin_dynamics,simulation_dynamics)
    burnin_swap_update,    sim_swap_update    = initialize_swap_update(input,model,burnin_dynamics,simulation_dynamics)

    ###############################
    ## INITIALIZE PRECONDITIONER ##
    ###############################

    preconditioner = initialize_preconditioner(input,model)

    ##########################################
    ## INITIALIZE GREENS FUNCTION ESTIMATOR ##
    ##########################################

    Gr = EstimateGreensFunction(model,input["measurements"]["num_random_vectors"])

    ###################################
    ## INITIALIZE FOUIER ACCELERATOR ##
    ###################################

    fa = initialize_fourieraccelerator(input, model)

    return (model, Gr, μ_tuner, sim_params, simulation_dynamics, burnin_dynamics,
            burnin_reflect_update, sim_reflect_update,
            burnin_swap_update, sim_swap_update,
            fa, preconditioner, container, burnin_start, sim_start, sim_stats)
end

###########################################################
## METHODS INITIALIZING SIMULATION THAT ARE NOT EXPORTED ##
###########################################################

"""
Initialize a hamiltonian.
"""
function initialize_model(input::Dict)

    if haskey(input,"holstein") && haskey(input,"ssh")
        error("Config file cannot include both ssh and holstein tables.")
    end

    if !haskey(input,"holstein") && !haskey(input,"ssh")
        error("No valid model table detected in config file.")
    end

    # initialize random number generator
    rng = initialize_rng(input)

    if haskey(input,"holstein")
        model = initialize_holstein_model(input,rng)
    elseif haskey(input,"ssh")
        model = initialize_ssh_model(input,rng)
    else
        error("Neither holstein or ssh model defined.")
    end

    # assign data folder to model
    assign_datafolder!(model,input["simulation"]["datafolder"])

    return model
end

"""
Initialize Holstein Model from config file.
"""
function initialize_holstein_model(input::Dict,rng::AbstractRNG)

    # define lattice geometry
    unit_cell = UnitCell(input["lattice"]["ndim"],
                         input["lattice"]["norbits"],
                         hcat(input["lattice"]["lattice_vectors"]...),
                         hcat(input["lattice"]["basis_vectors"]...))

    # define lattice
    lattice = Lattice(unit_cell, input["lattice"]["L"])

    # restart for GMRES solver
    if haskey(input["solver"],"restart")
        restart = input["solver"]["restart"]
    else
        restart = 20
    end

    # initialize holstein model
    holstein = HolsteinModel(lattice,
                             input["holstein"]["beta"],
                             input["holstein"]["dtau"],
                             is_complex      = false,
                             iterativesolver = input["solver"]["type"],
                             rng             = rng,
                             tol             = input["solver"]["tol"],
                             maxiter         = input["solver"]["maxiter"],
                             restart         = restart)

    # adding phonon frequencies
    if haskey(input["holstein"],"omega")
        for d in input["holstein"]["omega"]
            stddev = 0.0
            if haskey(d,"stddev")
                stddev = d["stddev"]
            end
            for orbit in d["orbit"]
                assign_ω!(holstein,d["val"],stddev,orbit)
            end
        end
    end

    # adding chemical potential
    if haskey(input["holstein"],"mu")
        for d in input["holstein"]["mu"]
            stddev = 0.0
            if haskey(d,"stddev")
                stddev = d["stddev"]
            end
            for orbit in d["orbit"]
                assign_μ!(holstein,d["val"],stddev,orbit)
            end
        end
    end

    # check in anharmic term defined
    if haskey(input["holstein"],"omega4")
        # adding anharmoic term to holstein model
        for d in input["holstein"]["omega4"]
            stddev = 0.0
            if haskey(d,"stddev")
                stddev = d["stddev"]
            end
            for orbit in d["orbit"]
                assign_ω₄!(holstein,d["val"],stddev,orbit)
            end
        end
    end

    # check if any hopping defined
    if haskey(input["holstein"],"t")
        for t in input["holstein"]["t"]
            stddev = 0.0
            if haskey(t,"stddev")
                stddev = t["stddev"]
            end
            assign_t!(holstein, t["val"], stddev, t["orbit"][1], t["orbit"][2], Vector{Int}(t["dL"]))
        end
    end

    # adding electron-phonon coupling
    if haskey(input["holstein"],"lambda")
        for d in input["holstein"]["lambda"]
            stddev = 0.0
            if haskey(d,"stddev")
                stddev = d["stddev"]
            end
            for orbit in d["orbit"]
                assign_λ!(holstein,d["val"],stddev,orbit)
            end
        end
    end

    # adding electron-phonon coupling
    if haskey(input["holstein"],"lambda2")
        for d in input["holstein"]["lambda2"]
            stddev = 0.0
            if haskey(d,"stddev")
                stddev = d["stddev"]
            end
            for orbit in d["orbit"]
                assign_λ₂!(holstein,d["val"],stddev,orbit)
            end
        end
    end

    # initialize model
    initialize_model!(holstein)

    return holstein
end

"""
Initialize SSH model from config file.
"""
function initialize_ssh_model(input::Dict,rng::AbstractRNG)

    # define lattice geometry
    unit_cell = UnitCell(input["lattice"]["ndim"],
                         input["lattice"]["norbits"],
                         hcat(input["lattice"]["lattice_vectors"]...),
                         hcat(input["lattice"]["basis_vectors"]...))

    # define lattice
    lattice = Lattice(unit_cell, input["lattice"]["L"])

    # restart for GMRES solver
    if haskey(input["solver"],"restart")
        restart = input["solver"]["restart"]
    else
        restart = 20
    end

    # initialize holstein model
    ssh = SSHModel(lattice,
                   input["ssh"]["beta"],
                   input["ssh"]["dtau"],
                   is_complex      = false,
                   iterativesolver = input["solver"]["type"],
                   rng             = rng,
                   tol             = input["solver"]["tol"],
                   maxiter         = input["solver"]["maxiter"],
                   restart         = restart)

    # adding chemical potential
    for d in input["ssh"]["mu"]
        stddev = 0.0
        if "stddev" in keys(d)
            stddev = d["stddev"]
        end
        for orbit in d["orbit"]
            assign_μ!(ssh,d["val"],stddev,orbit)
        end
    end

    # add hoppings
    if haskey(input["ssh"],"hopping")
        for d in input["ssh"]["hopping"]
            if haskey(d,"t_avg")
                t = d["t_avg"]
            else
                t = 0.0
            end
            if haskey(d,"t_std")
                σt = d["t_std"]
            else
                σt = 0.0
            end
            if haskey(d,"alpha_avg")
                α = d["alpha_avg"]
            else
                α = 0.0
            end
            if haskey(d,"alpha_std")
                σα = d["alpha_std"]
            else
                σα = 0.0
            end
            if haskey(d,"alpha2_avg")
                α₂ = d["alpha2_avg"]
            else
                α₂ = 0.0
            end
            if haskey(d,"alpha2_std")
                σα₂ = d["alpha2_std"]
            else
                σα₂ = 0.0
            end
            if haskey(d,"omega_avg")
                ω = d["omega_avg"]
            else
                ω = 0.0
            end
            if haskey(d,"omega_std")
                σω = d["omega_std"]
            else
                σω = 0.0
            end
            if haskey(d,"omega4_avg")
                ω₄ = d["omega4_avg"]
            else
                ω₄ = 0.0
            end
            if haskey(d,"omega4_std")
                σω₄ = d["omega4_std"]
            else
                σω₄ = 0.0
            end
            if haskey(d,"name")
                name = d["name"]
            else
                name = ""
            end
            o₁ = d["orbits"][1]
            o₂ = d["orbits"][2]
            dL = zeros(Int,3)
            dL[1:length(d["dL"])] .= d["dL"]
            assign_hopping!(ssh,t,σt,ω,σω,ω₄,σω₄,α,σα,α₂,σα₂,o₁,o₂,dL,name)
        end
    end

    # initialize model
    initialize_model!(ssh)

    return ssh
end

"""
Initialize Phonon Fields.
"""
function initialize_phonon_fields!(input::Dict,model::AbstractModel)

    # determine type of el-ph hamiltonian
    if haskey(input,"holstein")
        model_type = "holstein"
    else
        model_type= "ssh"
    end

    # intialize phonon field
    if !haskey(input[model_type],"read_phonon_config")
        input[model_type]["read_phonon_config"] = false
    end
    if input[model_type]["read_phonon_config"] # read in phonon field
        phononfile = input[model_type]["phonon_config_file"]
        read_phonons!(model, phononfile)
        cp(phononfile, joinpath(input["simulation"]["datafolder"],basename(phononfile)) )
    else # initialize to random phonon field
        init_phonons_half_filled!(model)
    end

    return nothing
end

"""
Initialize Preconditioner.
"""
function initialize_preconditioner(input::Dict,model::AbstractModel)

    # reconstruct preconditioner
    if haskey(input["solver"],"preconditioner")

        if haskey(input["solver"]["preconditioner"],"n")
            n = input["solver"]["preconditioner"]["n"]
        else
            n = 20
        end

        if haskey(input["solver"]["preconditioner"],"buf")
            buf = input["solver"]["preconditioner"]["buf"]
        else
            buf = 0.05
        end

        if haskey(input["solver"]["preconditioner"],"c1")
            c1 = input["solver"]["preconditioner"]["c1"]
        else
            c1 = 1.0
        end

        if haskey(input["solver"]["preconditioner"],"c2")
            c2 = input["solver"]["preconditioner"]["c2"]
        else
            c2 = 1.0
        end

        if lowercase(input["solver"]["type"])=="cg"
            preconditioner = SymmetricKPMPreconditioner(model,n,buf,c1,c2)
        else
            preconditioner = LeftRightKPMPreconditioner(model,n,buf,c1,c2)
        end
    else

        preconditioner = I
    end

    return preconditioner
end

"""
Initialize FourierAccelerator.
"""
function initialize_fourieraccelerator(input::Dict, model::AbstractModel)

    # defining FourierAccelerator type
    fa = FourierAccelerator(model)

    # set the mass used to construct fourier acceleration matrix
    for d in input["fourier_acceleration"]
        mass = d["mass"]
        if haskey(d,"c")
            c = d["c"]
        else
            c = 0.0
        end
        update_Q!(fa, model, d["omega_min"], d["omega_max"], mass)
        update_M!(fa, model, d["omega_min"], d["omega_max"], mass, c)
    end

    return fa
end

"""
Initialize SimulationParams.
"""
function initialize_simulation_params(input::Dict)

    if haskey(input,"hmc")
        meas_freq = input["hmc"]["meas_freq"]
        nsteps    = input["hmc"]["simulation_updates"]
        burnin    = input["hmc"]["burnin_updates"]
    else
        @assert input["langevin"]["burnin_timesteps"]%input["langevin"]["meas_freq"]==0
        meas_freq = input["langevin"]["meas_freq"]
        nsteps    = input["langevin"]["simulation_timesteps"]
        burnin    = input["langevin"]["burnin_timesteps"]
    end

    # set default checkpointing frequency in minutes
    if !haskey(input["simulation"],"checkpoint_freq")
        input["simulation"]["checkpoint_freq"]=10
    end

    # construct simulation parameters object.
    sim_params = SimulationParameters(burnin,
                                      nsteps,
                                      meas_freq,
                                      input["simulation"]["num_bins"],
                                      input["simulation"]["checkpoint_freq"],
                                      input["simulation"]["filepath"],
                                      input["simulation"]["foldername"],
                                      input["simulation"]["datafolder"])

    # make direcotory data will be written to
    if !isdir(sim_params.datafolder)
        mkdir(sim_params.datafolder)
    end

    # create log for simulation
    logfilename = joinpath(sim_params.datafolder, "$(sim_params.foldername).log")
    logio       = open(logfilename,"w+")
    logger      = SimpleLogger(logio)
    global_logger(logger)

    # log version info
    Pkg.status("ElPhDynamics",io=logger.stream)
    @info("Julia Version: $VERSION")
    flush(logger.stream)

    return sim_params
end

"""
Initialize Random Number Generator with Seed
"""
function initialize_rng(input::Dict)

    # initialize random number generator with seed
    if !haskey(input["simulation"],"random_seed")
        input["simulation"]["random_seed"] = abs(rand(Int))
    end
    seed = input["simulation"]["random_seed"]
    rng  = MersenneTwister(seed)

    # write random seed to log file
    @info("Random Seed: $seed")
    logger = global_logger()
    flush(logger.stream)

    return rng
end

"""
Initialize MuTuner.
"""
function initialize_mutuner(input::Dict,model::AbstractModel)

    if haskey(input,"tune_density")
        logfile        = joinpath(input["simulation"]["datafolder"],"mu_tuner_log.out")
        targed_density = input["tune_density"]["density"]
        memory         = input["tune_density"]["memory"]
        κ_min          = input["tune_density"]["kappa_min"]
        μ_tuner = MuTuner(true, mean(model.μ), targed_density*model.Nsites, model.Nsites, model.β, model.Δτ, memory, κ_min*model.Nsites, logfile)
    else
        μ_tuner = MuTuner(false, mean(model.μ), 1.0*model.Nsites, model.Nsites, model.β, model.Δτ, 0.75, 0.1, "")
    end

    return μ_tuner
end

"""
Initialize Dynamics.
"""
function initialize_dynamics(input::Dict,model::AbstractModel,burnin_start::Int=1,sim_start::Int=1)

    # check to bugs
    if haskey(input,"hmc") && haskey(input,"langevin")
        error("Config file cannot include both hmc and langevin tables.")
    end

    # number of degrees of freedom (phonon fields) to simulate
    NL = length(model)

    if haskey(input,"hmc")

        Δt              = input["hmc"]["dt"]
        tr              = input["hmc"]["trajectory_time"]
        α               = input["hmc"]["momentum_conservation_fraction"]
        Nb              = input["hmc"]["num_multitimesteps"]

        # log file instructions
        if haskey(input["hmc"],"log")
            log = input["hmc"]["log"]
        else
            log = false
        end
        if log && haskey(input["hmc"],"verbose")
            verbose = input["hmc"]["verbose"]
        else
            verbose = false
        end
        hmc_simulation_logfile = joinpath(input["simulation"]["datafolder"],"hmc_sim_log.out")
        hmc_burnin_logfile     = joinpath(input["simulation"]["datafolder"],"hmc_burnin_log.out")

        @assert 0.0 <= α < 1.0
        simulation_dynamics = HybridMonteCarlo(model, Δt, tr, α, Nb, log=log, verbose=verbose,
                                               logfilename=hmc_simulation_logfile,updates=sim_start)

        # defining burnin dynamics
        if haskey(input["hmc"], "burnin")
            if haskey(input["hmc"]["burnin"],"dt")
                Δt = input["hmc"]["burnin"]["dt"]
            end
            if haskey(input["hmc"]["burnin"],"trajectory_time")
                tr = input["hmc"]["burnin"]["trajectory_time"]
            end
            if haskey(input["hmc"]["burnin"],"momentum_conservation_fraction")
                α = input["hmc"]["burnin"]["momentum_conservation_fraction"]
            end
            if haskey(input["hmc"]["burnin"],"num_multitimesteps")
                Nb = input["hmc"]["burnin"]["num_multitimesteps"]
            end
            @assert 0.0 <= α < 1.0
        end
        burnin_dynamics = HybridMonteCarlo(simulation_dynamics, Δt, tr, α, Nb, log=log, verbose=verbose,
                                           logfilename=hmc_burnin_logfile, updates=burnin_start)

    elseif input["langevin"]["update_method"]==1

        Δt = input["langevin"]["dt"]
        simulation_dynamics = EulerDynamics(model,Δt)
        burnin_dynamics = simulation_dynamics

    elseif input["langevin"]["update_method"]==2

        Δt = input["langevin"]["dt"]
        simulation_dynamics = RungeKuttaDynamics(model,Δt)
        burnin_dynamics = simulation_dynamics

    elseif input["langevin"]["update_method"]==3

        Δt = input["langevin"]["dt"]
        simulation_dynamics = HeunsDynamics(model,Δt)
        burnin_dynamics = simulation_dynamics

    end

    return burnin_dynamics, simulation_dynamics
end

"""
Initialize Reflection Update.
"""
function initialize_reflect_update(input::Dict,model::AbstractModel,burnin_dynaimcs,simulation_dynamics)

    # define special update for simulation updates
    if haskey(input,"langevin")
        sim_reflect_update = NullUpdate()
    else
        if haskey(input,"holstein") && haskey(input["hmc"],"reflection_update")
            freq   = input["hmc"]["reflection_update"]["freq"]
            nsites = input["hmc"]["reflection_update"]["nsites"]
            sim_reflect_update = ReflectionUpdate(model,freq,nsites)
        else
            sim_reflect_update = NullUpdate()
        end
    end

    # define special update for burnin updates
    burnin_reflect_update = sim_reflect_update
    if haskey(input,"hmc") && haskey(input,"holstein")
        if haskey(input["hmc"],"burnin")
            if haskey(input["hmc"]["burnin"],"reflection_update")
                freq   = input["hmc"]["reflection_update"]["freq"]
                nsites = input["hmc"]["reflection_update"]["nsites"]
                burnin_reflect_update = ReflectionUpdate(model,freq,nsites)
            end
        end
    end

    return burnin_reflect_update, sim_reflect_update
end

"""
Initialize Swap Update.
"""
function initialize_swap_update(input::Dict,model::AbstractModel,burnin_dynaimcs,simulation_dynamics)

    # define special update for simulation updates
    if haskey(input,"langevin")
        sim_swap_update = NullUpdate()
    else
        if haskey(input,"holstein") && haskey(input["hmc"],"swap_update")
            freq   = input["hmc"]["swap_update"]["freq"]
            nbonds = input["hmc"]["swap_update"]["nbonds"]
            sim_swap_update = SwapUpdate(model,freq,nbonds)
        else
            sim_swap_update = NullUpdate()
        end
    end

    # define special update for burnin updates
    burnin_swap_update = sim_swap_update
    if haskey(input,"hmc") && haskey(input,"holstein")
        if haskey(input["hmc"],"burnin")
            if haskey(input["hmc"]["burnin"],"swap_update")
                freq   = input["hmc"]["swap_update"]["freq"]
                nbonds = input["hmc"]["swap_update"]["nbonds"]
                burnin_swap_update = SwapUpdate(model,freq,nbonds)
            end
        end
    end

    return burnin_swap_update, sim_swap_update
end

"""
Initialize Simulation Stats Dictionary.
"""
function initialize_sim_stats()::Dict

    sim_stats = Dict("simulation_time" => 0.0, "measurement_time" => 0.0, "write_time" => 0.0,
                     "iters" => 0.0, "acceptance_rate" => 0.0, "reflect_acceptance_rate" => 0.0,
                     "swap_acceptance_rate" => 0.0)

    return sim_stats
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
@coeff function G(reactor::Reactor)
  nu_0 = 1 + 1.5 * nu_T

  _gamma(value) = Gamma(nu_0) * Gamma(value) / Gamma(nu_0+value)

  cur_value = 1/nu_0 - 2*sqrt(epsilon)*_gamma(5/4) + epsilon*_gamma(3/2)

  cur_value *= (1 + nu_T) ^ 1.5

  cur_value
end


text:
Base.@kwdef struct DictFunction{S,T,D<:AbstractDict{S,T}}
    dict::D
    default_value::T = zero(eltype(values(dict)))
end

function DictFunction(dict::AbstractDict)
    DictFunction(dict, zero(eltype(values(dict))))
end

function DictFunction(keys, values; kwargs...)
    DictFunction(dict = SortedDict{eltype(keys),eltype(values)}(keys .=> values); kwargs...)
end

function (f::DictFunction)(x)
    get(f.dict, x, f.default_value)
end

Base.keys(dictfun::DictFunction) = Base.keys(dictfun.dict)
Base.values(dictfun::DictFunction) = Base.values(dictfun.dict)

Base.@kwdef struct DiscretizedDictFunction{D<:Discretizer,F<:DictFunction}
    discretizer::D
    dictfunction::F
end

Base.keys(dictfun::DiscretizedDictFunction) = Base.keys(dictfun.dictfunction)

Base.values(dictfun::DiscretizedDictFunction) = Base.values(dictfun.dictfunction)

function DiscretizedDictFunction(discretizer, values; kwargs...)
    DiscretizedDictFunction(
        discretizer,
        DictFunction(Base.values(discretizer), values; kwargs...),
    )
end

function (f::DiscretizedDictFunction)(x)
    f.dictfunction(f.discretizer(x))
end

function dictfun(discretizer::Discretizer, Zs::AbstractVector{<:EBayesSample}, f)
    skedasticity(Zs) == Homoskedastic() ||
        error(ArgumentError("Heteroskedastic likelihood not implemented."))
    Z = Zs[1]
    dictfun(discretizer, Z, f )
end
# TODO: Better semantics for this lifting functionality.
function dictfun(discretizer::Discretizer, Z::EBayesSample, f)
    interval_Zs = [@set Z.Z = _int for _int in discretizer.sorted_intervals]
    if !isa(f, AbstractVector)
        f = f.(interval_Zs)
    end
    DiscretizedDictFunction(discretizer, DictFunction(interval_Zs, f))
end


text:
#!/usr/bin/env julia

using LightXML, CSV, DataFrames, TimeZones, Dates, ArgParse, LiveServer
using Interpolations: LinearInterpolation
using Base.Threads, FileWatching

push!(LOAD_PATH, @__DIR__)
using DiscgolfRecord

function read_gpx(gpx_file)
    xdoc = parse_file(gpx_file)
    xroot = root(xdoc)  # an instance of XMLElement
    trkpts = xroot["trk"][1]["trkseg"][1]["trkpt"]

    df = DataFrame(time=Dates.DateTime[], lat=Float64[], lon=Float64[])
    dfmt = DateFormat("y-m-dTH:M:S\\Z")
    # traverse all its child nodes and print element names
    for (i,t) in enumerate(trkpts)
        ad = attributes_dict(t)
        tim = Dates.DateTime(content(t["time"][1]), dfmt)
        push!(df, [tim, parse(Float64,ad["lat"]), parse(Float64,ad["lon"])])
    end

    free(xdoc)
    df
end

function read_timestamp(timestamp_file)
    df = DataFrame(time=Dates.DateTime[], disc=String[])
    f=CSV.File(timestamp_file, delim=',', datarow=1, header=false)
    dfmt = DateFormat("yyyy-mm-ddTHH:MM:SSzzzz")
    for (i,r) in enumerate(f)
        t = DateTime(ZonedDateTime(r.Column1, dfmt),UTC)
        disc = r.Column2
        push!(df,[t, disc])
    end
    df
end

"""
    interpolate_gpx(gpx_file, times) -> Vector{Location}

Interpolates a gpx file at the given timestamps and writes out csv
"""
function interpolate_gpx(timestamp_file, gpx_file, out)

    ts_df = read_timestamp(timestamp_file)

    gpx_df = read_gpx(gpx_file)

    lat_i = LinearInterpolation(datetime2unix.(gpx_df.time), gpx_df.lat)
    lon_i = LinearInterpolation(datetime2unix.(gpx_df.time), gpx_df.lon)

    df = DataFrame(time=DateTime[], lat=Float64[], lon=Float64[], disc=String[])

    for r in eachrow(ts_df)
        lat = lat_i(datetime2unix(r.time))
        lon = lon_i(datetime2unix(r.time))
        push!(df, [r.time, lat, lon, r.disc])
    end

    df

    CSV.write(out, df)

end

function write_html(loc)
    html = joinpath(@__DIR__,"guess_n_check_template.html")
    cp(html, joinpath(loc,"index.html"), force=true)
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--csv"
            help = "Timestamp file"
            default = "timestamp.csv"
        "--gpx"
            help = "Location recording file"
            default = "recording.gpx"
        "--out"
            help = "Name of output file"
            default = "round_raw.csv"
    end
    return parse_args(s)
end

args = parse_commandline()

interpolate_gpx(args["csv"], args["gpx"], args["out"])
println("Round has been parsed and written to $(args["out"])")

r = parse_round_raw_csv(args["out"])
round_id = id(r)
mkdir(round_id)

cp(args["csv"],joinpath(round_id,args["csv"]), force=true)
cp(args["gpx"],joinpath(round_id,args["gpx"]), force=true)
cp(args["out"],joinpath(round_id,args["out"]), force=true)

guess_and_check(joinpath(round_id, args["out"]), joinpath(round_id, "check_round.json"))

write_html(round_id)

Threads.@spawn begin
    while true
        watch_file(joinpath(round_id,args["out"]))
        println("Detected change in $(joinpath(round_id,args["out"]))")
        r = parse_round_raw_csv(joinpath(round_id,args["out"]))
        guess_and_check(joinpath(round_id, args["out"]), joinpath(round_id, "check_round.json"))
    end
end

serve(port=8080, dir=round_id)


text:
### convlayers
function layerBackProp(
    cLayer::ConvLayer,
    model::Model,
    actFun::SoS,
    labels::AbstractArray,
) where {SoS<:Union{Type{softmax},Type{σ}}}

    lossFun = model.lossFun
    dlossFun = Symbol("d", lossFun)
    A = cLayer.A
    Y = labels
    dZ = eval(:($dlossFun($A, $Y)))

    return dZ
end #softmax or σ layerBackProp

function layerBackProp!(
    cLayer::CL,
    model::Model,
    Ai::AbstractArray = Array{Any,1}(undef,0),
    Ao::AbstractArray = Array{Any,1}(undef,0),
    dAo::AbstractArray = Array{Any,1}(undef,0);
    labels::AbstractArray = Array{Any,1}(undef,0),
    kwargs...
) where {CL <: ConvLayer}

    NNlib = getindex(kwargs, :NNlib; default=true)
    img2col = getindex(kwargs, :img2col; default=true)

    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    if length(Ao) == 0
        Ao = cLayer.A
    end
    m = size(Ao)[end]

    regulization, λ = model.regulization, model.λ

    actFun = cLayer.actFun

    dZ = []
    if cLayer.actFun == model.outLayer.actFun
        dZ = layerBackProp(cLayer, model, eval(:($actFun)), labels)
    elseif length(dAo) != 0
        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dAo .* eval(:($dActFun($Z)))
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dAo = []
        for nextLayer in cLayer.nextLayers
            try
                dAo .+= nextLayer.dA
            catch e
                dAo = nextLayer.dA #need to initialize dA
            end #try/catch
        end #for

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dAo .* eval(:($dActFun($Z)))
    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    if NNlib

        dNNConv!(cLayer, dZ, Ai)

    else

        cLayer.dA = similar(Ai) #the size before padding
        cLayer.dA .= 0

        if img2col
            dimg2colConvolve!(cLayer, Ai, cLayer.dA, dZ)
        else
            dconvolve!(cLayer, Ai, cLayer.dA, dZ)
        end #if img2col

    end #if NNlib

    cLayer.backCount += 1

    return nothing
end #function layerBackProp!(cLayer::Input

### Pooling Layers

#import only the needed parts not to have conflict
import NNlib.∇maxpool, NNlib.∇meanpool, NNlib.∇maxpool!, NNlib.∇meanpool!, NNlib.PoolDims

function layerBackProp!(
    cLayer::OneD,
    model::Model,
    Ai::AbstractArray = Array{Any,1}(undef,0),
    Ao::AbstractArray = Array{Any,1}(undef,0),
    dAo::AbstractArray = Array{Any,1}(undef,0);
    labels::AbstractArray = Array{Any,1}(undef,0),
    kwargs...
) where {OneD<:Union{MaxPool1D,AveragePool1D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    if length(Ao) == 0
        Ao = cLayer.A
    end

    padS = paddingSize(cLayer, Ai)
    cLayer.dA = similar(Ai) .= 0

    if length(dAo) == 0
        if all(
            i -> (i.backCount == cLayer.nextLayers[1].backCount),
            cLayer.nextLayers,
        )
            dAo = []
            for nextLayer in cLayer.nextLayers
                try
                    dAo .+= nextLayer.dA
                catch e
                    dAo = nextLayer.dA #need to initialize dA
                end #try/catch
            end #for
        else
            return nothing
        end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)
    end #if dA==nothing

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            ∇maxpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            ∇meanpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer

    else

        if cLayer.s == cLayer.f
            dfastPooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        else
            dpooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        end #if cLayer.s == cLayer.f
    end #if NNlib

    cLayer.backCount += 1
    return nothing

end #unction layerBackProp!(cLayer::OneD) where {OneD <: Union{MaxPool1D, AveragePool1D}}

function layerBackProp!(
    cLayer::TwoD,
    model::Model,
    Ai::AbstractArray = Array{Any,1}(undef,0),
    Ao::AbstractArray = Array{Any,1}(undef,0),
    dAo::AbstractArray = Array{Any,1}(undef,0);
    labels::AbstractArray = Array{Any,1}(undef,0),
    kwargs...
) where {TwoD<:Union{MaxPool2D,AveragePool2D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    if length(Ao) == 0
        Ao = cLayer.A
    end

    padS = paddingSize(cLayer, Ai)
    cLayer.dA = similar(Ai) .= 0

    if length(dAo) == 0
        if all(
            i -> (i.backCount == cLayer.nextLayers[1].backCount),
            cLayer.nextLayers,
        )
            dAo = []
            for nextLayer in cLayer.nextLayers
                try
                    dAo .+= nextLayer.dA
                catch e
                    dAo = nextLayer.dA #need to initialize dA
                end #try/catch
            end #for
        else
            return nothing
        end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)
    end #if dA==nothing

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            ∇maxpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            ∇meanpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer

    else

        if cLayer.s == cLayer.f
            dfastPooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        else
            dpooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        end #if cLayer.s == cLayer.f
    end #if NNlib

    cLayer.backCount += 1
    return nothing

end #function layerBackProp!(cLayer::TwoD) where {TwoD <: Union{MaxPool2D, AveragePool2D}}

function layerBackProp!(
    cLayer::ThreeD,
    model::Model,
    Ai::AbstractArray = Array{Any,1}(undef,0),
    Ao::AbstractArray = Array{Any,1}(undef,0),
    dAo::AbstractArray = Array{Any,1}(undef,0);
    labels::AbstractArray = Array{Any,1}(undef,0),
    kwargs...,
) where {ThreeD<:Union{MaxPool3D,AveragePool3D}}

    NNlib = getindex(kwargs, :NNlib) ? kwargs[:NNlib] : true
    if length(Ai) == 0
        Ai = cLayer.prevLayer.A
    end
    if length(Ao) == 0
        Ao = cLayer.A
    end

    padS = paddingSize(cLayer, Ai)
    cLayer.dA = similar(Ai) .= 0

    if length(dAo) == 0
        if all(
            i -> (i.backCount == cLayer.nextLayers[1].backCount),
            cLayer.nextLayers,
        )
            dAo = []
            for nextLayer in cLayer.nextLayers
                try
                    dAo .+= nextLayer.dA
                catch e
                    dAo = nextLayer.dA #need to initialize dA
                end #try/catch
            end #for
        else
            return nothing
        end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)
    end #if dA==nothing

    if NNlib
        pooldims = PoolDims(Ai, cLayer.f, stride = cLayer.s, padding = padS)
        if cLayer isa MaxPoolLayer
            ∇maxpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        elseif cLayer isa AveragePoolLayer
            ∇meanpool!(cLayer.dA, dAo, Ao, Ai, pooldims)
        end #if cLayer isa MaxPoolLayer

    else

        if cLayer.s == cLayer.f
            dfastPooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        else
            dpooling!(cLayer, Ai, cLayer.dA, Ao, dAo)
        end #if cLayer.s == cLayer.f
    end #if NNlib

    cLayer.backCount += 1
    return nothing

end #function layerBackProp!(cLayer::ThreeD) where {ThreeD <: Union{MaxPool3D, AveragePool3D}}


text:
##########################
### Points Conversions ###
##########################

"""
Create ECEF coordinates from a given `lla`
"""
function ECEF(lla::LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg, h = lla.lat, lla.lon, lla.alt
    d = datum

    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
    sinλ, cosλ = sind(λdeg), cosd(λdeg)

    N = d.a / sqrt(1 - d.e² * sinϕ^2)  # Radius of curvature (meters)

    x = (N + h) * cosϕ * cosλ
    y = (N + h) * cosϕ * sinλ
    z = (N * (1 - d.e²) + h) * sinϕ

    return ECEF(x, y, z)
end
"""
Create LLA coordinates from a given `ecef`
"""
function LLA(ecef::ECEF, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    x, y, z = ecef.x, ecef.y, ecef.z
    d = datum

    p = hypot(x, y)
    θ = atan(z*d.a, p*d.b)
    λ = atan(y, x)
    ϕ = atan(z + d.e′² * d.b * sin(θ)^3, p - d.e²*d.a*cos(θ)^3)

    N = d.a / sqrt(1 - d.e² * sin(ϕ)^2)  # Radius of curvature (meters)
    h = p / cos(ϕ) - N

    return LLA(rad2deg(ϕ), rad2deg(λ), h)
end

"""
Create ENU coordinates from given `ecef` coordinates
given a reference point `lla_ref` for linarization
"""
function ENU(ecef::ECEF, lla_ref::LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg = lla_ref.lat, lla_ref.lon

    ecef_ref = ECEF(lla_ref, datum)
    ∂x = ecef.x - ecef_ref.x
    ∂y = ecef.y - ecef_ref.y
    ∂z = ecef.z - ecef_ref.z

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [     -sinλ       cosλ  0.0
    #      -cosλ*sinϕ -sinλ*sinϕ cosϕ
    #       cosλ*cosϕ  sinλ*cosϕ sinϕ]
    #
    # east, north, up = R * [∂x, ∂y, ∂z]
    east  = ∂x * -sinλ      + ∂y * cosλ       + ∂z * 0.0
    north = ∂x * -cosλ*sinϕ + ∂y * -sinλ*sinϕ + ∂z * cosϕ
    up    = ∂x * cosλ*cosϕ  + ∂y * sinλ*cosϕ  + ∂z * sinϕ

    return ENU(east, north, up)
end

"""
Create ENU coordinates from given `ecef` coordinates
given a center of reference point `bounds` for linarization
"""
ENU(ecef::ECEF, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = ENU(ecef, OpenStreetMapX.center(bounds), datum)

"""
Create ECEF coordinates from given `enu` coordinates
and a reference point being center of `bounds` for linearization
"""
function ECEF(enu::ENU, lla_ref::LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg = lla_ref.lat, lla_ref.lon

    ecef_ref = ECEF(lla_ref, datum)

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [-sinλ -sinϕ*cosλ  cosϕ*cosλ
    #      	cosλ -sinϕ*sinλ  cosϕ*sinλ
    #        0.0       cosϕ       sinϕ]
    #
    # x,y,z = R * [east, north, up] + [ecef_ref.x, ecef_ref.y, ecef_ref.z]

    x  = -sinλ*enu.east + -sinϕ*cosλ*enu.north + cosϕ*cosλ*enu.up + ecef_ref.x
    y = cosλ*enu.east +  -sinϕ*sinλ*enu.north + cosϕ*sinλ*enu.up + ecef_ref.y
    z    = 0.0*enu.east + cosϕ*enu.north + sinϕ*enu.up + ecef_ref.z

    return ECEF(x, y, z)
end

"""
Create ECEF coordinates from given `enu` coordinates
and a reference point being center of `bounds` for linearization
"""
ECEF(enu::ENU, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = ECEF(enu, OpenStreetMapX.center(bounds), datum)

"""
Create ENU coordinates from given `lla` coordinates
and a reference point `lla_ref` for linearization
"""
ENU(lla::LLA, lla_ref::LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = ENU(ECEF(lla, datum), lla_ref, datum)

"""
Create ENU coordinates from given `lla` coordinates
and a reference point being center of `bounds` for linearization
"""
ENU(lla::LLA, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = ENU(ECEF(lla, datum), bounds, datum)

"""
Create LLA coordinates from given `enu` coordinates
and a reference point `lla_ref` for linearization
"""
LLA(enu::ENU, lla_ref::LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = LLA(ECEF(enu,lla_ref))

"""
Create LLA coordinates from given `enu` coordinates
and a reference point being center of `bounds` for linearization
"""
LLA(enu::ENU, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = LLA(ECEF(enu,bounds))

"""
Converts a dictionary of `LLA` `nodes` into a dictionary of `ECEF` values.
Uses a reference point `lla_ref` for linearization.
"""
function ECEF(nodes::Dict{Int,LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,ECEF}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = ECEF(node, datum)
    end

    return r
end

"""
Converts a dictionary of `ECEF` `nodes` into a dictionary of `LLA` values.
"""
function LLA(nodes::Dict{Int,ECEF}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,LLA}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = LLA(node, datum)
    end

    return r
end

"""
Converts a dictionary of `LLA` and `ECEF` `nodes` into a dictionary of `ENU` values.
Uses a reference point `lla_ref` for linearization.
"""
function ENU(nodes::Dict{Int,T}, lla_ref::LLA,
            datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) where T<:Union{LLA,ECEF}
    r = Dict{Int,ENU}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = ENU(node, lla_ref, datum)
    end

    return r
end

"""
Converts a dictionary of `LLA` and `ECEF` `nodes` into a dictionary of `ENU` values.
Uses the center of the given `bounds` for linearization.
"""
ENU(nodes::Dict{Int,T}, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) where T<:Union{LLA,ECEF} = ENU(nodes, OpenStreetMapX.center(bounds), datum)

"""
Converts a dictionary of `ENU` `nodes` into a dictionary of `ECEF` values.
Uses a reference point `lla_ref` for linearization.
"""
function ECEF(nodes::Dict{Int,ENU},lla_ref::LLA , datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,ECEF}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = ECEF(node, lla_ref, datum)
    end

    return r
end

"""
Converts a dictionary of `ENU` `nodes` into a dictionary of `ECEF` values.
Uses the center of the given `bounds` for linearization.
"""
ECEF(nodes::Dict{Int,ENU}, bounds::Bounds{LLA}, datum::Ellipsoid = WGS84) = ECEF(nodes, OpenStreetMapX.center(bounds), datum)

"""
Converts a dictionary of `ENU` `nodes` into a dictionary of `LLA` values.
Uses a reference point `lla_ref` for linearization.
"""
function LLA(nodes::Dict{Int,ENU},lla_ref::LLA , datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,LLA}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = LLA(node, lla_ref, datum)
    end

    return r
end

"""
Converts a dictionary of `ENU` `nodes` into a dictionary of `LLA` values.
Uses the center of the given `bounds` for linearization.
"""
LLA(nodes::Dict{Int,ENU}, bounds::Bounds{LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = LLA(nodes, OpenStreetMapX.center(bounds), datum)


text:
using Documenter, NarmaxLasso

makedocs(
    format = :html,
    sitename = "NarmaxLasso.jl",
    authors="Antonio H. Ribeiro",
    pages = Any[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Guide" => "guide.md",
        "Library" => "library.md"
    ]
)

deploydocs(
    repo = "github.com/antonior92/NarmaxLasso.jl.git",
    target = "build",
    julia = "0.6",
    osname = "linux",
    deps = nothing,
    make = nothing,
)


text:
import Base: eltype

"""
All networks in the package belong to the `AbstractEcologicalNetwork` type. They
all have a field `A` to represent interactions as a *matrix*, and a number of
fields for species. See the documentation for `AbstractBipartiteNetwork` and
`AbstractUnipartiteNetwork`, as well as `AllowedSpeciesTypes` for the allowed
types for species.

Note that *all* species in a network (including both levels of a bipartite
network) *must* have the same type. For example, `["a", :b, "c"]` is not a valid
array of species, as not all its elements have the same type.
"""
abstract type AbstractEcologicalNetwork end

"""
This abstract type groups all unipartite networks, regardless of the type of
information. Unipartite networks have *a single* field for species, named `S`,
which has the same number of elements as the size of the matrix.

Any unipartite network can be declared (we'll use the example of a binary
network) either using `UnipartiteNetwork(A, S)` (assuming `A` is a matrix of
interactions and `S` is a vector of species names), or `UnipartiteNetwork(A)`,
in which case the species will be named automatically.
"""
abstract type AbstractUnipartiteNetwork <: AbstractEcologicalNetwork end

"""
This abstract type groups all bipartite networks, regardless of the type of
information. Bipartite networks have *two* fields for species, named `T` (for
top, corresponding to matrix *rows*), and `B` (for bottom, matrix *columns*).

Any bipartite network can be declared (we'll use the example of a binary
network) either using `BipartiteNetwork(A, T, B)` (assuming `A` is a matrix of
interactions and `T` and `B` are vectors of species names for the top and bottom
level), or `BipartiteNetwork(A)`, in which case the species will be named
automatically.
"""
abstract type AbstractBipartiteNetwork <: AbstractEcologicalNetwork end

"""
A bipartite deterministic network is a matrix of boolean values.
"""
mutable struct BipartiteNetwork{Bool, ST} <: AbstractBipartiteNetwork
    edges::SparseMatrixCSC{Bool,Int64}
    T::Vector{ST}
    B::Vector{ST}
    function BipartiteNetwork{Bool, NT}(edges::M, T::Vector{NT}, B::Vector{NT}) where {M<:SparseMatrixCSC, NT}
        dropzeros!(edges)
        check_bipartiteness(edges, T, B)
        new{Bool,NT}(edges, T, B)
    end
end

function BipartiteNetwork(A::M, T::Union{Vector{TT},Nothing}=nothing, B::Union{Vector{TT},Nothing}=nothing) where {M <: AbstractMatrix{Bool}, TT}
    if isnothing(B)
        B = "b".*string.(1:size(A, 2))
    else
        _check_species_validity(TT)
    end
    if isnothing(T)
        T = "t".*string.(1:size(A, 1))
    else
        _check_species_validity(TT)
    end
    allunique(T) || throw(ArgumentError("All top-level species must be unique"))
    allunique(B) || throw(ArgumentError("All bottom-level species must be unique"))
    allunique(vcat(B,T)) || throw(ArgumentError("Bipartite networks cannot share species across levels"))
    isequal(length(T))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(B))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return BipartiteNetwork{Bool,eltype(T)}(sparse(A), T, B)
end

"""
An unipartite deterministic network is a matrix of boolean values.
"""
mutable struct UnipartiteNetwork{Bool, ST} <: AbstractUnipartiteNetwork
    edges::SparseMatrixCSC{Bool,Int64}
    S::Vector{ST}
    function UnipartiteNetwork{Bool, ST}(edges::M, S::Vector{ST}) where {M<:SparseMatrixCSC, ST}
        check_unipartiteness(edges, S)
        dropzeros!(edges)
        new{Bool,ST}(edges, S)
    end
end

function UnipartiteNetwork(A::M, S::Union{Vector{TT},Nothing}=nothing) where {M <: AbstractMatrix{Bool}, TT}
    if isnothing(S)
        S = "s".*string.(1:size(A, 1))
    else
        _check_species_validity(TT)
    end
    allunique(S) || throw(ArgumentError("All species must be unique"))
    isequal(length(S))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(S))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return UnipartiteNetwork{Bool,eltype(S)}(sparse(A), S)
end

"""
A bipartite probabilistic network is a matrix of floating point numbers, all of
which must be between 0 and 1.
"""
mutable struct BipartiteProbabilisticNetwork{IT <: AbstractFloat, ST} <: AbstractBipartiteNetwork
    edges::SparseMatrixCSC{IT}
    T::Vector{ST}
    B::Vector{ST}
    function BipartiteProbabilisticNetwork(edges, T, B)
        dropzeros!(edges)
        allunique(T) || throw(ArgumentError("All top-level species must be unique"))
        allunique(B) || throw(ArgumentError("All bottom-level species must be unique"))
        allunique(vcat(B,T)) || throw(ArgumentError("Bipartite networks cannot share species across levels"))
        check_probability_values(edges)
        _check_species_validity(eltype(T))
        _check_species_validity(eltype(B))
        check_bipartiteness(edges, T, B)
        new{eltype(edges),eltype(T)}(edges, T, B)
    end
    function BipartiteProbabilisticNetwork{IT,ST}(edges::SparseMatrixCSC{IT}, T::Vector{ST}, B::Vector{ST}) where {IT <: Real, ST}
        return BipartiteProbabilisticNetwork(edges, T, B)
    end
end

function BipartiteProbabilisticNetwork(A::Array{IT,2}, T::Union{Vector{TT},Nothing}=nothing, B::Union{Vector{TT},Nothing}=nothing) where {IT <: AbstractFloat, TT}
    isnothing(B) ? (B = "b".*string.(1:size(A, 2))) : _check_species_validity(TT)
    isnothing(T) ? (T = "t".*string.(1:size(A, 1))) : _check_species_validity(TT)
    isequal(length(T))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(B))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return BipartiteProbabilisticNetwork(sparse(A), T, B)
end

"""
A bipartite quantitative network is matrix of numbers. It is assumed that the
interaction strength are *positive*.
"""
mutable struct BipartiteQuantitativeNetwork{IT <: Number, ST} <: AbstractBipartiteNetwork
    edges::SparseMatrixCSC{IT}
    T::Vector{ST}
    B::Vector{ST}
    function BipartiteQuantitativeNetwork(edges, T, B)
        dropzeros!(edges)
        allunique(T) || throw(ArgumentError("All top-level species must be unique"))
        allunique(B) || throw(ArgumentError("All bottom-level species must be unique"))
        allunique(vcat(B,T)) || throw(ArgumentError("Bipartite networks cannot share species across levels"))
        _check_species_validity(eltype(T))
        _check_species_validity(eltype(B))
        check_bipartiteness(edges, T, B)
        new{eltype(edges),eltype(T)}(edges, T, B)
    end
    function BipartiteQuantitativeNetwork{IT,ST}(edges::SparseMatrixCSC{IT}, T::Vector{ST}, B::Vector{ST}) where {IT <: Real, ST}
        return BipartiteQuantitativeNetwork(edges, T, B)
    end
end

function BipartiteQuantitativeNetwork(A::Array{IT,2}, T::Union{Vector{TT},Nothing}=nothing, B::Union{Vector{TT},Nothing}=nothing) where {IT <: Number, TT}
    isnothing(B) ? (B = "b".*string.(1:size(A, 2))) : _check_species_validity(TT)
    isnothing(T) ? (T = "t".*string.(1:size(A, 1))) : _check_species_validity(TT)
    isequal(length(T))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(B))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return BipartiteQuantitativeNetwork(sparse(A), T, B)
end

"""
A unipartite probabilistic network is a square matrix of floating point numbers,
all of which must be between 0 and 1.
"""
mutable struct UnipartiteProbabilisticNetwork{IT <: AbstractFloat, ST} <: AbstractUnipartiteNetwork
    edges::SparseMatrixCSC{IT}
    S::Vector{ST}
    function UnipartiteProbabilisticNetwork(edges, S)
        dropzeros!(edges)
        allunique(S) || throw(ArgumentError("All species must be unique"))
        check_probability_values(edges)
        _check_species_validity(eltype(S))
        check_unipartiteness(edges, S)
        new{eltype(edges),eltype(S)}(edges, S)
    end
    function UnipartiteProbabilisticNetwork{IT,ST}(edges::SparseMatrixCSC{IT}, S::Vector{ST}) where {IT <: Real, ST}
        return UnipartiteProbabilisticNetwork(edges, S)
    end
end

function UnipartiteProbabilisticNetwork(A::Array{IT,2}, S::Union{Vector{TT},Nothing}=nothing) where {IT <: AbstractFloat, TT}
    isnothing(S) ? (S = "s".*string.(1:size(A, 2))) : _check_species_validity(TT)
    isequal(length(S))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(S))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return UnipartiteProbabilisticNetwork(sparse(A), S)
end

"""
A unipartite quantitative network is a square matrix of numbers.
"""
mutable struct UnipartiteQuantitativeNetwork{IT <: Number, ST} <: AbstractUnipartiteNetwork
    edges::SparseMatrixCSC{IT}
    S::Vector{ST}
    function UnipartiteQuantitativeNetwork(edges, S)
        dropzeros!(edges)
        allunique(S) || throw(ArgumentError("All species must be unique"))
        _check_species_validity(eltype(S))
        check_unipartiteness(edges, S)
        new{eltype(edges),eltype(S)}(edges, S)
    end
    function UnipartiteQuantitativeNetwork{IT,ST}(edges::SparseMatrixCSC{IT}, S::Vector{ST}) where {IT <: Number, ST}
        return UnipartiteQuantitativeNetwork(edges, S)
    end
end

function UnipartiteQuantitativeNetwork(A::Array{IT,2}, S::Union{Vector{TT},Nothing}=nothing) where {IT <: Number, TT}
    isnothing(S) ? (S = "s".*string.(1:size(A, 2))) : _check_species_validity(TT)
    isequal(length(S))(size(A,1)) || throw(ArgumentError("The matrix has the wrong number of top-level species"))
    isequal(length(S))(size(A,2)) || throw(ArgumentError("The matrix has the wrong number of bottom-level species"))
    return UnipartiteQuantitativeNetwork(sparse(A), S)
end

"""
This is a union type for both Bipartite and Unipartite probabilistic networks.
Probabilistic networks are represented as arrays of floating point values ∈
[0;1].
"""
ProbabilisticNetwork = Union{BipartiteProbabilisticNetwork, UnipartiteProbabilisticNetwork}

"""
This is a union type for both Bipartite and Unipartite deterministic networks.
All networks from these class have adjacency matrices represented as arrays of
Boolean values.
"""
BinaryNetwork = Union{BipartiteNetwork, UnipartiteNetwork}

"""
This is a union type for both unipartite and bipartite quantitative networks.
All networks of this type have adjancency matrices as two-dimensional arrays of
numbers.
"""
QuantitativeNetwork = Union{BipartiteQuantitativeNetwork, UnipartiteQuantitativeNetwork}

"""
All non-probabilistic networks
"""
DeterministicNetwork = Union{BinaryNetwork, QuantitativeNetwork}


text:
using Distributions
using Turing

x = Float64[1 2]
priors = 0

@model gauss(x) = begin
  priors = TArray{Float64}(2)
  priors[1] ~ InverseGamma(2,3)         # s
  priors[2] ~ Normal(0, sqrt(priors[1])) # m
  for i in 1:length(x)
    x[i] ~ Normal(priors[2], sqrt(priors[1]))
  end
  priors
end

chain = sample(gauss(x), PG(10, 10))
chain = sample(gauss(x), SMC(10))