text:
using Logging
using DataFrames
using ZipFile
using CSV
using Serialization
using JLD2
using Feather
using Arrow
using Parquet
using JDF

export file_format, compression_formats, compress, uncompress, df_formats, df_write, df_read

"""
    file_format(file::AbstractString)::Symbol

Returns the file format from the file suffix.

# Arguments
- `file::AbstractString`: file path
"""
file_format(file::AbstractString)::Symbol = Symbol(splitext(file)[2][2:end])

"""
    compression_formats()::Vector{Symbol}

Returns the file formats supported for compressed archives.
"""
compression_formats()::Vector{Symbol} = [:zip]

"""
    compress(f::Function, file::AbstractString, subfile::AbstractString)

Creates a compressed file.  The type of compression is inferred from the
output file name.  subfile is the name of the file within the compressed
archive.  f is a function which writes to the archive.

The typical use case is:

    compress(file, subfile) do io
        # put here the body of function f(io) to write to io
    end

# Arguments
- `f::Function`: function to execute on the file contents
- `file::AbstractString`: compressed file path
- `subfile::AbstractString`: subfile within the archive
"""
function compress(f::Function, file::AbstractString, subfile::AbstractString)
    format = file_format(file)
    @debug "BEGIN compress: $format $file"
    t = time()

    if format == :zip
        local zfiles = ZipFile.Writer(file)

        try
            local zf = ZipFile.addfile(zfiles, subfile)

            try
                f(zf)
            finally
                close(zf)
            end
        finally
            close(zfiles)
        end
    else
        throw(ErrorException("Unsupported compression format: format=$format"))
    end

    @debug "END compress: $format $(time() - t)"
end

"""
    uncompress(f::Function, file::AbstractString)

Uncompresses a compressed file.  The type of compression is inferred from the
output file name.  f is a function which reads from the archive.

The typical use case is:

    uncompress(file) do io
        # put here the body of function f(io) to write to io
    end

# Arguments
- `f::Function`: function to execute on the file contents
- `file::AbstractString`: file path
"""
function uncompress(f::Function, file::AbstractString)
    format = file_format(file)
    @debug "BEGIN uncompress: $format $file"
    t = time()

    if format == :zip

        zfiles = ZipFile.Reader(file)

        try
            if size(zfiles.files,1) != 1
                throw(ErrorException("Zip file does not contain exactly one file: zipfile=$file contents=$(zfiles.files)"))
            end

            zf = zfiles.files[1]

            try
                f(zf.name, zf)
            finally
                close(zf)
            end
        finally
            close(zfiles)
        end
    else
        throw(ErrorException("Unsupported compression format: format=$format"))
    end

    @debug "END uncompress: $format $(time() - t)"
end

"""
    df_formats()::Vector{Symbol}

Returns the file formats supported for DataFrames.
"""
df_formats()::Vector{Symbol} = [:csv, :ser, :jld2, :jld2c, :feather, :arrow, :arrow_lz4, :arrow_zstd, :parquet, :jdf]

"""
    df_write(file::AbstractString, df::DataFrame; subformat::Union{Nothing,Symbol}=nothing, dictencode::Bool=true)

Writes a DataFrame to a file.  The file suffix determines how the DataFrame is serialized.
If the file has a compressed suffix, subformat determines how the DataFrame
is serialized in the compressed archive.

# Arguments
- `file::AbstractString`: file path
- `df::DataFrame`: dataframe to export
- `subformat::Union{Nothing,Symbol}=nothing`: format to used within a compressed archive.
- `dictencode::Bool=true`: true to use a dictionary encoding, if possible, and false otherwise.
"""
function df_write(file::AbstractString, df::DataFrame; subformat::Union{Nothing,Symbol}=nothing, dictencode::Bool=true)
    format = file_format(file)

    if in(format, compression_formats())
        if subformat == nothing
            throw(ErrorException("Compressed archives must specify the DataFrame serialization format (subformat)."))
        end

        subfile = "db." * String(subformat)

        compress(file, subfile) do f
            _df_write(f, subformat, df, dictencode=dictencode)
        end

        return
    end

    _df_write(file, format, df, dictencode=dictencode)
end

"""
Writes a DataFrame.  File can be a file path or an IO.
- `dictencode::Bool=true`: true to use a dictionary encoding, if possible, and false otherwise.
"""
function _df_write(file, format::Symbol, df::DataFrame; dictencode::Bool=true)

    @debug "BEGIN df_write: $format $file"
    t = time()

    if format == :csv
        CSV.write(file, df)
    elseif format == :ser
        serialize(file, df)
    elseif format == :jld2
        if !isa(file, AbstractString)
            throw(ErrorException("JLD2 is not supported in compressed archives."))
        end

        JLD2.@save file df=df
    elseif format == :jld2c
        if !isa(file, AbstractString)
            throw(ErrorException("JLD2 is not supported in compressed archives."))
        end

        JLD2.@save file {compress=true} df=df
    elseif format == :feather
        Feather.write(file, df)
    elseif format == :arrow
        Arrow.write(file, df; dictencode=true)
    elseif format == :arrow_lz4
        Arrow.write(file, df; compress=:lz4, dictencode=dictencode)
    elseif format == :arrow_zstd
        Arrow.write(file, df; compress=:zstd, dictencode=dictencode)
    elseif format == :parquet
        if !isa(file, AbstractString)
            throw(ErrorException("Parquet is not supported in compressed archives."))
        end

        Parquet.write_parquet(file, df)
    elseif format == :jdf
        if !isa(file, AbstractString)
            throw(ErrorException("JDF is not supported in compressed archives."))
        end

        JDF.savejdf(file, df)
    else
        throw(ErrorException("Unsupported dataframe format: format=$format"))
    end

    @debug "END df_write: $format $(time() - t)"
end

"""
    df_read(file::AbstractString; dates_as_strings::Bool=true, missing_type::Type=String, missing_types::Dict{String,Type}=Dict{String,Type}())::DataFrame

Reads a DataFrame from a file.  The file suffix determines how the DataFrame is deserialized.

# Arguments
- `file::AbstractString`: file path
- `dates_as_strings::Bool=true`: true to parse dates as strings; false to parse dates as dates.
- `missing_type::Type=String`: type to use if all values are missing,
- `missing_types::Dict{String,Type}=Dict{String,Type}()`: type to use for missing columns
"""
function df_read(file::AbstractString; dates_as_strings::Bool=true, missing_type::Type=String, missing_types::Dict{String,Type}=Dict{String,Type}())::DataFrame
    format = file_format(file)

    if in(format, compression_formats())
        local df

        uncompress(file) do name, f
            format = file_format(name)
            df = _df_read(f, format; dates_as_strings=dates_as_strings, missing_type=missing_type, missing_types=missing_types)
        end

        return df
    end

    return _df_read(file, format; dates_as_strings=dates_as_strings, missing_type=missing_type, missing_types=missing_types)
end

"""
Reads a DataFrame.  File can be a file path or an IO.
"""
function _df_read(file, format::Symbol; dates_as_strings::Bool=true, missing_type::Type=String, missing_types::Dict{String,Type}=Dict{String,Type}())::DataFrame
    @debug "BEGIN df_read $format: $file"
    t = time()

    if format == :csv
        if dates_as_strings
            df = CSV.File(read(file); dateformat="---DON'T PARSE DATES---") |> DataFrame
        else
            df = CSV.read(file, DataFrame)
        end
    elseif format == :ser
        df = deserialize(file)
    elseif format == :jld2 || format == :jld2c
        if !isa(file, AbstractString)
            throw(ErrorException("JLD2 is not supported in compressed archives."))
        end

        JLD2.@load file df
    elseif format == :feather
        df = Feather.read(file)
    elseif format == :arrow  || format == :arrow_lz4 || format == :arrow_zstd
        df = DataFrame(Arrow.Table(file))
    elseif format == :parquet
        if !isa(file, AbstractString)
            throw(ErrorException("Parquet is not supported in compressed archives."))
        end

        parquetfile = Parquet.File(file)

        try
            local cc = BatchedColumnsCursor(parquetfile)
            local batchvals, state = iterate(cc)

            df = DataFrame()

            for key in keys(batchvals)
                df[!,key] = batchvals[key]
            end

        finally
            close(parquetfile)
        end

    elseif format == :jdf
        if !isa(file, AbstractString)
            throw(ErrorException("JDF is not supported in compressed archives."))
        end

        df = DataFrame(JDF.loadjdf(file));
    else
        throw(ErrorException("Unsupported dataframe format: format=$format"))
    end

    # Some serialization methods can not handle columns of type Missing.
    for n in names(df)
        if eltype(df[!,n]) == Missing
            T = haskey(missing_types, n) ? missing_types[n] : missing_type
            @debug "Column with all missing values retyped as Vector{Union{Missing, $T}}: name=$n"
            df[!,n] = Vector{Union{Missing, T}}(missing, size(df, 1))
        end
    end

    @debug "END df_read $format: $file $(time() - t)"
    return df
end


text:
#The m-by-n Hilbert matrix has matrix elements
# H_{ij} = 1/(i+j-1)
export Hilbert

struct Hilbert{T}
    m :: Int
    n :: Int
end
Hilbert(T::Type, m::Integer, n::Integer) = Hilbert{T}(m, n)
Hilbert(m::Integer, n::Integer) = Hilbert{Rational{Int}}(m, n)
Hilbert(n::Integer) = Hilbert(n, n)

function Matrix(H::Hilbert{T}) where T
    Hf = zeros(T, H.m, H.n)
    for i=1:H.m, j=1:H.n
        Hf[i, j] = one(T)/(i+j-1)
    end
    Hf
end

function inv(H::Hilbert{T}) where T
    invH=zeros(T,H.m,H.m)
    for i=1:H.m, j=1:H.m
        invH[i,j]=(-1)^(i+j)*(i+j-1)*binomial(H.m+i-1,H.m-j)*
                  binomial(H.m+j-1,H.m-i)*binomial(i+j-2,i-1)^2
    end
    invH
end

ext: jl


text:
export
    b1900,
    b1950,
    badkpv,
    bltfrm,
    bodc2n,
    bodc2s,
    boddef,
    bodfnd,
    bodn2c,
    bods2c,
    bodvcd,
    bodvrd,
    brcktd,
    brckti,
    bschoc,
    bschoi,
    bsrchc,
    bsrchd,
    bsrchi

"""
    b1900()

Returns the Julian Date corresponding to Besselian date 1900.0.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/b1900_c.html)
"""
function b1900()
    ccall((:b1900_c, libcspice), SpiceDouble, ())
end

"""
    b1950()

Returns the Julian Date corresponding to Besselian date 1950.0.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/b1950_c.html)
"""
function b1950()
    ccall((:b1950_c, libcspice), SpiceDouble, ())
end

"""
    badkpv(caller, name, comp, size, divby, typ)

Determine if a kernel pool variable is present and if so that it has the correct size and type.

### Arguments ###

- `caller`: Name of the routine calling this routine
- `name`: Name of a kernel pool variable
- `comp`: Comparison operator
- `size`: Expected size of the kernel pool variable
- `divby`: A divisor of the size of the kernel pool variable
- `type`: Expected type of the kernel pool variable

### Output ###

The function returns `false` if the kernel pool variable is OK otherwise an exception is thrown.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/badkpv_c.html)
"""
function badkpv(caller, name, comp, size, divby, typ)
    res = ccall((:badkpv_c, libcspice), SpiceBoolean,
                (Cstring, Cstring, Cstring, SpiceInt, SpiceInt, SpiceChar),
                caller, name, comp, size, divby, first(typ))
    handleerror()
    Bool(res)
end

"""
    bltfrm(frmcls)

Return a SPICE set containing the frame IDs of all built-in frames of a specified class.

### Arguments ###

- `frmcls`: Frame class

### Output ###

- `idset`: Set of ID codes of frames of the specified class

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bltfrm_c.html)
"""
function bltfrm(frmcls)
    idset = SpiceIntCell(256)
    data = idset.data
    GC.@preserve data begin
        ccall((:bltfrm_c, libcspice), Cvoid, (SpiceInt, Ref{Cell{SpiceInt}}), frmcls, idset.cell)
    end
    handleerror()
    idset
end

"""
    bodc2n(code)

Translate the SPICE integer code of a body into a common name for that body.

### Arguments ###

- `code`: Integer ID code to be translated into a name

### Output ###

A common name for the body identified by code or `nothing` if none was found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodc2n_c.html)
"""
function bodc2n(code)
    len = 36
    name = Array{SpiceChar}(undef, len)
    found = Ref{SpiceBoolean}()
    ccall((:bodc2n_c, libcspice), Cvoid, (SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceBoolean}),
          code, len, name, found)
    Bool(found[]) || return nothing
    chararray_to_string(name)
end

"""
    bodc2s(code)

Translate a body ID code to either the corresponding name or if no name to ID code mapping exists,
the string representation of the body ID value.

### Arguments ###

- `code`: Integer ID code to translate to a string

### Output ###

Returns a string corresponding to `code`

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodc2n_c.html)
"""
function bodc2s(code)
    len = 36
    name = Array{SpiceChar}(undef, len)
    ccall((:bodc2s_c, libcspice), Cvoid, (SpiceInt, SpiceInt, Ref{SpiceChar}),
          code, len, name)
    chararray_to_string(name)
end

"""
    boddef(name, code)

Define a body name/ID code pair for later translation via [`bodn2c`](@ref) or [`bodc2n`](@ref).

### Arguments ###

- `name`: Common name of some body
- `code`: Integer code for that body

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/boddef_c.html)
"""
function boddef(name, code)
    length(name) > 35 && throw(SpiceError("The maximum allowed length for a name is 35 characters."))
    ccall((:boddef_c, libcspice), Cvoid, (Cstring, SpiceInt), name, code)
    handleerror()
end

"""
    bodfnd(body, item)

Determine whether values exist for some item for any body in the kernel pool.

### Arguments ###

- `body`: ID code of body
- `item`: Item to find ("RADII", "NUT_AMP_RA", etc.)

### Output ###

Returns `true` if the item is in the kernel pool and `false` if it is not.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodfnd_c.html)
"""
function bodfnd(body, item)
    Bool(ccall((:bodfnd_c, libcspice), SpiceBoolean, (SpiceInt, Cstring), body, item))
end

"""
    bodn2c(name)

Translate the name of a body or object to the corresponding SPICE integer ID code.

### Arguments ###

- `name`: Body name to be translated into a SPICE ID code

### Output ###

Return the SPICE integer ID code for the named body or `nothing` if none was found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodn2c_c.html)
"""
function bodn2c(name)
    code = Ref{SpiceInt}()
    found = Ref{SpiceBoolean}()
    ccall((:bodn2c_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}, Ref{SpiceBoolean}),
          name, code, found)
    Bool(found[]) || return nothing
    Int(code[])
end

"""
    bods2c(name)

Translate a string containing a body name or ID code to an integer code.

### Arguments ###

- `name`: String to be translated to an ID code

### Output ###

Retunrs the integer ID code corresponding to `name` or `nothing` if none as found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bods2c_c.html)
"""
function bods2c(name)
    code = Ref{SpiceInt}()
    found = Ref{SpiceBoolean}()
    ccall((:bods2c_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}, Ref{SpiceBoolean}),
          name, code, found)
    Bool(found[]) || return nothing
    Int(code[])
end

"""
    bodvcd(bodyid, item)

Fetch from the kernel pool the double precision values of an item associated with a body, where the body
is specified by an integer ID code.

### Arguments ###

- `bodyid`: Body ID code
- `item`: Item for which values are desired. (`"RADII"`, `"NUT_PREC_ANGLES"`, etc.)
- `maxn`: Maximum number of values that may be returned (default: 100)

### Output ###

Returns the requested values.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodvcd_c.html)
"""
function bodvcd(bodyid, item, maxn=100)
    values = Array{SpiceDouble}(undef, maxn)
    dim = Ref{SpiceInt}()
    ccall((:bodvcd_c, libcspice), Cvoid, (SpiceInt, Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceDouble}),
          bodyid, item, maxn, dim, values)
    handleerror()
    values[1:dim[]]
end

"""
    bodvrd(bodynm, item)

Fetch from the kernel pool the double precision values of an item associated with a body.

### Arguments ###

- `bodynm`: Body name
- `item`: Item for which values are desired. (`"RADII"`, `"NUT_PREC_ANGLES"`, etc.)
- `maxn`: Maximum number of values that may be returned (default: 100)

### Output ###

- `values`: Values

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/bodvrd_c.html)
"""
function bodvrd(bodynm, item, maxn=100)
    values = Array{SpiceDouble}(undef, maxn)
    dim = Ref{SpiceInt}()
    ccall((:bodvrd_c, libcspice), Cvoid, (Cstring, Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceDouble}),
          bodynm, item, maxn, dim, values)
    handleerror()
    values[1:dim[]]
end

function _brcktd(number, e1, e2)
    ccall((:brcktd_c, libcspice), SpiceDouble, (SpiceDouble, SpiceDouble, SpiceDouble),
        number, e1, e2)
end

"""
    brcktd(number, e1, e2)

!!! warning "Deprecated"
    Use `clamp` from Julia's standard library instead.
"""
brcktd

@deprecate brcktd clamp

function _brckti(number, e1, e2)
    ccall((:brckti_c, libcspice), SpiceInt, (SpiceInt, SpiceInt, SpiceInt),
        number, e1, e2)
end

"""
    brckti(number, e1, e2)

!!! warning "Deprecated"
    Use `clamp` from Julia's standard library instead.
"""
brckti

@deprecate brckti clamp

function _bschoc(value, array, order)
    array_, m, n = chararray(array)
    order_ = SpiceInt.(order .- 1)
    idx = ccall((:bschoc_c, libcspice), SpiceInt,
        (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceInt}),
        value, m, n, array_, order_)
    handleerror()
    return idx == -1 ? idx : idx + 1
end

"""
    bschoc(value, array, order)

!!! warning "Deprecated"
    Use `findfirst(array .== value)` instead.
"""
bschoc

@deprecate bschoc(value, array, order) findfirst(array .== value)

function _bschoi(value, array, order)
    array_ = SpiceInt.(array)
    n = length(array)
    order_ = SpiceInt.(order .- 1)
    idx = ccall((:bschoi_c, libcspice), SpiceInt,
        (SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}),
        value, n, array_, order_)
    handleerror()
    return idx == -1 ? idx : idx + 1
end

"""
    bschoi(value, array, order)

!!! warning "Deprecated"
    Use `findfirst(array .== value)` instead.
"""
bschoi

@deprecate bschoi(value, array, order) findfirst(array .== value)

function _bsrchc(value, array)
    array_, m, n = chararray(array)
    idx = ccall((:bsrchc_c, libcspice), SpiceInt,
        (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
        value, m, n, array_)
    handleerror()
    return idx == -1 ? idx : idx + 1
end

"""
    bsrchc(value, array)

!!! warning "Deprecated"
    Use `findfirst(array .== value)` instead.
"""
bsrchc

@deprecate bsrchc(value, array) findfirst(array .== value)

function _bsrchd(value, array)
    n = length(array)
    idx = ccall((:bsrchd_c, libcspice), SpiceInt,
        (SpiceDouble, SpiceInt, Ref{SpiceDouble}),
        value, n, array)
    handleerror()
    return idx == -1 ? idx : idx + 1
end

"""
    bsrchd(value, array)

!!! warning "Deprecated"
    Use `findfirst(array .== value)` instead.
"""
bsrchd

@deprecate bsrchd(value, array) findfirst(array .== value)

function _bsrchi(value, array)
    array_ = SpiceInt.(array)
    n = length(array)
    idx = ccall((:bsrchi_c, libcspice), SpiceInt,
        (SpiceInt, SpiceInt, Ref{SpiceInt}),
        value, n, array_)
    handleerror()
    return idx == -1 ? idx : idx + 1
end

"""
    bsrchi(value, array)

!!! warning "Deprecated"
    Use `findfirst(array .== value)` instead.
"""
bsrchi

@deprecate bsrchi(value, array) findfirst(array .== value)


text:
###
### Sampler states
###

struct Emcee{space, E<:AMH.Ensemble} <: InferenceAlgorithm
    ensemble::E
end

function Emcee(n_walkers::Int, stretch_length=2.0)
    # Note that the proposal distribution here is just a Normal(0,1)
    # because we do not need AdvancedMH to know the proposal for
    # ensemble sampling.
    prop = AMH.StretchProposal(nothing, stretch_length)
    ensemble = AMH.Ensemble(n_walkers, prop)
    return Emcee{(), typeof(ensemble)}(ensemble)
end

struct EmceeState{V<:AbstractVarInfo,S}
    vi::V
    states::S
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::Model,
    spl::Sampler{<:Emcee};
    resume_from = nothing,
    kwargs...
)
    if resume_from !== nothing
        state = loadstate(resume_from)
        return AbstractMCMC.step(rng, model, spl, state; kwargs...)
    end

    # Sample from the prior
    n = spl.alg.ensemble.n_walkers
    vis = [VarInfo(rng, model, SampleFromPrior()) for _ in 1:n]

    # Update the parameters if provided.
    if haskey(kwargs, :init_params)
        for vi in vis
            initialize_parameters!(vi, kwargs[:init_params], spl)

            # Update log joint probability.
            model(rng, vi, SampleFromPrior())
        end
    end

    # Compute initial transition and states.
    transition = map(vis) do vi
        Transition(vi)
    end
    state = EmceeState(
        vis[1],
        map(vis) do vi
            # Transform to unconstrained space.
            DynamicPPL.link!(vi, spl)
            AMH.Transition(vi[spl], getlogp(vi))
        end
    )

    return transition, state
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Emcee},
    state::EmceeState;
    kwargs...
)
    # Generate a log joint function.
    vi = state.vi
    densitymodel = AMH.DensityModel(gen_logπ(vi, SampleFromPrior(), model))

    # Compute the next states.
    states = last(AbstractMCMC.step(rng, densitymodel, spl.alg.ensemble, state.states))

    # Compute the next transition and state.
    transition = map(states) do _state
        vi[spl] = _state.params
        DynamicPPL.invlink!(vi, spl)
        t = Transition(tonamedtuple(vi), _state.lp)
        DynamicPPL.link!(vi, spl)
        return t
    end
    newstate = EmceeState(vi, states)

    return transition, newstate
end

function AbstractMCMC.bundle_samples(
    samples::Vector{<:Vector},
    model::AbstractModel,
    spl::Sampler{<:Emcee},
    state::EmceeState,
    chain_type::Type{MCMCChains.Chains};
    save_state = false,
    kwargs...
)
    # Convert transitions to array format.
    # Also retrieve the variable names.
    params_vec = map(_params_to_array, samples)

    # Extract names and values separately.
    nms = params_vec[1][1]
    vals_vec = [p[2] for p in params_vec]

    # Get the values of the extra parameters in each transition.
    extra_vec = map(get_transition_extras, samples)

    # Get the extra parameter names & values.
    extra_params = (internals = extra_vec[1][1],)
    extra_values_vec = [e[2] for e in extra_vec]

    # Extract names & construct param array.
    nms = [nms; extra_params]
    parray = map(x -> hcat(x[1], x[2]), zip(vals_vec, extra_values_vec))
    parray = cat(parray..., dims=3)

    # Get the average or final log evidence, if it exists.
    le = getlogevidence(samples, state, spl)

    # Set up the info tuple.
    if save_state
        info = (model = model, sampler = spl, samplerstate = state)
    else
        info = NamedTuple()
    end

    # Concretize the array before giving it to MCMCChains.
    parray = MCMCChains.concretize(parray)

    # Chain construction.
    return MCMCChains.Chains(
        parray,
        nms,
        extra_params;
        evidence=le,
        info=info,
    ) |> sort
end


text:
# helper function
function readvalue(io, s::String, dtype=Int)
    l = strip(readline(io))
    occursin(s, l) || throw(InputError(filepath))
    parse(dtype, split(l)[end])
end

function readnettntp(filepath::String;
                     _dtypes = Dict(:init_node => Int, :term_node => Int, :capacity => Float64,
                                    :length => Float64, :free_flow_time => Float64, :power => Float64,
                                    :b => Float64, :speed => Float64, :toll => Float64, :link_type => Int))
    # reading metadata
    nzones, nnodes, ftnode, nlinks, headerrow = open(filepath) do io
        nzones = readvalue(io, "<NUMBER OF ZONES>")
        nnodes = readvalue(io, "<NUMBER OF NODES>")
        ftnode = readvalue(io, "<FIRST THRU NODE>")
        nlinks = readvalue(io, "<NUMBER OF LINKS>")

        headerrow = 4
        while(headerrow += 1; ((line = strip(readline(io))) == "") || !startswith(line, "~"))
            if eof(io)
                throw(InputError("Check network data."))
            end
        end

        return (nzones, nnodes, ftnode, nlinks, headerrow)
    end

    # reading network-link data
    linkdf = read(filepath, skipto=headerrow+1, header=headerrow, types=_dtypes)
    linkdf = linkdf[:,2:(end-1)]

    return (nzones, nnodes, ftnode, nlinks, linkdf)
end

function readtriptntp(filepath::String)
    nzones, totalflow, trips = open(filepath) do io
        # reading metadata
        nzones = readvalue(io, "<NUMBER OF ZONES>")
        totalflow = readvalue(io, "<TOTAL OD FLOW>", Float64)

        # reading trip data
        while (line = strip(readline(io)); ~occursin("Origin", line)) end

        trips = zeros(nzones, nzones)
        i = parse(Int, split(line)[end])
        while !eof(io)
            line = strip(readline(io))
            if occursin("Origin", line)
                i = parse(Int, split(line)[end])
            elseif (line == "")

            else
                line = line[1:(end - (line[end] == ';'))]
                for traw in split(line, ";")
                    ttmp = split(traw)
                    trips[i, parse(Int, ttmp[1])] = parse(Float64, ttmp[end])
                end
            end
        end

        return (nzones, totalflow, trips)
    end

    return (nzones, totalflow, trips)
end

function readnodetntp(filepath::String)
    data = read(filepath)
    data = data[:, 1:(end-1)]
    return data
end

readflowtntp(filepath::String) = read(filepath)

function readtntpdata(folderpath::String)
    contents = filter((f) -> endswith(f, ".tntp"), readdir(folderpath))
    filetypes = [:net, :trips, :flow, :node]

    filepaths = Dict()
    for ft in filetypes
        files = filter(x -> endswith(x, "_$ft.tntp"), contents)
        if length(files) > 1
            throw("More than one $ft files found in the folder!")
        elseif length(files) == 0
            throw("No $ft files found in the folder!")
        else
            f = files[1]
        end
        filepaths[ft] = folderpath*f
    end

    nzones, nnodes, ftnode, nlinks, linkdf = readnettntp(filepaths[:net])
    nzones2, totalflow, trips = readtriptntp(filepaths[:trips])

    geometry = nothing
    if :node in keys(filepaths)
        geometry = readnodetntp(filepaths[:node])
    end

    bestsolution = nothing
    if :flow in keys(filepaths)
        bestsolution = readflowtntp(filepaths[:flow])
        @assert nlinks == nrow(bestsolution) "Number of links inconsistent with flow file"
    end

    @assert nzones == nzones2 "Number of zones inconsistent"
    @assert totalflow ≈ sum(trips) "Number of trips inconsistent"
    @assert nlinks == nrow(linkdf) "Number of links inconsistent"
    @assert all((1 .<= linkdf.init_node .<= nnodes) .& (1 .<= linkdf.term_node .<= nnodes)) "Number of nodes inconsistent"
    @assert ftnode <= nnodes "First through node inconsistent"

    return (nnodes = nnodes,
            ftnode = ftnode,
            linkdf = linkdf,
            trips = trips,
            geometry = geometry,
            bestsolution = bestsolution)
end


text:
#####
##### Weighted Essentially Non-Oscillatory (WENO) schemes
#####

struct WENO{K, FT, C, B, Γ} <: AbstractAdvectionScheme
    C :: C
    B :: B
    γ :: Γ
    ε :: FT
    η :: FT
end

function WENO(FT, n)
    (n < 3 || iseven(n)) &&
        error("WENO schemes are only defined for order n = 3, 5, 7, 9, 11, ...")

    @warn "Generated WENO schemes are experimental! They may be very slow (allocate tons of memory) and may not work on GPUs."

    n >= 7 &&
        @info "Computing WENO-$n smoothness indicator coefficients. This could take a while for n >= 7..."

    k = Int((n + 1) / 2)
    C = eno_coefficients_matrix(FT, k)
    B = β_coefficients(FT, k)
    γ = optimal_weno_weights(FT, k)

    ε = 1e-6  # To avoid dividing by zero
    η = 2     # WENO exponent

    return WENO{k, FT, typeof(C), typeof(B), typeof(γ)}(C, B, γ, ε, η)
end

WENO(n) = WENO(Float64, n)

#####
##### ENO interpolants
#####

@inline px(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.C[ℓ+1, r+1] * ϕ[i-r-1+ℓ, j, k] for ℓ in 0:K-1)
@inline py(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.C[ℓ+1, r+1] * ϕ[i, j-r-1+ℓ, k] for ℓ in 0:K-1)
@inline pz(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.C[ℓ+1, r+1] * ϕ[i, j, k-r-1+ℓ] for ℓ in 0:K-1)

#####
##### Jiang & Shu (1996) WENO smoothness indicators
#####

@inline βx(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.B[m+1, n+1, r+1] * ϕ[i-r-1+m, j, k] * ϕ[i-r-1+n, j, k] for m in 0:K-1 for n in 0:m)
@inline βy(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.B[m+1, n+1, r+1] * ϕ[i, j-r-1+m, k] * ϕ[i, j-r-1+n, k] for m in 0:K-1 for n in 0:m)
@inline βz(i, j, k, weno::WENO{K}, r, ϕ) where K = @inbounds sum(weno.B[m+1, n+1, r+1] * ϕ[i, j, k-r-1+m] * ϕ[i, j, k-r-1+n] for m in 0:K-1 for n in 0:m)

#####
##### WENO raw weights
#####

@inline αx(i, j, k, weno, r, ϕ) = @inbounds weno.γ[r+1] / (βx(i, j, k, weno, r, ϕ) + weno.ε)^weno.η
@inline αy(i, j, k, weno, r, ϕ) = @inbounds weno.γ[r+1] / (βy(i, j, k, weno, r, ϕ) + weno.ε)^weno.η
@inline αz(i, j, k, weno, r, ϕ) = @inbounds weno.γ[r+1] / (βz(i, j, k, weno, r, ϕ) + weno.ε)^weno.η

#####
##### WENO normalized weights
#####

@inline function weno_weights_x(i, j, k, weno::WENO{K}, ϕ) where K
    α = [αx(i, j, k, weno, r, ϕ) for r in 0:K-1]
    return α ./ sum(α)
end

@inline function weno_weights_y(i, j, k, weno::WENO{K}, ϕ) where K
    α = [αy(i, j, k, weno, r, ϕ) for r in 0:K-1]
    return α ./ sum(α)
end

@inline function weno_weights_z(i, j, k, weno::WENO{K}, ϕ) where K
    α = [αz(i, j, k, weno, r, ϕ) for r in 0:K-1]
    return α ./ sum(α)
end

#####
##### WENO flux reconstruction
#####

@inline function weno_flux_x(i, j, k, weno::WENO{K}, ϕ) where K
    w = weno_weights_x(i, j, k, weno, ϕ)
    return sum(w[r+1] * px(i, j, k, weno, r, ϕ) for r in 0:K-1)
end

@inline function weno_flux_y(i, j, k, weno::WENO{K}, ϕ) where K
    w = weno_weights_y(i, j, k, weno, ϕ)
    return sum(w[r+1] * py(i, j, k, weno, r, ϕ) for r in 0:K-1)
end

@inline function weno_flux_z(i, j, k, weno::WENO{K}, ϕ) where K
    w = weno_weights_z(i, j, k, weno, ϕ)
    return sum(w[r+1] * pz(i, j, k, weno, r, ϕ) for r in 0:K-1)
end

#####
##### Momentum advection fluxes
#####

@inline momentum_flux_uu(i, j, k, grid, weno::WENO, u)    = ℑxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * weno_flux_x(i, j, k, weno, u)
@inline momentum_flux_uv(i, j, k, grid, weno::WENO, u, v) = ℑxᶠᵃᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * weno_flux_y(i, j, k, weno, u)
@inline momentum_flux_uw(i, j, k, grid, weno::WENO, u, w) = ℑxᶠᵃᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * weno_flux_z(i, j, k, weno, u)

@inline momentum_flux_vu(i, j, k, grid, weno::WENO, u, v) = ℑyᵃᶠᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * weno_flux_x(i, j, k, weno, v)
@inline momentum_flux_vv(i, j, k, grid, weno::WENO, v)    = ℑyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * weno_flux_y(i, j, k, weno, v)
@inline momentum_flux_vw(i, j, k, grid, weno::WENO, v, w) = ℑyᵃᶠᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * weno_flux_z(i, j, k, weno, v)

@inline momentum_flux_wu(i, j, k, grid, weno::WENO, u, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) * weno_flux_x(i, j, k, weno, w)
@inline momentum_flux_wv(i, j, k, grid, weno::WENO, v, w) = ℑzᵃᵃᶠ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) * weno_flux_y(i, j, k, weno, w)
@inline momentum_flux_ww(i, j, k, grid, weno::WENO, w)    = ℑzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, w) * weno_flux_z(i, j, k, weno, w)

#####
##### Advective tracer fluxes
#####

@inline advective_tracer_flux_x(i, j, k, grid, weno::WENO, u, c) = Ax_ψᵃᵃᶠ(i, j, k, grid, u) * weno_flux_x(i, j, k, weno, c)
@inline advective_tracer_flux_y(i, j, k, grid, weno::WENO, v, c) = Ay_ψᵃᵃᶠ(i, j, k, grid, v) * weno_flux_y(i, j, k, weno, c)
@inline advective_tracer_flux_z(i, j, k, grid, weno::WENO, w, c) = Az_ψᵃᵃᵃ(i, j, k, grid, w) * weno_flux_z(i, j, k, weno, c)

#####
##### Need to advect momentum like tracers
#####

@inline div_ũu(i, j, k, grid, advection::WENO, U) = div_uc(i, j, k, grid, advection, U, U.u)
@inline div_ũv(i, j, k, grid, advection::WENO, U) = div_uc(i, j, k, grid, advection, U, U.v)
@inline div_ũw(i, j, k, grid, advection::WENO, U) = div_uc(i, j, k, grid, advection, U, U.w)


text:
"""
```julia
BayesianSVM()
```

The [Bayesian SVM](https://arxiv.org/abs/1707.05532) is a Bayesian interpretation of the classical SVM.
```math
p(y|f) ∝ exp(2 max(1-yf,0))
````

---
For the analytic version of the likelihood, it is augmented via:
```math
p(y|f,ω) = 1/(sqrt(2πω) exp(-0.5((1+ω-yf)^2/ω))
```
where ``ω ∼ 𝟙[0,∞)`` has an improper prior (his posterior is however has a valid distribution, a Generalized Inverse Gaussian). For reference [see this paper](http://ecmlpkdd2017.ijs.si/papers/paperID502.pdf)
"""
struct BayesianSVM{T<:Real} <: ClassificationLikelihood{T}
    ω::AbstractVector{T}
    θ::AbstractVector{T}
    function BayesianSVM{T}() where {T<:Real}
        new{T}()
    end
    function BayesianSVM{T}(ω::AbstractVector{<:Real},θ::AbstractVector{<:Real}) where {T<:Real}
        new{T}(ω,θ)
    end
end

function BayesianSVM()
    BayesianSVM{Float64}()
end

function init_likelihood(likelihood::BayesianSVM{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer,nFeatures::Integer) where T
    BayesianSVM{T}(abs.(rand(T,nSamplesUsed)),zeros(T,nSamplesUsed))
end
function pdf(l::BayesianSVM,y::Real,f::Real)
    svmlikelihood(y*f)
end

function Base.show(io::IO,model::BayesianSVM{T}) where T
    print(io,"Bayesian SVM")
end

"""Return likelihood equivalent to SVM hinge loss"""
function svmlikelihood(f::Real)
    pos = svmpseudolikelihood(f)
    return pos./(pos.+svmpseudolikelihood(-f))
end

"""Return the pseudo likelihood of the SVM hinge loss"""
function svmpseudolikelihood(f::Real)
    return exp(-2.0*max.(1.0-f,0))
end

function compute_proba(l::BayesianSVM{T},μ::Vector{T},σ²::Vector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    sig_pred = zeros(T,N)
    for i in 1:N
        x = pred_nodes.*sqrt(max(σ²[i],zero(T))).+μ[i]
        pred[i] =  dot(pred_weights,svmlikelihood.(x))
        sig_pred[i] = max(dot(pred_weights,svmlikelihood.(x).^2)-pred[i]^2,zero(T))
    end
    return pred, sig_pred
end

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

function local_updates!(l::BayesianSVM{T},y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    l.ω .= abs2.(one(T) .- y.*μ) + diag_cov
    l.θ .= one(T)./sqrt.(l.ω)
end

@inline ∇E_μ(l::BayesianSVM{T},::AOptimizer,y::AbstractVector) where {T} = (y.*(l.θ.+one(T)),)
@inline ∇E_Σ(l::BayesianSVM{T},::AOptimizer,y::AbstractVector) where {T} = (0.5.*l.θ,)

function expec_log_likelihood(l::BayesianSVM{T},i::AnalyticVI,y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    tot = -(0.5*length(y)*logtwo)
    tot += dot(μ,y)
    tot += -0.5*dot(l.θ,diag_cov)+dot(l.θ,abs2.(one(T).-y.*μ))
    return tot
end

AugmentedKL(l::BayesianSVM,::AbstractVector) = GIGEntropy(l)

function GIGEntropy(l::BayesianSVM)
    return 0.5*sum(log.(l.ω))+sum(log.(2.0*besselk.(0.5,sqrt.(l.ω))))-0.5*sum(sqrt.(l.ω))
end


text:
"""
`InitialGuessPODE`: Initial guess for partitioned ordinary differential equations

### Fields

* `int`: interpolation structure
* `v`:   vector field for q
* `f`:   vector field for p
* `Δt`:  time step
* `s`:   number of extrapolation stages (for initialisation)
"""
mutable struct InitialGuessPODE{DT, TT, VT, FT, IT <: Interpolator}
    int::IT
    v::VT
    f::FT
    Δt::TT
    s::Int

    function InitialGuessPODE{DT,TT,VT,FT,IT}(interp, v, f, Δt) where {DT,TT,VT,FT,IT}
        new(interp, v, f, Δt, get_config(:ig_extrapolation_stages))
    end
end

function InitialGuessPODE(interp, equ::PODE{DT,TT,VT,FT}, Δt::TT) where {DT,TT,VT,FT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

function InitialGuessPODE(interp, equ::IODE{DT,TT,ϑT,FT,GT,VT}, Δt::TT) where {DT,TT,ϑT,FT,GT,VT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

function InitialGuessPODE(interp, equ::VODE{DT,TT,ϑT,FT,GT,VT}, Δt::TT) where {DT,TT,ϑT,FT,GT,VT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

function InitialGuessPODE(interp, equ::IDAE{DT,TT,ϑT,FT,UT,GT,ϕT,VT}, Δt::TT) where {DT,TT,ϑT,FT,UT,GT,ϕT,VT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

function InitialGuessPODE(interp, equ::PDAE{DT,TT,VT,FT,UT,GT,ϕT}, Δt::TT) where {DT,TT,VT,FT,UT,GT,ϕT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

function InitialGuessPODE(interp, equ::VDAE{DT,TT,θT,FT,GT,G̅T,ϕT,ψT,VT}, Δt::TT) where {DT,TT,θT,FT,GT,G̅T,ϕT,ψT,VT}
    InitialGuessPODE{DT,TT,VT,FT,interp}(interp(zero(DT), one(DT), Δt, equ.d), equ.v, equ.f, Δt)
end

"Initialise initial guess, i.e., given t₀, t₁, q₁ compute q₀, v₀, v₁."
function initialize!(ig::InitialGuessPODE{DT,TT,VT,IT},
                t₁::TT,
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                f₁::SolutionVector{DT},
                t₀::TT,
                q₀::SolutionVector{DT},
                p₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                f₀::SolutionVector{DT}) where {DT,TT,VT,IT}

    midpoint_extrapolation(ig.v, ig.f, t₁, t₀, q₁, q₀, p₁, p₀, ig.s)

    ig.v(t₀, q₀, p₀, v₀)
    ig.v(t₁, q₁, p₁, v₁)
    ig.f(t₀, q₀, v₀, f₀)
    ig.f(t₁, q₁, v₁, f₁)
end

function update!(ig::InitialGuessPODE{DT,TT}, t₁::TT,
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::Vector{DT},
                f₁::Vector{DT}) where {DT,TT}
    ig.v(t₁, q₁, p₁, v₁)
    ig.f(t₁, q₁, v₁, f₁)
end

function CommonFunctions.evaluate!(ig::InitialGuessPODE{DT,TT},
                q₀::SolutionVector{DT},
                p₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                f₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                f₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                c_q::TT) where {DT,TT}

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁."
        guess_q .= q₁
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c_q, guess_q)
    end
end

function CommonFunctions.evaluate!(ig::InitialGuessPODE{DT,TT},
                q₀::SolutionVector{DT},
                p₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                f₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                f₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                guess_v::SolutionVector{DT},
                c_q::TT) where {DT,TT}
    @assert length(guess_q) == length(guess_v)

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁ and v=0."
        guess_q .= q₁
        guess_v .= 0
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c_q, guess_q, guess_v)
    end
end

function CommonFunctions.evaluate!(ig::InitialGuessPODE{DT,TT},
                q₀::SolutionVector{DT},
                p₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                f₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                f₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                guess_p::SolutionVector{DT},
                c_q::TT, c_p::TT) where {DT,TT}
    @assert length(guess_q) == length(guess_p)

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁."
        guess_q .= q₁
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c_q, guess_q)
    end

    if p₀ == p₁
        @warn "p₀ and p₁ in initial guess are identical! Setting p=p₁."
        guess_p .= p₁
    else
        evaluate!(ig.int, p₀, p₁, f₀, f₁, one(TT)+c_p, guess_p)
    end
end

function CommonFunctions.evaluate!(ig::InitialGuessPODE{DT,TT},
                q₀::SolutionVector{DT},
                p₀::SolutionVector{DT},
                v₀::SolutionVector{DT},
                f₀::SolutionVector{DT},
                q₁::SolutionVector{DT},
                p₁::SolutionVector{DT},
                v₁::SolutionVector{DT},
                f₁::SolutionVector{DT},
                guess_q::SolutionVector{DT},
                guess_p::SolutionVector{DT},
                guess_v::SolutionVector{DT},
                guess_f::SolutionVector{DT},
                c_q::TT, c_p::TT) where {DT,TT}
    @assert length(guess_q) == length(guess_v)
    @assert length(guess_p) == length(guess_f)
    @assert length(guess_q) == length(guess_p)

    if q₀ == q₁
        @warn "q₀ and q₁ in initial guess are identical! Setting q=q₁ and v=0."
        guess_q .= q₁
        guess_v .= 0
    else
        evaluate!(ig.int, q₀, q₁, v₀, v₁, one(TT)+c_q, guess_q, guess_v)
    end

    if p₀ == p₁
        @warn "p₀ and p₁ in initial guess are identical! Setting p=p₁ and f=0."
        guess_p .= p₁
        guess_f .= 0
    else
        evaluate!(ig.int, p₀, p₁, f₀, f₁, one(TT)+c_p, guess_p, guess_f)
    end
end


text:
using DifferentialEquations
const linear_bigα = big(1.01)
f(u,p,t) = (linear_bigα*u)

# Add analytical solution so that errors are checked
f_analytic(u0,p,t) = u0*exp(linear_bigα*t)
ff = ODEFunction(f,analytic=f_analytic)
prob = ODEProblem(ff,big(0.5),(0.0,1.0))
sol = solve(prob,Feagin14(),dt=1//16,adaptive=false);

println(sol.errors)

eps(Float64)

sol =solve(prob,Feagin14());
println(sol.errors); print("The length was $(length(sol))")

using DiffEqDevTools
dts = 1.0 ./ 2.0 .^(10:-1:4)
sim = test_convergence(dts,prob,Feagin14())

using Plots
gr()
plot(sim)

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using PlotlyJS

export drawMolGMM, drawMolGMMs, drawmol, drawPharmacophoreGMMs, plotdrawing

const atom_colors =
    Dict(:C  => "#383838",   # dark grey
         :H  => "#b5b5b5",   # light grey
         :O  => "#d62728",   # red
         :N  => "#1f77b4",   # blue
         :S  => "#cbd123",   # yellow
         :Cl => "#2ca02c",   # green
         :P  => "#e29522")   # orange

const feature_colors =
    Dict(:donor         => "#ff00ff",  # magenta
         :acceptor      => "#00ff00",  # green
         :cation        => "#ff0000",  # red
         :anion         => "#0000ff",  # blue
         :hydrophobe    => "#00ffff",  # cyan
         :rings         => "#ffa500",  # orange
         :aromaticrings => "#964b00")  # brown

function drawatoms(atoms::AbstractVector{SDFileAtom}; colordict=atom_colors, markersize=10)
    if isempty(atoms)
        return
    end
    color = colordict[atoms[1].symbol]  # all atoms assumed to be of the same element
    xs = [a.coords[1] for a in atoms]
    ys = [a.coords[2] for a in atoms]
    zs = [a.coords[3] for a in atoms]
    return scatter3d(;x=xs, y=ys, z=zs,
                     mode="markers", marker=attr(color=color, size=markersize, symbol="circle"),
                     showlegend=false, showscale=false, name="", hoverinfo="skip",
                     hovertemplate="")
end

function drawbonds(nodes::AbstractVector{SDFileAtom}, edges, edgeidxs, bonds::AbstractVector{SDFileBond}, planedir=SVector{3}([0.,0.,1.]); name = "bonds", linewidth=3, separation=0.2)
    xs, ys, zs = Float64[], Float64[], Float64[]
    for eidx in edgeidxs
        edge = edges[eidx]
        bondorder = bonds[eidx].order
        pos1, pos2 = nodes[edge[1]].coords, nodes[edge[2]].coords
        sepvec = cross(planedir, pos1.-pos2)
        sepvec = sepvec/norm(sepvec)
        if isodd(bondorder)
            d=0.
        else
            d=0.5
        end
        for i=1:Int(ceil(bondorder/2))
            dvec = (d+1-i)*separation*sepvec
            p1, p2 = pos1 .+ dvec, pos2 .+ dvec
            push!(xs, p1[1], p2[1], NaN)
            push!(ys, p1[2], p2[2], NaN)
            push!(zs, p1[3], p2[3], NaN)
            if (d+1-i)>0
                p1, p2 = pos1 .- dvec, pos2 .- dvec
                push!(xs, p1[1], p2[1], NaN)
                push!(ys, p1[2], p2[2], NaN)
                push!(zs, p1[3], p2[3], NaN)
            end
        end
    end

    return scatter3d(;x=xs, y=ys, z=zs,
                      mode="lines", line=attr(color="black", width=linewidth),
                      showlegend=!isnothing(name), name=isnothing(name) ? "" : name, hoverinfo="skip",
                      hovertemplate="")
end

drawbonds(mol::UndirectedGraph, planedir=SVector{3}([0.,0.,1.]); kwargs...
    ) = drawbonds(nodeattrs(mol), mol.edges, 1:length(mol.edges), edgeattrs(mol), planedir; kwargs...)

drawbonds(sg::SubgraphView, planedir=SVector{3}([0.,0.,1.]); kwargs...
    ) = drawbonds(nodeattrs(sg), sg.graph.edges, sg.edges, edgeattrs(sg), planedir; kwargs...)

function drawmol(mol::Union{GraphMol,UndirectedGraph,SubgraphView}; colordict=atom_colors, atoms=true, bonds=true, hydrogens=false, markersize=10, linewidth=1)
    traces = AbstractTrace[]
    if !hydrogens == true
        if typeof(mol) <: GraphMol
            mol = removehydrogens(mol)
        elseif typeof(mol) <: SubgraphView
            mol = nodesubgraph(removehydrogens(mol.graph), nodeset(removehydrogens(mol.graph)) ∩ nodeset(mol))
        end
    end

    formula = isa(mol,SubgraphView) ? molecularformula(mol.graph) : molecularformula(mol)
    elements = Symbol[]
    for idx in eachindex(formula)
        if isletter(formula[idx])
            # the last character should always be a number
            if idx < length(formula)
                if isletter(formula[idx+1])
                    push!(elements, Symbol(formula[idx:idx+1]))
                else
                    push!(elements, Symbol(formula[idx]))
                end
            else
                push!(elements, Symbol(formula[idx]))
            end
        end
    end

    coordmat = fill(NaN, 3, length(nodeset(mol)))
    # each atom type (CHONS, in most cases) gets its own trace to differentiate between colors, while keeping the number of traces low
    for el in elements
        elatoms = SDFileAtom[]
        for idx in nodeset(mol)
            if nodeattr(mol,idx).symbol == el
                coordmat[:,idx] = nodeattr(mol,idx).coords
                push!(elatoms, nodeattr(mol,idx))
            end
        end
        if atoms && !isempty(elatoms)
            push!(traces, drawatoms(elatoms; colordict=colordict, markersize=markersize))
        end
    end

    if bonds
        push!(traces, drawbonds(mol; linewidth=linewidth))
    end
    return traces
end

function drawMolGMM(molgmm::MolGMM; color=GaussianMixtureAlignment.default_colors[1], atoms=true, bonds=true, hydrogens=false, markersize=10, linewidth=1, kwargs...)
    sg = nodesubgraph(molgmm.graph, molgmm.nodes)
    traces = drawmol(sg; atoms=atoms, bonds=bonds, hydrogens=hydrogens, markersize=markersize, linewidth=linewidth)
    push!(traces, drawIsotropicGMM(molgmm; color=color, kwargs...)...)
    return traces
end

function drawMolGMMs(molgmms::AbstractVector{<:MolGMM}; colors=GaussianMixtureAlignment.default_colors, kwargs...)
    traces = AbstractTrace[]
    for (i,molgmm) in enumerate(molgmms)
        color = colors[mod(i-1, length(colors))+1]
        push!(traces, drawMolGMM(molgmm; color=color, kwargs...)...)
    end
    return traces
end

function drawPharmacophoreGMM(fgmm::PharmacophoreGMM; colordict=feature_colors, atoms=true, bonds=true, hydrogens=false, markersize=10, linewidth=1, kwargs...)
    sg = nodesubgraph(fgmm.graph, fgmm.nodes)
    traces = drawmol(sg; atoms=atoms, bonds=bonds, hydrogens=hydrogens, markersize=markersize, linewidth=linewidth)
    push!(traces, drawMultiGMM(fgmm; colordict=colordict, kwargs...)...)
    return traces
end

function drawPharmacophoreGMMs(fgmms::AbstractVector{<:PharmacophoreGMM}; colordict=feature_colors, kwargs...)
    traces = AbstractTrace[]
    for fgmm in fgmms
        push!(traces, drawPharmacophoreGMM(fgmm; colordict=colordict, kwargs...)...)
    end
    return traces
end


text:
############################################################################################
#                              DYNAMICS CONSTRAINTS										   #
############################################################################################

abstract type AbstractDynamicsConstraint <: CoupledConstraint end
sense(::AbstractDynamicsConstraint) = Equality()
state_dim(con::AbstractDynamicsConstraint) = size(con.model)[1]
control_dim(con::AbstractDynamicsConstraint) = size(con.model)[2]
Base.length(con::AbstractDynamicsConstraint) = size(con.model)[1]

""" $(TYPEDEF)
An equality constraint imposed by the discretized system dynamics. Links adjacent time steps.
Supports both implicit and explicit integration methods. Can store values internally for
more efficient computation of dynamics and dynamics Jacobians over the entire trajectory,
particularly for explicit methods. These constraints are used in Direct solvers, where
the dynamics are explicit stated as constraints in a more general optimization method.

# Constructors
```julia
DynamicsConstraint{Q}(model::AbstractModel, N)
```
where `N` is the number of knot points and `Q<:QuadratureRule` is the integration method.
"""
struct DynamicsConstraint{Q<:QuadratureRule,L<:AbstractModel,N,M,NM,T} <: AbstractDynamicsConstraint
	model::L
    fVal::Vector{SVector{N,T}}
    xMid::Vector{SVector{N,T}}
    ∇f::Vector{SizedMatrix{N,NM,T,2,Matrix{T}}}
	A::Vector{SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	B::Vector{SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	grad::Vector{GradientExpansion{T,N,M}}
	∇fMid::Vector{SizedMatrix{N,NM,T,2,Matrix{T}}}
	Am::Vector{SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
	Bm::Vector{SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
end

function DynamicsConstraint{Q}(model::L, N) where {Q,L}
	T = Float64  # TODO: get this from somewhere
	n,m = size(model)
	fVal = [@SVector zeros(n) for k = 1:N]
	xMid = [@SVector zeros(n) for k = 1:N]
	∇f   = [SizedMatrix{n,n+m}(zeros(n,n+m)) for k = 1:N]
	∇fm  = [SizedMatrix{n,n+m}(zeros(n,n+m)) for k = 1:3]
	ix,iu = 1:n, n .+ (1:m)
	A  = [view(∇f[k].data, ix,ix) for k = 1:N]
	B  = [view(∇f[k].data, ix,iu) for k = 1:N]
	Am = [view(∇fm[k].data,ix,ix) for k = 1:3]
	Bm = [view(∇fm[k].data,ix,iu) for k = 1:3]
	grad  = [GradientExpansion{T}(n,m) for k = 1:3]
	DynamicsConstraint{Q,L,n,m,n+m,T}(model, fVal, xMid, ∇f, A, B,
		grad, ∇fm, Am, Bm)
end

@inline DynamicsConstraint(model, N) = DynamicsConstraint{DEFAULT_Q}(model, N)
integration(::DynamicsConstraint{Q}) where Q = Q

widths(con::DynamicsConstraint{<:Any,<:Any,N,M},n::Int=N,m::Int=M) where {N,M} = (n+m,n+m)
widths(con::DynamicsConstraint{<:Explicit,<:Any,N,M},n::Int=N,m::Int=M) where {N,M} = (n+m,n)

get_inds(con::DynamicsConstraint{<:Explicit}, n, m) = (1:n+m, (n+m) .+ (1:n))

# Explict
function evaluate(con::DynamicsConstraint{Q}, z1::AbstractKnotPoint, z2::AbstractKnotPoint) where Q <: Explicit
	RobotDynamics.discrete_dynamics(Q, con.model, z1) - state(z2)
end

function jacobian!(∇c, con::DynamicsConstraint{Q,L},
		z::AbstractKnotPoint, z2::AbstractKnotPoint{<:Any,n}, i=1) where {Q,L,n}
	if i == 1
		RobotDynamics.discrete_jacobian!(Q, ∇c, con.model, z)
		return L <: RD.LinearModel
	elseif i == 2
		for i = 1:n
			∇c[i,i] = -1
		end
		return true   # is constant
	end
	# return nothing
end

function ∇jacobian!(G, con::DynamicsConstraint{<:Explicit},
	z::AbstractKnotPoint, z2::AbstractKnotPoint, λ, i=1)
	if i == 1
		dyn(x) = evaluate(con, StaticKnotPoint(z,x), z2)'λ
		G .+= ForwardDiff.hessian(dyn, z.z)
	elseif i == 2
		nothing
	end
	return false
end


text:
using Combinatorics
using DataStructures

# implementation of J. De Leeuw Multivariate Cumulants in R (2012)

function outer(a::Vector{T}, b::Vector{T}) where T<: AbstractFloat
    sa = size(a,1)
    sb = size(b,1)
    z = zeros(T, sa*sb)
    k = 1
    for i = 1:sb
        for j = 1:sa
            @inbounds z[k] = a[j]*b[i]
            k += 1
        end
    end
    return z
end

function setparts(n)
    ret = zeros(Int, n, length(partitions(collect(1:n))))
    for (i, partition) in enumerate(partitions(collect(1:n)))
        group = 1
        for p in partition
            for j in p
                ret[j,i] = group
            end
            group += 1
        end
    end
    ret
end

function raw_moments_upto_p(x, p=4)
    n, m = size(x)
    if p==1
        return vcat(1, mean(x, axis=2))
    end
    y = zeros(eltype(x), (m+1)^p)
    for i in 1:n
        xi = vcat(1, x[i,:])
        #z = kron([xi for _ in 1:p]...)
        z = xi
        for j in 2:p
          z = outer(xi, z)
        end
        y += z
    end
    reshape(y, repeat([m+1],p)...)./n
end

function cumulants_from_raw_moments(raw)
        dimr = size(raw)
        nvar::Int64 = dimr[1]
        cumu = zeros(eltype(raw), dimr...)
        nele = prod(dimr)
        ldim = length(dimr)
        spp = Array{Array{Int64,2}}(undef, ldim)
        qpp::Array{Int64} = Array{Int64}(undef, ldim)
        rpp = Array{Any}(undef, ldim)

        function one_cumulant_from_raw_moments(jnd, raw)
            jnd = [jnd[find(jnd.!=1)]...] - 1
            nnd = length(jnd)
            ndr::Int64 = size(raw)[1]
            nrt = length(size(raw))
            raw = rpp[nnd]
            nvar = ndr - 1
            nraw = max(1, length(size(raw)))
            sp = spp[nraw]
            _, nbell = size(sp)
            sterm = 0.0
            for i in 1:nbell
                ind = sp[:, i]
                und = unique(ind)
                term = qpp[length(und)]
                for j in und
                    knd = jnd[find(ind.==j)] + 1
                    lnd =  vcat(knd, repeat([1], nraw - length(knd)))
                    term *= raw[lnd...]
                end
            sterm +=  term
            end
            return sterm
        end

        for i in 1:ldim
            spp[i] = setparts(i)
            qpp[i] = factorial(i)
            if mod(i,2)==1
                qpp[i] = -qpp[i]
            end
            rpp[i] = raw[hcat([collect(1:nvar) for i in 1:i])...]
        end
        qpp = vcat(1, qpp)
        for i in 2:nele
            ind = ind2sub(dimr, i)
            cumu[i] = one_cumulant_from_raw_moments(ind, raw)
        end
        return cumu
end

function cumulants_upto_p(x, p = 4)
    return cumulants_from_raw_moments(raw_moments_upto_p(x, p))
end

function first_four_cumulants(x)
    cumu = cumulants_upto_p(x)
    nsel = 2:size(cumu)[1]
    OrderedDict(:c1 => cumu[1, 1, 1, nsel],
                :c2 => cumu[1, 1, nsel, nsel],
                :c3 => cumu[1, nsel, nsel, nsel],
                :c4 => cumu[nsel, nsel, nsel, nsel]
                )
end


text:
        # boundingbox
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 0)
        𝘅₃ = Point2D{T}(1, 1)
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        bb = boundingbox(q)
        @test bb.xmin ≈ 0
        @test bb.ymin ≈ 0
        @test bb.xmax ≈ 2
        @test bb.ymax ≈ 1
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 2)
        𝘅₃ = Point2D{T}(1, 1)
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        bb = boundingbox(q)
        @test bb.xmin ≈ 0
        @test bb.ymin ≈ 0
        @test bb.xmax ≈ 2
        @test bb.ymax ≈ 2
        𝘅₁ = Point2D{T}(0, 0)
        𝘅₂ = Point2D{T}(2, 0)
        𝘅₃ = Point2D{T}(2.1, 1)
        q = QuadraticSegment2D(𝘅₁, 𝘅₂, 𝘅₃)
        bb = boundingbox(q)
        @test bb.xmin ≈ 0
        @test bb.ymin ≈ 0
        @test bb.xmax ≈ 2.3272727272727276
        @test bb.ymax ≈ 1

            # union
            aab = union(AABox2D(Point2D{T}(0,0), Point2D{T}(2, 2)),
                        AABox2D(Point2D{T}(1,1), Point2D{T}(3, 3)))
            @test aab ≈ AABox2D(Point2D{T}(0,0), Point2D{T}(3, 3))

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
## Exports

export AbstractSLAM
export ManageSolveSettings, HandlerStateMachine, SolverStateMachine
export manageSolveTree!, stopManageSolveTree!, blockSolvingInProgress
export SLAMCommonHelper, SLAMWrapperLocal
export checkSolveStrideTrigger!, checkSolveStride, triggerSolve!, blockProgress

## types

abstract type AbstractSLAM end

@enum HandlerStateMachine HSMReady HSMHandling HSMOverlapHandling HSMBlocking
@enum SolverStateMachine SSMReady SSMConsumingSolvables SSMSolving

"""
$(TYPEDEF)

DevNotes
- TODO consolidate with SLAMWrapperLocal
"""
mutable struct SLAMWrapper <: AbstractSLAM
  fg::AbstractDFG
  tree::NothingUnion{IncrementalInference.AbstractBayesTree}
  lndmidx::Int
  lastposesym::Symbol
  lastpose::SE3
  lbl2usrid::Dict{Symbol,Int}
  usrid2lbl::Dict{Int, Symbol}
  SLAMWrapper() = new()
  SLAMWrapper(a,b,c) = new(a,b,c,:x1, SE3(0), Dict{Symbol,Int}(), Dict{Int,Symbol}()) # TODO to be deprecated
end

"""
    $TYPEDEF

Container for settings and synchronization tools for use with manageSolveTree.
"""
mutable struct ManageSolveSettings
  solveStride
  loopSolver
  solvables
  solveInProgress
  poseSolveToken
  canTakePoses::Condition
  drtCurrent
  # constructor
  ManageSolveSettings(;solveStride=10,
                       loopSolver=true,
                       solvables=Channel{Vector{Symbol}}(100),
                       solveInProgress=SSMReady,
                       poseSolveToken=Channel{Int}(2), #  ensure only one solve can occur at a time (see blockProgress)
                       canTakePoses::Condition=Condition(),
                       drtCurrent=(:null,:null)) = new(solveStride, loopSolver, solvables, solveInProgress, poseSolveToken, canTakePoses, drtCurrent)
end

mutable struct SLAMCommonHelper
  lastPoseOdomBuffer::SE3 # common helper
  SLAMCommonHelper(lpo::SE3=SE3(0)) = new(lpo)
end

mutable struct SLAMWrapperLocal{G <: AbstractDFG} <: AbstractSLAM
  dfg::G
  poseCount::Int
  frameCounter::Int
  poseStride::Int # pose every frameStride (naive pose trigger)
  helpers::SLAMCommonHelper
  solveSettings::ManageSolveSettings
end

SLAMWrapperLocal(;dfg::G=initfg(),
                  poseCount::Int=0,
                  frameCounter::Int=0,
                  poseStride::Int=10,
                  helpers::SLAMCommonHelper=SLAMCommonHelper(),
                  solveSettings::ManageSolveSettings=ManageSolveSettings() ) where {G <: AbstractDFG}= SLAMWrapperLocal{G}(dfg,
                      poseCount, frameCounter, poseStride, helpers, solveSettings)
#

"""
    $SIGNATURES

Trigger a factor graph `solveTree!(slam.dfg,...)` after clearing the solvable buffer `slam.??` (assuming the `manageSolveTree!` task is already running).

Notes
- Used in combination with `manageSolveTree!`
"""
triggerSolve!(slam::SLAMWrapperLocal) = put!(slam.solveSettings.poseSolveToken, slam.poseCount)

function checkSolveStride(slam::SLAMWrapperLocal)
  slam.poseCount % slam.solveSettings.solveStride == 0
end

"""
    $SIGNATURES

Check and trigger a solve if `slam.poseCount` reached the solve stride `slam.solveSettings.solveStride`.

Notes
- Used in combination with `manageSolveTree!`

Related

triggerSolve!
"""
function checkSolveStrideTrigger!(slam::SLAMWrapperLocal; force::Bool=false)
  if force || checkSolveStride(slam)
    @info "trigger a solve, $(length(slam.solveSettings.poseSolveToken.data)) ====================================="
    triggerSolve!(slam)
    @info "after slam solve trigger"
    # also update DRT containers

    return true
  end
  return false
end

"""
    $SIGNATURES

Block progress on current task until SLAMWrapper is ready to continue taking on more data.

Notes:
- Use to block addition of more data into slam.dfg object, usually used for preventing duplicate solutions being invoked concurrently.
  - Consider dfg segments: 1. already busy solving, 2. recently added; and now attempting to launch solve on 1&2 before solve on 1. finishes.
- Also see dead reckon tether (DRT).
- Assumes unsolved poses are always added and never removed before being solved.

Related

manageSolveTree!, SLAMWrapperLocal
"""
function blockProgress(slam::SLAMWrapperLocal)
  mss = slam.solveSettings
  # determine if solve is in progress
  # mss.solveInProgress

  # isready means solve is in progress
  if 1 < length(mss.poseSolveToken.data) && checkSolveStride(slam)
    @warn "blockProgress called and forced to wait on previous slam solve not having completed."
    wait(mss.canTakePoses)
  end
end

"""
    $SIGNATURES

Stops a `manageSolveTree!` session.  Usually up to the user to do so as a SLAM process comes to completion.

Related

manageSolveTree!
"""
function stopManageSolveTree!(slam::SLAMWrapperLocal)
  slam.solveSettings.loopSolver = false
end

"""
    $SIGNATURES

Block the progression of calling task if `::SLAMWrapperLocal` is being solved by a presumed `manageSolveTree!` task.
"""
function blockSolvingInProgress(slam::SLAMWrapperLocal)
  if slam.solveSettings.solveInProgress != SSMReady
    @info "blockSolvingInProgress on ::SLAMWrapperLocal"
    wait(slam.solveSettings.canTakePoses)
    @info "blockSolvingInProgress notified of completion."
  end
end

"""
    $SIGNATURES

Asynchronous solver manager that can run concurrently while other Tasks are modifying a common distributed factor graph object.

Notes
- When adding Variables and Factors, use `solvable=0` to disable the new fragments until ready for inference.
  - e.g. `addVariable!(fg, :x45, Pose2, solvable=0)`
  - These parts of the factor graph can simply be activated for solving `setSolvable!(fg, :x45, 1)`
"""
function manageSolveTree!(dfg::AbstractDFG,
                          mss::ManageSolveSettings;
                          dbg::Bool=false,
                          timinglog=Base.stdout,
                          limitfixeddown::Bool=true  )
  #
  @info "logpath=$(getLogPath(dfg))"

  # allow async process
  # getSolverParams(dfg).async = true

  # prep with empty tree
  tree = emptyBayesTree()

  # needs to run asynchronously
  ST = @async begin
    while @show length(ls(dfg, :x0, solvable=1)) == 0
      "waiting for prior on x0" |> println
      sleep(1)
    end
    solvecycle = 0
    # keep solving
    while mss.loopSolver
      t0 = time_ns()
      solvecycle += 1
      # add any newly solvables (atomic)
      while !isready(mss.solvables) && mss.loopSolver && !isready(mss.poseSolveToken)
        @info "waiting on !isready(mss.solvables)=$(!isready(mss.solvables)) and mss.loopSolver=$(mss.loopSolver)"
        sleep(0.2)
      end
      dt_wait = (time_ns()-t0)/1e9

      #add any new solvables
      while isready(mss.solvables) && mss.loopSolver
        mss.solveInProgress = SSMConsumingSolvables
        @show tosolv = take!(mss.solvables)
        for sy in tosolv
          # setSolvable!(dfg, sy, 1) # see DFG #221
          # TODO temporary workaround
          getfnc = occursin(r"f", string(sy)) ? getFactor : getVariable
          getfnc(dfg, sy).solvable = 1
        end
      end
      dt_solvable = (time_ns()-t0)/1e9

      @info "Ensure all new variables initialized"
      ensureAllInitialized!(dfg)
      dt_init = (time_ns()-t0)/1e9
      dt_disengage = 0.0
      dt_save1 = 0.0
      dt_solve = 0.0

      # mss.solveInProgress = SSMReady

      # solve only every 10th pose
      @show length(mss.poseSolveToken.data)
      if 0 < length(mss.poseSolveToken.data)
      # if 10 <= mss[:poseStride]
        @info "reduce problem size by disengaging older parts of factor graph"
        setSolvableOldPoses!(dfg, youngest=getSolverParams(dfg).qfl+round(Int,mss.solveStride/2), oldest=100, solvable=0)
        dt_disengage = (time_ns()-t0)/1e9

        # set up state machine flags to allow overlapping or block
        mss.solveInProgress = SSMSolving # obsolete
        # mss[:poseStride] = 0

        # do the actual solve (with debug saving)
        lasp = getLastPoses(dfg, filterLabel=r"x\d", number=1)[1]
        !dbg ? nothing : saveDFG(dfg, joinpath(getLogPath(dfg), "fg_before_$(lasp)"))
        dt_save1 = (time_ns()-t0)/1e9
        # constrain solve with the latest pose at the top
        # @show latestPose = intersect(getLastPoses(dfg, filterLabel=r"x\d", number=12), ls(dfg, r"x\d", solvable=1))[end]
        tree, smt, hist = solveTree!(dfg, tree) # , variableConstraints=[latestPose;]
        dt_solve = (time_ns()-t0)/1e9
        !dbg ? nothing : saveDFG(dfg, joinpath(getLogPath(dfg), "fg_after_$(lasp)"))

        # unblock LCMLog reader for next STRIDE segment
        mss.solveInProgress = SSMReady

        # adjust latest RTT after solve, latest solved -- hard coded pose stride 10
        lastList = sortDFG(ls(dfg, r"x\d+9\b|x9\b", solvable=1))
        if 0 < length(lastList)
          lastSolved = lastList[end]
          mss.drtCurrent = (lastSolved, Symbol("drt_"*string(lastSolved)[2:end]))
        end

        # remove a token to allow progress to continue
        gotToken = take!(mss.poseSolveToken)

        # notify any processes that might be waiting on current solve to complete
        # lock(mss.canTakePoses)  # uncomment when upgrading to Threads.Condition
        notify(mss.canTakePoses) # = HSMHandling
        # unlock(mss.canTakePoses)

        "end of solve cycle, token=$gotToken" |> println
      else
        "sleep a solve cycle" |> println
        mss.solveInProgress = SSMReady
        sleep(0.2)
      end

      dt_finish = (time_ns() - t0)/1e9

      # store timing results
      println(timinglog, "$solvecycle, $t0, $dt_wait, $dt_solvable, $dt_init, $dt_disengage, $dt_save1, $dt_solve, $(tree.buildTime), $dt_finish")
    end
  end
  return ST
end

#


text:
n = 10
EmError = ones(5)
νm = [456812, 456808,456811, 456802]
Em = ones(length(νm))
σ = 20.0
γRate = 1.0e8
source = LightSource(n,Em,νm,σ,γRate)

# test LightSource assertion errors
@test_throws AssertionError LightSource(n,EmError,νm,σ,γRate)
@test_throws AssertionError LightSource(n,Em,νm,σ,456900,456800,γRate)

# test derived quantity functions
@test ν0(source) == 456808.25
@test Δm(source) == [3.75, -0.25, 2.75, -6.25]
@test lineShifts(source) == [4, 1, 10, 3, 6, 9]

# test eField construction
field = eField(source,42)
ν1 = field.νn[1]
ϕ11 = field.ϕmn[1,1]
ν1test = 456797.12946247705
ϕ11test = 5.89027451092873
@test ν1≈ν1test
@test ϕ11≈ϕ11test

@test intensity(0.5) == 0.25
@test intensity(0.01,field) ≈ 6.060748222475236

t = 10.0 #nanoseconds
dt = 0.1 #nanoseconds
nbar = 100.0
intensityVec = map(time->intensity(time,field),0:dt:t)
γVec = γIntensity(nbar,intensityVec)
@test sum(γVec.γvec)≈nbar

bs =  Beamsplitter(1,1)
@test bs.r^2+bs.t^2 ≈ 1.0


text:
function optimisedTauLeap!(model::Model, p::Path, tf::Float64, tauOpt::TauOpt, tauVar::TauVar)

    # ---- Main loop
    fly = zeros(model.nr)
    while p.t < tf
        moleculeLimit!(p, tf) && break
        model.F(p)
        if sum(p.a) <= 0.0
            break
        end
        tr = tf - p.t

        τ, switch = computeLeap(model, p, tr, tauOpt, tauVar)

        # ---- Switch to direct method?
        if switch == true
            nDirectMethod!(p, model.F, model.ν, tf, model.nr, tauOpt.nd)

        # ---- Else accept step
        else
            for i in eachindex(p.x)
                @inbounds p.x[i] += tauVar.dx[i]
            end

            @simd for i in eachindex(fly)
                @inbounds fly[i] += p.a[i]*τ
            end
            p.t += τ
        end
    end
    p.fly += fly./p.θ[1:4]
end

optimisedTauLeap = optimisedTauLeap!

function gUpdate!(hor::Vector{Int}, g::Vector{Float64}, x::Vector{Int})
    for i in eachindex(hor)
        if hor[i] == 22
            g[i] = 2.0 + 1.0 / (x[i] - 1.0)
            continue
        elseif hor[i] == 32
            g[i] = 1.5 * (2.0 + 1.0 / (x[i] - 1.0))
            continue
        elseif hor[i] == 33
            g[i] = 3.0 + 1.0 / (x[i] - 1.0) + 2.0 / (x[i] - 2.0)
        end
    end
end

function computeLeap(model::Model, p::Path, tr::Float64, tauOpt::TauOpt, tauVar::TauVar)::Tuple{Float64,Bool}
    τ = -Inf

    # ---- index of non-critical reactions
    Jncr = Vector{Int}(undef, 0)
    Jcr = Vector{Int}(undef, 0)

    for j in eachindex(tauVar.l)
        tauVar.l[j] = minimum(p.x ./ abs.(vec(tauVar.neg_nu[j, :])))
        tauVar.l[j] >= tauOpt.nc ? push!(Jncr, j) : push!(Jcr, j)
    end

    Acr = view(p.a, Jcr)
    Ancr = view(p.a, Jncr)

    if length(Jncr) > 0

        # ---- No need to update hor[i] == 1, hor[i] == 2 as is constant
        gUpdate!(model.hor, tauVar.g, p.x)
        μ = zeros(Float64, length(p.x))
        σ2 = zeros(Float64, length(p.x))
        for i in tauVar.irs
            @inbounds for j in Jncr
                μ[i] += model.ν[j, i] * p.a[j]
                σ2[i] += model.ν[j, i]^2 * p.a[j]
            end
        end

        maxTerm = max(maximum(view(tauOpt.ϵ * (p.x ./ tauVar.g), tauVar.irs)), 1.0)
        leftTerm = (maxTerm ./ abs.(μ))
        rightTerm = (maxTerm.^2 ./ σ2)
        τ1 = minimum([leftTerm; rightTerm])
    else
        τ1 = Inf
    end

    switchCondition = (tauOpt.dtf / sum(p.a))::Float64

    # ---- total propensity of critical reactions
    sumAcr = sum(Acr)::Float64

    # ---- Loop until suitable τ that leads to non-negative states
    while true

        # ---- Check if we should switch to direct method
        if τ1 < switchCondition
            return τ, true
        end

        # ---- Time to next critical reaction
        τ2 = randexp()/sumAcr#rand(Exponential(1.0 / sumAcr))#-log(rand()) / sumAcr

        # ---- Case 1
        if τ1 <= τ2
            τ = τ1
            τ > tr && (τ = tr)
            @inbounds for s in eachindex(tauVar.k)
                tauVar.k[s] = zero(eltype(tauVar.k))
            end
            poissonSampler!(Jncr, tauVar.k, Ancr, τ)

        # ---- Case 2
        else
            τ = τ2
            τ > tr && (τ = tr)
            rxn = selectReaction(Acr, sumAcr, length(Acr))::Int
            @inbounds for s in eachindex(tauVar.k)
                tauVar.k[s] = zero(eltype(tauVar.k))
            end
            view(tauVar.k, Jcr)[rxn] = 1
            poissonSampler!(Jncr, tauVar.k, Ancr, τ)
        end

        # ---- Check proposed step, recompute τ if negative species, or break
        for i in eachindex(tauVar.dx)
            tauVar.dx[i] = 0
            for j in eachindex(tauVar.k)
                tauVar.dx[i] += tauVar.k[j] * model.ν[j, i]
            end
        end
        jump = true
        for i in eachindex(p.x)
            if p.x[i] + tauVar.dx[i] < zero(eltype(p.x))
                τ1 /= 2.0
                jump = false
                break
            end
        end
        if jump
            @simd for i in eachindex(p.ra)
                @inbounds p.ra[i] += tauVar.k[i]
            end
            break
        end
    end
    return τ, false
end

function poissonSampler!(Jncr::Vector{Int}, k::Vector{Int}, Ancr::SubArray, τ::Float64)
    for i in eachindex(Ancr)
        view(k, Jncr)[i] = Int(rand(Poisson(Ancr[i] * τ)))
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Printf

function print_progress_header(; detailed_output::Bool)
    @printf(
        "%8s  %9s %9s %13s %13s %13s %9s %8s",
        "time",
        "visited",
        "pending",
        "obj",
        "primal-bound",
        "dual-bound",
        "gap",
        "lp-iter"
    )
    if detailed_output
        @printf(
            " %6s %6s %-24s %6s %6s %6s",
            "node",
            "parent",
            "branch-var",
            "b-val",
            "depth",
            "iinfes"
        )
    end
    println()
    flush(stdout)
end

function print_progress(
    pool::NodePool,
    node::Node;
    time_elapsed::Float64,
    print_interval::Int,
    detailed_output::Bool,
    primal_update::Bool,
)::Nothing
    prefix = primal_update ? "*" : " "
    if (pool.processed % print_interval == 0) || isempty(pool.pending) || primal_update
        @printf(
            "%8.2f %1s%9d %9d %13.6e %13.6e %13.6e %9.2e %8d",
            time_elapsed,
            prefix,
            pool.processed,
            length(pool.processing) + length(pool.pending),
            node.obj * node.mip.sense,
            pool.primal_bound * node.mip.sense,
            pool.dual_bound * node.mip.sense,
            pool.gap,
            pool.mip.lp_iterations,
        )
        if detailed_output
            if isempty(node.branch_variables)
                branch_var_name = "---"
                branch_value = "---"
            else
                branch_var_name = name(node.mip, last(node.branch_variables))
                L = min(24, length(branch_var_name))
                branch_var_name = branch_var_name[1:L]
                branch_value = @sprintf("%.2f", last(node.branch_values))
            end
            @printf(
                " %6d %6s %-24s %6s %6d %6d",
                node.index,
                node.parent === nothing ? "---" : @sprintf("%d", node.parent.index),
                branch_var_name,
                branch_value,
                length(node.branch_variables),
                length(node.fractional_variables)
            )
        end
        println()
        flush(stdout)
    end
end


text:
export  PlaneStressPlasticityLawBased, getStress

mutable struct PlaneStressPlasticityLawBased
    H::Array{Float64} # tangent matrix for plane stress
    E::Float64 # Young's modulus
    ν::Float64 # Poisson's ratio
    ρ::Float64 # density
    a::Array{Float64}  # explicit hardening law parameters
    σ0::Array{Float64} # stress at last time step
    σ0_::Array{Float64} # σ0 to be updated in `commitHistory`
    ε0::Array{Float64} # stress at last time step
    ε0_::Array{Float64} # σ0 to be updated in `commitHistory`
end

@doc """
    The yield function
        \bar{σ} - σY(\bar{\epsilon} ,  a)

    σY is described by a hardening law, for example
    The power-law:
            σY = a[1] + a[2] ε_equiv^a[3]

    To modify the law, you need to modify both function f(σ, ε, a) and fε(ε, a)

""" ->
function f(σ, ε, a)
    # von Mises equivalent strain
    ε_equiv = sqrt(4.0/9.0*(ε[1]^2 + ε[2]^2 - ε[1]*ε[2]) + 1.0/3.0*ε[3]^2)
    σY = a[1] + a[2]*ε_equiv^a[3]
    return sqrt(σ[1]^2-σ[1]*σ[2]+σ[2]^2+3*σ[3]^2)-σY
end

function fσ(σ)
    σ1, σ2, σ3 = σ[1], σ[2], σ[3]
    J2 = sqrt(σ1^2-σ1*σ2+σ2^2+3*σ3^2)
    z = [(σ1 - σ2/2)/J2;
        (-σ1/2 + σ2)/J2;
        3*σ3/J2]
end

function fε(ε, a)
    ε1, ε2, ε3 = ε[1], ε[2], ε[3]
    ε_equiv = sqrt(4.0/9.0*(ε[1]^2 + ε[2]^2 - ε[1]*ε[2]) + 1.0/3.0*ε[3]^2)

    z = [(4.0/9.0 * ε1 - 2.0/9.0 * ε2)/ε_equiv;
         (4.0/9.0 * ε2 - 2.0/9.0 * ε1)/ε_equiv;
         1.0/3.0 * ε3/ε_equiv] * (a[2]*a[3]*ε_equiv^(a[3] - 1.0))

    return z
end

function fσσ(σ)
    σ1, σ2, σ3 = σ[1], σ[2], σ[3]
    J2 = sqrt(σ1^2-σ1*σ2+σ2^2+3*σ3^2)
    [     (-σ1 + σ2/2)*(σ1 - σ2/2)/J2^3 + 1/J2 (σ1/2 - σ2)*(σ1 - σ2/2)/J2^3 - 1/(2*J2)                                   -3*σ3*(σ1 - σ2/2)/J2^3;
    (-σ1 + σ2/2)*(-σ1/2 + σ2)/J2^3 - 1/(2*J2)    (-σ1/2 + σ2)*(σ1/2 - σ2)/J2^3 + 1/J2                                  -3*σ3*(-σ1/2 + σ2)/J2^3;
    3*σ3*(-σ1 + σ2/2)/J2^3                                                        3*σ3*(σ1/2 - σ2)/J2^3 -9*σ3^2/J2^3 + 3/J2]
end

function PlaneStressPlasticityLawBased(prop::Dict{String, Any})
    E = prop["E"]; ν = prop["nu"]; ρ = prop["rho"]; a = prop["args"]
    H = zeros(3,3)
    H[1,1] = E/(1. -ν*ν)
    H[1,2] = H[1,1]*ν
    H[2,1] = H[1,2]
    H[2,2] = H[1,1]
    H[3,3] = E/(2.0*(1.0+ν))
    σ0 = zeros(3); σ0_ = zeros(3); ε0 = zeros(3); ε0_ = zeros(3)
    PlaneStressPlasticityLawBased(H, E, ν, ρ, a, σ0, σ0_, ε0, ε0_)
end

@doc """
    For debugging pathdependent materials, we an build arbitrary pathdependent material law

    sigma = sigma0 + (eps - eps0)**2

""" ->
function getStress(self::PlaneStressPlasticityLawBased,  strain::Array{Float64},  Dstrain::Array{Float64}, Δt::Float64 = 0.0)
    Newton_maxiter = 10
    Newton_Abs_Err = 1e-8
    Newton_Rel_Err = 1e-5

    local dΔσdΔε
    ε = strain
    ε0 = Dstrain
    σ0 = self.σ0

    E = self.E
    H = self.H
    a = self.a

    Δγ = 0.0

    σ = σ0 + H*(ε-ε0)
    r2 = f(σ, ε, a)
    #@show r2, α, K, σY
    if r2 < 0
        #@show "Elasticity"
        σ = σ0 + H*(ε-ε0)
        dΔσdΔε = H
        #@show "elastic", α, α0, σ, σ0,ε, ε0, σ0 + H*(ε-ε0) , H*ε
    else
        # solve for Δγ
        function compute(σ, Δγ)
            r1 = σ - (σ0 + H*(ε - ε0) - Δγ * H* fσ(σ))
            r2 = f(σ, ε, a)
            J = [UniformScaling(1.0)+Δγ*H*fσσ(σ)  H*fσ(σ)
                 reshape(fσ(σ),1,3)                  0]
            return [r1;r2], J
        end

        function compute_sensitivity(σ, Δγ)
            J = [UniformScaling(1.0)+Δγ*H*fσσ(σ)  H*fσ(σ)
                 reshape(fσ(σ),1,3)                  0]

            δ = J\[H ;reshape(fε(ε, a),1,3)]

            return δ[1:3,:]
        end
        res0, _ = compute(σ, Δγ)
        for i = 1:Newton_maxiter
            res, J = compute(σ, Δγ)
            δ = -J\res
            σ += δ[1:3]; Δγ += δ[4]
            # #@show norm(res)/norm(res0)
            if norm(res)/norm(res0) < Newton_Rel_Err || norm(res) < Newton_Abs_Err
                break
            end
            if i == Newton_maxiter
                function f(∂∂u)
                    res, J = compute(∂∂u[1:3],∂∂u[4])
                end
                gradtest(f, [σ;Δγ])
                error("Newton in GetStress does not converge ", res, res0, norm(res)/norm(res0))
            end
        end

        # if f(σ, α, σY, K)>0
        #     @show sqrt(σ[1]^2-σ[1]*σ[2]+σ[2]^2+3*σ[3]^2), α, K
        # end

        dΔσdΔε = compute_sensitivity(σ, Δγ)

    end
    # #@show Δγ
    self.σ0_ = σ[:]
    self.ε0_ = ε
    # self.σ0_ = self.σ0
    #@show σ, dΔσdΔε

    return σ, dΔσdΔε
end

function getTangent(self::PlaneStressPlasticityLawBased)
    error("Not implemented")
end

function commitHistory(self::PlaneStressPlasticityLawBased)
    self.σ0 = self.σ0_
    self.ε0 = self.ε0_
end


text:
# MIT License

# Copyright (c) 2021 Dr. William A. Pisani

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    angles(bonds::Array{Array{Int64,N} where N,1}) -> Array{Tuple{Int64,Int64,Int64},1}

Determines the angles of a molecule or molecular system from the array of bonds given as input.
"""
function angles(bonds::Array{Array{Int64,N} where N,1})  # returns a set of tuples
    println("Determining angles... ")
    angle_set = Set(Tuple{Int64,Int64,Int64}[])
    count = 1
    for i in 1:length(bonds)
        atom1 = bonds[i][1]
        atom2 = bonds[i][2]
        atom12set = Set([atom1,atom2])
        for j in 1:length(bonds)
            atom3 = bonds[j][1]
            atom4 = bonds[j][2]
            atom34set = Set([atom3,atom4])
            atom1234set = Set([atom1,atom2,atom3,atom4])
            if length(atom1234set) == 3
                centeratom = collect(setdiff(atom12set,symdiff(atom12set,atom34set)))
                outsideatoms = collect(symdiff(atom12set,atom34set))
                angle = (outsideatoms[1],centeratom[1],outsideatoms[2])
                if ∉(angle,angle_set)
                    angle = reverse(angle)
                    push!(angle_set,angle)
                    count += 1
                end
            end
        end
    end
    println("Angles determined!")
    return collect(angle_set) # Return Array{Tuple{Int64,Int64,Int64},1}
end

max_stars_repo_name: wapisani/CompMatSciTools

max_issues_repo_name: wapisani/CompMatSciTools

max_forks_repo_name: wapisani/CompMatSciTools


text:
using ImageDraw
using Test

@testset "Polygon" begin
    vert=CartesianIndex{2}[]
    push!(vert, CartesianIndex(1,1))
    push!(vert, CartesianIndex(1,3))
    push!(vert, CartesianIndex(1,5))
    push!(vert, CartesianIndex(3,3))
    push!(vert, CartesianIndex(3,1))

    img = @inferred draw(zeros(Gray{Bool},5,5), Polygon(vert))
    @test all(x->x==true, img[1,:])==true
    @test all(x->x==true, img[1:3,1])==true
    @test all(x->x==true, img[3,1:3])==true
    @test all(x->x==false, img[2,2:3])==true
    @test img[2,4]==true

    img=zeros(Gray{Bool},5,5)
    @inferred draw!(img, Polygon(vert))
    @test all(x->x==true, img[1,:])==true
    @test all(x->x==true, img[1:3,1])==true
    @test all(x->x==true, img[3,1:3])==true
    @test all(x->x==false, img[2,2:3])==true
    @test img[2,4]==true

    img=zeros(RGB{N0f8},5,5)
    @inferred draw!(img, Polygon(vert), RGB(1,0,0))
    @test all(x->x==RGB{N0f8}(1,0,0), img[1,:])==true
    @test all(x->x==RGB{N0f8}(1,0,0), img[1:3,1])==true
    @test all(x->x==RGB{N0f8}(1,0,0), img[3,1:3])==true
    @test all(x->x==RGB{N0f8}(0,0,0), img[2,2:3])==true
    @test img[2,4]==RGB{N0f8}(1,0,0)

    poly_tuples = [(1,1),(3,1),(5,1),(3,3),(1,3)]
    img = @inferred draw(zeros(Gray{Bool},5,5), Polygon(poly_tuples))
    @test all(x->x==true, img[1,:])==true
    @test all(x->x==true, img[1:3,1])==true
    @test all(x->x==true, img[3,1:3])==true
    @test all(x->x==false, img[2,2:3])==true
    @test img[2,4]==true
end

@testset "Path" begin
    vert=CartesianIndex{2}[]
    push!(vert, CartesianIndex(1,1))
    push!(vert, CartesianIndex(1,3))
    push!(vert, CartesianIndex(1,5))
    push!(vert, CartesianIndex(3,3))
    push!(vert, CartesianIndex(3,1))

    img=zeros(RGB{N0f8},5,5)
    @inferred draw!(img, Path(vert))
    @test all(x->x==RGB{N0f8}(1,1,1), img[1,:])==true
    @test all(x->x==RGB{N0f8}(1,1,1), img[3,1:3])==true
    @test all(x->x==RGB{N0f8}(0,0,0), img[2,1:3])==true
    @test img[2,4]==RGB{N0f8}(1,1,1)

    poly_tuples = [(1,1),(3,1),(5,1),(3,3),(1,3)]

    img=zeros(RGB{N0f8},5,5)
    @inferred draw!(img, Path(poly_tuples))
    @test all(x->x==RGB{N0f8}(1,1,1), img[1,:])==true
    @test all(x->x==RGB{N0f8}(1,1,1), img[3,1:3])==true
    @test all(x->x==RGB{N0f8}(0,0,0), img[2,1:3])==true
    @test img[2,4]==RGB{N0f8}(1,1,1)

    # No longer invalid, the program will draw what it can on the image

    # img = zeros(Gray{N0f8},5,5)
    # invalid_points = [(6,1), (4,4), (1,7), (7,6)]
    # @test_throws ErrorException draw!(img, Path(invalid_points))
    #
    # invalid_points = [(4,4), (1,7), (7,6)]
    # @test_throws ErrorException draw!(img, Path(invalid_points))
end

@testset "RegularPolygon" begin
    img = zeros(Gray, 10, 10)
    expected = copy(img)
    expected[3, 2:6] .= Gray(1)
    expected[7, 2:6] .= Gray(1)
    expected[3:7, 2] .= Gray(1)
    expected[3:7, 6] .= Gray(1)
    @test all(expected .== @inferred draw(img, RegularPolygon(CartesianIndex(5,4), 4, 4, π/4))) == true
end


text:
function score(x, y)
    circles = [
        (1, 10), # inner
        (5, 5),  # middle
        (10, 1), # outer
    ]

    # distance from origin = (0, 0)
    d = sqrt(x^2 + y^2)

    for (radius, score) in circles
        if d <= radius
            return score
        end
    end
    return 0
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using ObaraSaika
using Test
using PyCall
using QuadGK

import ObaraSaika: ERI, compute_integrals, boys, compute_normalisation

@testset "ObaraSaika.jl" begin
    @testset "Normalisation" begin
        f(x) = exp(-1.5 * abs(x)^2)
        integral = 4π * quadgk(x -> f(x) * f(x) * x * x, 0, 100)[1]
        @test compute_normalisation(1.5, [0]) ≈ 1 / sqrt(integral)
    end

    @testset "Boys" begin
        boystest(m, t) = quadgk(x -> x^(2m) * exp(-t * x^2), 0, 1)[1]
        @test boystest(0, 1.5)   ≈ boys(0, 1.5)
        @test boystest(0, 0.375) ≈ boys(0, 0.375)
        @test boystest(1, 1.5)   ≈ boys(1, 1.5)
        @test boystest(1, 0.375) ≈ boys(1, 0.375)
    end

    @testset "H 1s" begin
        bas = pyimport("pyscf.gto").parse("H S\n1.5 1.0\nH S\n0.9 1.0")
        mol  = pyimport("pyscf.gto").M(basis=bas, atom="H 0 0 0", spin=1, unit="Bohr")
        ref = mol.intor("int2e")

        #                    c    zeta  P         am
        basis_functions = [[(1.0, 1.5, [0, 0, 0], 0)],
                           [(1.0, 0.9, [0, 0, 0], 0)]]
        res = compute_integrals(ERI(), basis_functions)
        @test res ≈ ref
    end

    @testset "H₂ 1s" begin
        bas = pyimport("pyscf.gto").parse("H S\n1.0 1.0\n")
        mol  = pyimport("pyscf.gto").M(basis=bas, atom="H 0 0 0; H 0 1 0", unit="Bohr")
        ref = mol.intor("int2e")

        #                    c    zeta  P         am
        basis_functions = [[(1.0, 1.0, [0, 0, 0], 0)],
                           [(1.0, 1.0, [0, 1, 0], 0)]]
        res = compute_integrals(ERI(), basis_functions)
        @test res ≈ ref
    end

    @testset "H sto-3g" begin
        mol  = pyimport("pyscf.gto").M(basis="sto-3g", atom="H 0 0 0", spin=1, unit="Bohr")
        ref = mol.intor("int2e")

        #                    c           zeta         P         am
        basis_functions = [[(0.15432897, 3.42525091, [0, 0, 0], 0),
                            (0.53532814, 0.62391373, [0, 0, 0], 0),
                            (0.44463454, 0.16885540, [0, 0, 0], 0)]]
        res = compute_integrals(ERI(), basis_functions)
        @test res ≈ ref atol=1e-7  # TODO That's not a great agreement
    end

    @testset "H₂ sto-3g" begin
        mol  = pyimport("pyscf.gto").M(basis="sto-3g", atom="H 0 0 0; H 0 1 0", unit="Bohr")
        ref = mol.intor("int2e")

        #                    c           zeta         P         am
        basis_functions = [[(0.15432897, 3.42525091, [0, 0, 0], 0),
                            (0.53532814, 0.62391373, [0, 0, 0], 0),
                            (0.44463454, 0.16885540, [0, 0, 0], 0)],
                           [(0.15432897, 3.42525091, [0, 1, 0], 0),
                            (0.53532814, 0.62391373, [0, 1, 0], 0),
                            (0.44463454, 0.16885540, [0, 1, 0], 0)]]
        res = compute_integrals(ERI(), basis_functions)
        @test res ≈ ref atol=1e-7  # TODO That's not a great agreement
    end

    @testset "Be 3-21G" begin
        mol  = pyimport("pyscf.gto").M(basis="3-21g", atom="Be 0 0 0", unit="Bohr")
        ref = mol.intor("int2e")
        #                    c           zeta         P         am
        basis_functions = [[( 0.0644263, 71.8876000, [0, 0, 0], 0),
                            ( 0.3660960, 10.7289000, [0, 0, 0], 0),
                            ( 0.6959340,  2.2220500, [0, 0, 0], 0)],
                           [(-0.4210640,  1.2954800, [0, 0, 0], 0),
                            ( 1.2240700,  0.2688810, [0, 0, 0], 0)],
                           [(-0.4210640,  1.2954800, [0, 0, 0], 1),
                            ( 1.2240700,  0.2688810, [0, 0, 0], 1)],
                           [( 1.0000000,  0.0773500, [0, 0, 0], 0)],
                           [( 1.0000000,  0.0773500, [0, 0, 0], 1)]]
        res = compute_integrals(ERI(), basis_functions)
        @test res ≈ ref atol=1e-7  # TODO That's not a great agreement
    end
end


text:
using JuMP

function TAP(demands, fcoeffs)

    #load OD pair labels
    odPairLabel = readall("od_pair_label_dict_Sioux_simp_refined.json")
    odPairLabel = JSON.parse(odPairLabel)

    odPairLabel_ = readall("od_pair_label_dict__Sioux_simp_refined.json")
    odPairLabel_ = JSON.parse(odPairLabel_)

    #load link labels
    linkLabel = readall("link_label_dict_Sioux_simp.json")
    linkLabel = JSON.parse(linkLabel)

    linkLabel_ = readall("link_label_dict_Sioux_simp_.json")
    linkLabel_ = JSON.parse(linkLabel_)

    #load node-link incidence
    nodeLink = readall("node_link_incidence_Sioux_simp.json")
    nodeLink = JSON.parse(nodeLink)

    start_node = ta_data.start_node

    demandsVec = zeros(length(odPairLabel_))

    for i = 1:length(demandsVec)
        demandsVec[i] = demands[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])]
    end

    # m = Model(solver=GurobiSolver(OutputFlag=false))
    m = Model()

    @defVar(m, linkFlow[1:numLinks])

    @defVar(m, odLinkFlow[keys(demands), 1:numLinks])
    for (s,t) = keys(demands)
        for k = 1:numLinks
            @addConstraint(m, odLinkFlow[(s,t), k] >= 0)
        end
    end

    for k = 1:numLinks
        linkFlowPartial = sum([odLinkFlow[(s,t), k] for (s,t) = keys(demands)])
        @addConstraint(m, linkFlowPartial == linkFlow[k])
    end

    for l = 1:numNodes
        for (s,t) = keys(demands)
            if s != t
                odLinkFlowPartial = sum([nodeLink["$(l-1)-$(k-1)"] * odLinkFlow[(s,t), k] for k = 1:numLinks])
                if (l == s)
                    @addConstraint(m, odLinkFlowPartial + demands[(s,t)] == 0)
                elseif (l == t)
                    @addConstraint(m, odLinkFlowPartial - demands[(s,t)] == 0)
                else
                    @addConstraint(m, odLinkFlowPartial == 0)
                end
            end
        end
    end

    @defNLExpr(g[a=1:numLinks], sum{fcoeffs[i] * (linkFlow[a]^i) / (capacity[a]^(i-1)) / i,
        i = 1:length(fcoeffs)})

    @defNLExpr(f, sum{free_flow_time[a] * g[a], a = 1:numLinks})

#    @defNLExpr(f, sum{free_flow_time[a]*linkFlow[a] + .03*free_flow_time[a]*((linkFlow[a])^5)/((capacity[a])^4),
#        a = 1:numLinks})
#     @defNLExpr(f, sum{free_flow_time[a]*linkFlow[a] + .15*free_flow_time[a]*((linkFlow[a])^2)/(capacity[a]),
#         a = 1:numLinks})

    @setNLObjective(m, Min, f)

    TT = STDOUT # save original STDOUT stream
    redirect_stdout()
    solve(m)
    redirect_stdout(TT) # restore STDOUT

    tapFlowVect = getValue(linkFlow)

    tapFlows = Dict{(Int64,Int64),Float64}()

    for i = 1:length(ta_data.start_node)
        key = (ta_data.start_node[i], ta_data.end_node[i])
        tapFlows[key] = tapFlowVect[i]
    end

#    getObjectiveValue(m)

    return tapFlows, tapFlowVect
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic