text:
"""
    dn = datenum_gregjulian(year,month,day,gregorian::Bool)

Days since 1858-11-17 according to the Gregorian (`gregorian` is `true`) or
Julian Calendar (`gregorian` is `false`) based on the Algorithm of
Jean Meeus [1] modified to handle dates prior to 100 BC.

The year -1, correspond to 1 BC. The year 0 does not exist in the
Gregorian or Julian Calendar.

[1] Meeus, Jean (1998) Astronomical Algorithms (2nd Edition).
Willmann-Bell,  Virginia. p. 63
"""
function datenum_gregjulian(year,month,day,gregorian::Bool,has_year_zero = false)
    # turn year equal to -1 (1 BC) into year = 0
    if (year < 0) && !has_year_zero
        year = year+1
    end

    if gregorian
        # bring year in range of 1601 to 2000
        ncycles = (2000 - year) ÷ 400
        year = year + 400 * ncycles
        # 146_097 = days in 400 years = 365.2425 * 400
        return datenum_ac(year,month,day,gregorian) - ncycles*146_097
    else
        return datenum_ac(year,month,day,gregorian)
    end

end

# Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,  Virginia. p. 63
# However, the algorithm does not work for -100:03:01 and before in
# the proleptic Gregorian Calendar

function datenum_ac(year,month,day,gregorian::Bool)

    if month <= 2
        # if the date is January or February, it is considered
        # the 13rth or 14th month of the preceeding year
        year = year - 1
        month = month + 12
    end

    B =
        if gregorian
            A = year ÷ 100
            2 - A + A ÷ 4
        else
            0
        end

    # benchmark shows that it is 40% faster replacing
    # trunc(Int64,365.25 * (year + 4716))
    # by
    # (1461 * (year + 4716)) ÷ 4
    #
    # and other floating point divisions

    # Z is the Julian Day plus 0.5
    # 1461/4 is 365.25
    # 153/5 is 30.6

    # why 153/5 (or 30.6001 ?)
    # month+1 varies between 4 (March), 5 (April), .. 14 (December),
    # 15 (January), 16 (February)

    # cm = 153 * (4:16) ÷ 5; cm[2:end]-cm[1:end-1]
    #
    # length of each month
    # --------------------
    # 31  March
    # 30  April
    # 31  May
    # 30  June
    # 31  July
    # 31  August
    # 30  September
    # 31  October
    # 30  November
    # 31  December
    # 31  January
    # 30  February (wrong, but not used, since it is the last month)

    Z = (1461 * (year + 4716)) ÷ 4 + (153 * (month+1)) ÷ 5 + day + B - 2401525
    # Modified Julan Day
    return Z + DATENUM_OFFSET
end

"""
    year, month, day = datetuple_gregjulian(Z::Integer,gregorian::Bool)

Compute year, month and day from Z which is the Modified Julian Day
for the Gregorian (true) or Julian (false) calendar.
using the algorithm of Meeus (1998) modified to handle dates prior to 300 AC.

For example:
Z = 0 for the 1858 November 17 00:00:00

Algorithm:

Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,
Virginia. p. 63
"""
function datetuple_gregjulian(Z0::T,gregorian::Bool,has_year_zero = false) where T

    # promote to at least Int64
    Z = promote_type(T, Int64)(Z0)

    # Z is Julian Day plus 0.5
    Z = Z + (2_400_001 - DATENUM_OFFSET)

    A = Z

    if gregorian
        # 1867216.25 - 0.5 corresponds to the date 400-02-29T18:00:00
        # lets magic happen
        #α = trunc(Int64, (Z - 1867_216.25)/36524.25)
        # 400 years = 146097 days = (400 * 365 + 100 - 4 + 1) days
        # α number of centuries since 400-02-29
        α = (4*Z - 7468865) ÷ 146097

        β = Z + (-7468865 + 146097) ÷ 4
        if (β <= 0) && (β % 146097 == 0)
            # correction for 300-03-01 AC, 101-03-01 BC, 501-03-01 BC, ...
            A += 1
        end

        # +α: add leap days for 1700, 1800, 1900, 2000, 2100,
        # -(α ÷ 4): remove leap days for 2000, 2400, ... (already included)
        # so that Julian and Gregorian calendar coincide
        A += 1 + α - (α ÷ 4)
    end

    # even more magic...
    B = A + 1524
    #C = trunc(Int64, (B - 122.1) / 365.25)
    C = (20*B - 2442) ÷ 7305
    #D = trunc(Int64, 365.25 * C)
    # 1461 = 3*365 + 366
    D = 1461 * C ÷ 4
    #E = trunc(Int64, (B-D)/30.6001)
    E = (10000 * (B-D)) ÷ 306001

    #day = B - D - trunc(Int64,30.6001 * E)
    day = B - D - (306001 * E) ÷ 10000

    month = (E < 14 ? E-1 : E-13)
    y = (month > 2 ? C - 4716 : C - 4715)

    # turn year 0 into year -1 (1 BC)
    if (y <= 0) && !has_year_zero
        y = y-1
    end
    return y,month,day
end

datetuple_ymd(::Type{DateTimeProlepticGregorian},Z::Number) =
    datetuple_gregjulian(Z,true,_hasyear0(DateTimeProlepticGregorian))

datetuple_ymd(::Type{DateTimeJulian},Z::Number) =
    datetuple_gregjulian(Z,false,_hasyear0(DateTimeJulian))

datetuple_ymd(::Type{DateTimeStandard},Z::Number) =
    datetuple_gregjulian(Z,Z >= DN_GREGORIAN_CALENDAR,_hasyear0(DateTimeStandard))

datenum(::Type{DateTimeProlepticGregorian},y,m,d) =
    datenum_gregjulian(y,m,d,true,_hasyear0(DateTimeProlepticGregorian))

datenum(::Type{DateTimeJulian},y,m,d) =
    datenum_gregjulian(y,m,d,false,_hasyear0(DateTimeJulian))

datenum(::Type{DateTimeStandard},y,m,d) =
    datenum_gregjulian(y,m,d,(y,m,d) >= GREGORIAN_CALENDAR,_hasyear0(DateTimeStandard))


text:
if isempty(ARGS)
    using Bukdu
    @info :BUKDU_VERSION Bukdu.BUKDU_VERSION
end

using Jive
runtests(@__DIR__, node1=["bukdu/plugs/loggers"])

ext: jl


text:
function start_webserver(host=SERVER_HOST, port=SERVER_PORT)
    @info "Listening on $SERVER_HOST:$SERVER_PORT"
    HTTP.serve(host, port) do request::HTTP.Request
        try
            @debug "Processing request."
            response = handle_http_request_method(request)
            if response != nothing
                return response
            end
            # check and parse raw http request
            check_content_type(request)
            r = parse_content_body(request)
            # parse input parameters (with validation)
            parse_parameters!(r)
            # calculate response
            response = run_ocal(r)
            return build_response(200, response)
        catch e
            if isa(e, ArgumentError)
                return build_response(400, Dict(:error => e.msg))
            else
                return build_response(500, Dict(:error => "Unknown error $e."))
            end
        end
    end
end


text:
struct PosteriorLTISDE{Tprior<:LTISDE, Tdata} <: AbstractGP
    prior::Tprior
    data::Tdata
end

# Avoids method ambiguity.
posterior(fx::FiniteLTISDE, y::AbstractVector) = _posterior(fx, y)
posterior(fx::FiniteLTISDE, y::AbstractVector{<:Real}) = _posterior(fx, y)

_posterior(fx, y) = PosteriorLTISDE(fx.f, (y=y, x=fx.x, Σy=fx.Σy))

const FinitePosteriorLTISDE = FiniteGP{<:PosteriorLTISDE}

function AbstractGPs.cov(fx::FinitePosteriorLTISDE)
    @error "Intentionally not implemented. Please don't try to explicitly compute this cov. matrix."
end

function AbstractGPs.marginals(fx::FinitePosteriorLTISDE)
    if fx.x != fx.f.data.x
        x, y, σ²s, pr_indices = build_inference_data(fx.f, fx.x)

        model = build_lgssm(fx.f.prior(x, σ²s))
        σ²s_pr_full = build_prediction_obs_vars(pr_indices, x, fx.Σy.diag)
        model_post = replace_observation_noise_cov(posterior(model, y), σ²s_pr_full)
        return map(marginals, marginals(model_post)[pr_indices])
    else
        f = Zygote.literal_getfield(fx, Val(:f))
        prior = Zygote.literal_getfield(f, Val(:prior))
        x = Zygote.literal_getfield(fx, Val(:x))
        data = Zygote.literal_getfield(f, Val(:data))
        Σy = Zygote.literal_getfield(data, Val(:Σy))
        Σy_diag = Zygote.literal_getfield(Σy, Val(:diag))
        y = Zygote.literal_getfield(data, Val(:y))

        Σy_new = Zygote.literal_getfield(fx, Val(:Σy))
        Σy_new_diag = Zygote.literal_getfield(Σy_new, Val(:diag))

        model = build_lgssm(AbstractGPs.FiniteGP(prior, x, Σy))
        model_post = replace_observation_noise_cov(posterior(model, y), Σy_new_diag)
        return map(marginals, marginals(model_post))
    end
end

function AbstractGPs.mean_and_var(fx::FinitePosteriorLTISDE)
    ms = marginals(fx)
    return map(mean, ms), map(var, ms)
end

AbstractGPs.mean(fx::FinitePosteriorLTISDE) = mean_and_var(fx)[1]

AbstractGPs.var(fx::FinitePosteriorLTISDE) = mean_and_var(fx)[2]

function AbstractGPs.rand(rng::AbstractRNG, fx::FinitePosteriorLTISDE)
    x, y, σ²s, pr_indices = build_inference_data(fx.f, fx.x)

    model = build_lgssm(fx.f.prior(x, σ²s))
    σ²s_pr_full = build_prediction_obs_vars(pr_indices, x, fx.Σy.diag)
    model_post = replace_observation_noise_cov(posterior(model, y), σ²s_pr_full)
    return rand(rng, model_post)[pr_indices]
end

AbstractGPs.rand(fx::FinitePosteriorLTISDE) = rand(Random.GLOBAL_RNG, fx)

function AbstractGPs.logpdf(fx::FinitePosteriorLTISDE, y_pr::AbstractVector{<:Real})
    x, y, σ²s, pr_indices = build_inference_data(fx.f, fx.x, fx.Σy.diag, y_pr)

    σ²s_pr_full = build_prediction_obs_vars(pr_indices, x, fx.Σy.diag)
    y_pr_full = build_prediction_obs(pr_indices, x, y_pr)

    model = build_lgssm(fx.f.prior(x, σ²s))
    model_post = replace_observation_noise_cov(posterior(model, y), σ²s_pr_full)
    return logpdf(model_post, y_pr_full)
end

# Join the dataset used to construct `f` and the one specified by `x_pred`, `σ²s_pred`, and
# `y_pred`. This is used in all of the inference procedures. Also provide a collection of
# indices that can be used to obtain the requested prediction locations from the joint data.

# This is the most naive way of going about this.
# The present implementation assumes that there are no overlapping data, which will lead to
# numerical issues if violated.
function build_inference_data(
    f::PosteriorLTISDE,
    x_pred::AbstractVector{<:Real},
    σ²s_pred::AbstractVector,
    y_pred::AbstractVector,
)

    # Pull out the input data.
    x_cond = f.data.x
    x_raw = vcat(x_cond, x_pred)

    # Pull out the observations and create arbitrary fake observations at prediction locs.
    y_cond = f.data.y
    y_raw = vcat(y_cond, fill(missing, length(x_pred)))

    # Pull out obs. noise variance and make it really large for prediction locations.
    σ²s_cond = diag(f.data.Σy)
    σ²s_raw = vcat(σ²s_cond, σ²s_pred)

    # Put all of the data in order.
    idx = sortperm(x_raw)
    x = x_raw[idx]
    y = y_raw[idx]
    σ²s = σ²s_raw[idx]

    # The last length(x_pred) indices belong to the predictions.
    pr_indices = sortperm(idx)[end-length(x_pred) + 1:end]

    return x, y, σ²s, pr_indices
end

# If no observations or variances are provided, make the observations arbitrary and the
# variances very large to simulate missing data.
function build_inference_data(f::PosteriorLTISDE, x_pred::AbstractVector{<:Real})
    σ²s_pred = fill(convert(eltype(f.data.Σy), _large_var_const()), length(x_pred))
    y_pred = fill(missing, length(x_pred))
    return build_inference_data(f, x_pred, σ²s_pred, y_pred)
end

# Functions that make predictions at new locations require missings to be placed at the
# locations of the training data.
function build_prediction_obs_vars(
    pr_indices::AbstractVector{<:Integer},
    x_full::AbstractVector{<:Real},
    σ²s_pr::AbstractVector,
)
    σ²s_pr_full = zeros(length(x_full))
    σ²s_pr_full[pr_indices] .= σ²s_pr
    return σ²s_pr_full
end

function build_prediction_obs(
    pr_indices::AbstractVector{<:Integer},
    x_full::AbstractVector{<:Real},
    y_pr::AbstractVector,
)
    y_pr_full = Vector{Union{Missing, eltype(y_pr)}}(undef, length(x_full))
    y_pr_full[pr_indices] .= y_pr
    y_pr_full[setdiff(eachindex(x_full), pr_indices)] .= missing
    return y_pr_full
end


text:
using ARCTR
using Test
using NLPModels
using JuMP

# first test the LDLt function
include("../src/Utilities/testLDLt.jl")

# test all solvers with the well known Woods test function
include("woods.jl")
nlp = MathProgNLPModel(woods(), name="woods")

solver = ALL_solvers[21]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose=true)

reset!(nlp);


text:
get_sources(g::AbstractGraph{T}) where T <: Integer = vertices(g)[indegree(g) .== 0]
get_sinks(g::AbstractGraph{T})  where T <: Integer = vertices(g)[outdegree(g) .== 0]

function add_source_target!(g::AbstractGraph{T}) where T <: Integer
    sinks = get_sinks(g)
    sources = get_sources(g)

    add_vertex!(g)
    for source in sources
        add_edge!(g, vertices(g)[end], source)
    end

    add_vertex!(g)
    for sink in sinks
        add_edge!(g, sink, vertices(g)[end])
    end
end

function N⁻(g::Graphs.SimpleDiGraph{T}, vseqt::Vector{T}) where T <: Integer
    nV = length(vseqt)
    val = zeros(nV)
    val[vseqt[1]] = 1.0
    for v in vseqt
        nb = outneighbors(g, v)
        val[nb] .+= val[v]
    end
    return val
end

function N⁺(g::Graphs.SimpleDiGraph{T}, vseqt::Vector{T}) where T <: Integer
    nV = length(vseqt)
    val = zeros(nV)
    val[vseqt[nV]] = 1.0
    for v in reverse(vseqt)
        nb = inneighbors(g, v)
        val[nb] .+= val[v]
    end
    return val
end

function weights_matrix(g::AbstractGraph{T}, weights::AbstractVector{U}) where T <: Integer where U <: Real
    senders = [src(e) for e in edges(g)]
    receivers = [dst(e) for e in edges(g)]
    return sparse(senders, receivers, weights, nv(g), nv(g))
end

function _weights_spc_raw(g::AbstractGraph{T}) where T <: Integer
    g = Graphs.SimpleDiGraph(g)
    add_source_target!(g)
    vseqt = Graphs.topological_sort_by_dfs(g)
    st = vseqt[[1, end]]
    N_m = N⁻(g, vseqt)
    N_p = N⁺(g, vseqt)
    idx = [!(v in st) for v in vertices(g)]
    rem_vertices!(g, st)
    N_m[idx], N_p[idx], N_p[st[1]]
end

function _weights_spc_vertices(g, normalize)
    N_m, N_p, tf = _weights_spc_raw(g)
    vw = N_m .* N_p

    if normalize == :none
        return vw
    elseif normalize == :log
        return log.(vw)
    elseif normalize == :totalflow
        return vw /= tf
    else
        error("normalize needs to be :none, :log, or :totalflow.")
    end
end

function _weights_spc_edges(g, normalize)
    N_m, N_p, tf = _weights_spc_raw(g)

    ew = map(edges(g)) do e
        N_m[src(e)] * N_p[dst(e)]
    end

    if normalize == :none
        ew = ew
    elseif normalize == :log
        ew = log.(ew)
    elseif normalize == :totalflow
        ew = ew /= tf
    else
        error("normalize needs to be :none, :log, or :totalflow.")
    end

    weights_matrix(g, ew)
end

abstract type MainPathWeight end

abstract type MainPathEdgeWeight <: MainPathWeight end
abstract type MainPathVertexWeight <: MainPathWeight end

"""
    SPCEdge(normalize=:log)(g)

Struct representing Search Path Cout (SPC) edge weights, as defined in Batagelj (2003).
`normalize=:totalflow` indicates that weights should be normalized relative to the
total SPC flow, `normalize=:log` indicates that the logarithm of weights is returned.
Set `normalize` to :none to indicate that weights should not be transformed.

The struct can be called to compute weights directly or
can be passed to the `mainpath` function for dispatch.
"""
Base.@kwdef struct SPCEdge <: MainPathEdgeWeight
    normalize::Symbol = :log
end

(x::SPCEdge)(g) = _weights_spc_edges(g, x.normalize)

"""
    SPCVertex(normalize=:log)(g)

Struct representing Search Path Cout (SPC) vertex weights, as defined in Batagelj (2003).
`normalize=:totalflow` indicates that weights should be normalized relative to the
total SPC flow, `normalize=:log` indicates that the logarithm of weights is returned.
Set `normalize` to :none to indicate that weights should not be transformed.

The struct can be called to compute weights directly or
can be passed to the `mainpath` function for dispatch.
"""
Base.@kwdef struct SPCVertex <: MainPathVertexWeight
    normalize::Symbol = :log
end

(x::SPCVertex)(g) = _weights_spc_vertices(g, x.normalize)


text:
@enum CachingOptimizerState NO_OPTIMIZER EMPTY_OPTIMIZER ATTACHED_OPTIMIZER
@enum CachingOptimizerMode MANUAL AUTOMATIC

# TODO: Benchmark to check if CachingOptimizer should be parameterized on the ModelLike type.

"""
    CachingOptimizer

`CachingOptimizer` is an intermediate layer that stores a cache of the model
and links it with an optimizer. It supports incremental model
construction and modification even when the optimizer doesn't.

A `CachingOptimizer` may be in one of three possible states (`CachingOptimizerState`):

- `NO_OPTIMIZER`: The CachingOptimizer does not have any optimizer.
- `EMPTY_OPTIMIZER`: The CachingOptimizer an empty optimizer. The optimizer is not synchronized with the cached model.
- `ATTACHED_OPTIMIZER`: The CachingOptimizer a optimizer, and it is synchronized with the cached model.

A `CachingOptimizer` has two modes of operation (`CachingOptimizerMode`):

- `MANUAL`: The only methods that change the state of the `CachingOptimizer` are [`reset_optimizer`](@ref), [`drop_optimizer`](@ref), and [`attach_optimizer`](@ref). Attempting to perform an operation in the incorrect state results in an error.
- `AUTOMATIC`: The `CachingOptimizer` changes its state when necessary. For example, `optimize!` will automatically call `attach_optimizer` (an optimizer must have been previously set). Attempting to add a constraint or perform a modification not supported by the optimizer results in a drop to `EMPTY_OPTIMIZER` mode.
"""
mutable struct CachingOptimizer{OptimizerType, ModelType<:MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::Union{Nothing, OptimizerType}
    model_cache::ModelType
    state::CachingOptimizerState
    mode::CachingOptimizerMode
    model_to_optimizer_map::IndexMap
    optimizer_to_model_map::IndexMap
    # CachingOptimizer externally uses the same variable and constraint indices
    # as the model_cache. model_to_optimizer_map maps from the model_cache indices to the
    # optimizer indices.
end

function CachingOptimizer(model_cache::MOI.ModelLike, mode::CachingOptimizerMode)
    CachingOptimizer{MOI.AbstractOptimizer, typeof(model_cache)}(nothing, model_cache, NO_OPTIMIZER, mode, IndexMap(), IndexMap())
end

"""
    CachingOptimizer(model_cache::MOI.ModelLike, optimizer::AbstractOptimizer)

Creates an `CachingOptimizer` in `AUTOMATIC` mode, with the optimizer `optimizer`.
The model_cache manager returned behaves like an `AbstractOptimizer` as long as no
`CachingOptimizer`-specific functions (e.g. `reset_optimizer`) are called on it.
The type of the optimizer returned is `CachingOptimizer{typeof(optimizer),
typeof(model_cache)}` so it does not support the function
`reset_optimizer(::CachingOptimizer, new_optimizer)` if the type of
`new_optimizer` is different from the type of `optimizer`.
"""
function CachingOptimizer(model_cache::MOI.ModelLike, optimizer::MOI.AbstractOptimizer)
    @assert MOI.is_empty(model_cache)
    @assert MOI.is_empty(optimizer)
    CachingOptimizer{typeof(optimizer), typeof(model_cache)}(optimizer, model_cache, ATTACHED_OPTIMIZER, AUTOMATIC, IndexMap(), IndexMap())
end

## Methods for managing the state of CachingOptimizer.

"""
    state(m::CachingOptimizer)::CachingOptimizerState

Returns the state of the CachingOptimizer `m`. See [`CachingOptimizer`](@ref).
"""
state(m::CachingOptimizer) = m.state

"""
    mode(m::CachingOptimizer)::CachingOptimizerMode

Returns the operating mode of the CachingOptimizer `m`. See [`CachingOptimizer`](@ref).
"""
mode(m::CachingOptimizer) = m.mode

"""
    reset_optimizer(m::CachingOptimizer, optimizer::MOI.AbstractOptimizer)

Sets or resets `m` to have the given empty optimizer. Can be called
from any state. The `CachingOptimizer` will be in state `EMPTY_OPTIMIZER` after the call.
"""
function reset_optimizer(m::CachingOptimizer, optimizer::MOI.AbstractOptimizer)
    @assert MOI.is_empty(optimizer)
    m.optimizer = optimizer
    m.state = EMPTY_OPTIMIZER
    return
end

"""
    reset_optimizer(m::CachingOptimizer)

Detaches and empties the current optimizer. Can be called from `ATTACHED_OPTIMIZER`
or `EMPTY_OPTIMIZER` state. The `CachingOptimizer` will be in state `EMPTY_OPTIMIZER`
after the call.
"""
function reset_optimizer(m::CachingOptimizer)
    m.state == EMPTY_OPTIMIZER && return
    @assert m.state == ATTACHED_OPTIMIZER
    MOI.empty!(m.optimizer)
    m.state = EMPTY_OPTIMIZER
    return
end

"""
    drop_optimizer(m::CachingOptimizer)

Drops the optimizer, if one is present. Can be called from any state.
The `CachingOptimizer` will be in state `NO_OPTIMIZER` after the call.
"""
function drop_optimizer(m::CachingOptimizer)
    m.optimizer = nothing
    m.state = NO_OPTIMIZER
    return
end

"""
    attach_optimizer(model::CachingOptimizer)

Attaches the optimizer to `model`, copying all model data into it. Can be called
only from the `EMPTY_OPTIMIZER` state. If the copy succeeds, the
`CachingOptimizer` will be in state `ATTACHED_OPTIMIZER` after the call,
otherwise an error is thrown; see [`copy_to`](@ref) for more details on which
errors can be thrown.
"""
function attach_optimizer(model::CachingOptimizer)
    @assert model.state == EMPTY_OPTIMIZER
    # We do not need to copy names because name-related operations are handled by `m.model_cache`
    indexmap = MOI.copy_to(model.optimizer, model.model_cache, copy_names=false)
    model.state = ATTACHED_OPTIMIZER
    # MOI does not define the type of index_map, so we have to copy it into a
    # concrete container. Also load the reverse map.
    model.model_to_optimizer_map = IndexMap()
    model.optimizer_to_model_map = IndexMap()
    for k in keys(indexmap)
        model.model_to_optimizer_map[k] = indexmap[k]
        model.optimizer_to_model_map[indexmap[k]] = k
    end
end

function MOI.copy_to(m::CachingOptimizer, src::MOI.ModelLike; kws...)
    return automatic_copy_to(m, src; kws...)
end
function supports_default_copy_to(model::CachingOptimizer, copy_names::Bool)
    return supports_default_copy_to(model.model_cache, copy_names)
end

function MOI.empty!(m::CachingOptimizer)
    MOI.empty!(m.model_cache)
    if m.state == ATTACHED_OPTIMIZER
        MOI.empty!(m.optimizer)
    end
    if m.state == EMPTY_OPTIMIZER && m.mode == AUTOMATIC
        m.state = ATTACHED_OPTIMIZER
    end
    m.model_to_optimizer_map = IndexMap()
    m.optimizer_to_model_map = IndexMap()
end
MOI.is_empty(m::CachingOptimizer) = MOI.is_empty(m.model_cache)

# Optimizing and adding/modifying constraints and variables.

function MOI.optimize!(m::CachingOptimizer)
    if m.mode == AUTOMATIC && m.state == EMPTY_OPTIMIZER
        attach_optimizer(m)
    end
    # TODO: better error message if no optimizer is set
    @assert m.state == ATTACHED_OPTIMIZER
    MOI.optimize!(m.optimizer)
end

function MOI.add_variable(m::CachingOptimizer)
    if m.state == ATTACHED_OPTIMIZER
        if m.mode == AUTOMATIC
            try
                vindex_optimizer = MOI.add_variable(m.optimizer)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            vindex_optimizer = MOI.add_variable(m.optimizer)
        end
    end
    vindex = MOI.add_variable(m.model_cache)
    if m.state == ATTACHED_OPTIMIZER
        m.model_to_optimizer_map[vindex] = vindex_optimizer
        m.optimizer_to_model_map[vindex_optimizer] = vindex
    end
    return vindex
end

function MOI.add_variables(m::CachingOptimizer, n)
    if m.state == ATTACHED_OPTIMIZER
        if m.mode == AUTOMATIC
            try
                vindices_optimizer = MOI.add_variables(m.optimizer, n)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            vindices_optimizer = MOI.add_variables(m.optimizer, n)
        end
    end
    vindices = MOI.add_variables(m.model_cache, n)
    if m.state == ATTACHED_OPTIMIZER
        for (vindex, vindex_optimizer) in zip(vindices, vindices_optimizer)
            m.model_to_optimizer_map[vindex] = vindex_optimizer
            m.optimizer_to_model_map[vindex_optimizer] = vindex
        end
    end
    return vindices
end

function MOI.supports_constraint(m::CachingOptimizer, F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    MOI.supports_constraint(m.model_cache, F, S) && (m.state == NO_OPTIMIZER || MOI.supports_constraint(m.optimizer, F, S))
end

function MOI.add_constraint(m::CachingOptimizer, func::MOI.AbstractFunction, set::MOI.AbstractSet)
    if m.state == ATTACHED_OPTIMIZER
        if m.mode == AUTOMATIC
            try
                cindex_optimizer = MOI.add_constraint(m.optimizer, mapvariables(m.model_to_optimizer_map, func), set)
            catch err
                if err isa MOI.NotAllowedError
                    # It could be MOI.AddConstraintNotAllowed{F', S'} with F' != F
                    # or S' != S if, e.g., the `F`-in-`S` constraint is bridged
                    # to other constraints in `m.optimizer`
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            cindex_optimizer = MOI.add_constraint(m.optimizer, mapvariables(m.model_to_optimizer_map, func), set)
        end
    end
    cindex = MOI.add_constraint(m.model_cache, func, set)
    if m.state == ATTACHED_OPTIMIZER
        m.model_to_optimizer_map[cindex] = cindex_optimizer
        m.optimizer_to_model_map[cindex_optimizer] = cindex
    end
    return cindex
end

function MOI.modify(m::CachingOptimizer, cindex::CI, change::MOI.AbstractFunctionModification)
    if m.state == ATTACHED_OPTIMIZER
        cindex_optimizer = m.model_to_optimizer_map[cindex]
        change_optimizer = mapvariables(m.model_to_optimizer_map, change)
        if m.mode == AUTOMATIC
            try
                MOI.modify(m.optimizer, cindex_optimizer, change_optimizer)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.modify(m.optimizer, cindex_optimizer, change_optimizer)
        end
    end
    MOI.modify(m.model_cache, cindex, change)
end

# This function avoids duplicating code in the MOI.set methods for
# ConstraintSet and ConstraintFunction methods, but allows us to strongly type
# the third and fourth arguments of the set methods so that we only support
# setting the same type of set or function.
function replace_constraint_function_or_set(m::CachingOptimizer, attr, cindex, replacement)
    if m.state == ATTACHED_OPTIMIZER
        if m.mode == AUTOMATIC
            try
                MOI.set(m.optimizer, attr, m.model_to_optimizer_map[cindex], replacement)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.set(m.optimizer, attr, m.model_to_optimizer_map[cindex], replacement)
        end
    end
    MOI.set(m.model_cache, attr, cindex, replacement)
end

function MOI.set(m::CachingOptimizer, ::MOI.ConstraintSet, cindex::CI{F,S}, set::S) where {F,S}
    replace_constraint_function_or_set(m, MOI.ConstraintSet(), cindex, set)
end

function MOI.set(m::CachingOptimizer, ::MOI.ConstraintFunction, cindex::CI{F,S}, func::F) where {F,S}
    replace_constraint_function_or_set(m, MOI.ConstraintFunction(), cindex, func)
end

function MOI.modify(m::CachingOptimizer, obj::MOI.ObjectiveFunction, change::MOI.AbstractFunctionModification)
    if m.state == ATTACHED_OPTIMIZER
        change_optimizer = mapvariables(m.model_to_optimizer_map, change)
        if m.mode == AUTOMATIC
            try
                MOI.modify(m.optimizer, obj, change_optimizer)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.modify(m.optimizer, obj, change_optimizer)
        end
    end
    MOI.modify(m.model_cache, obj, change)
end

MOI.is_valid(m::CachingOptimizer, index::MOI.Index) = MOI.is_valid(m.model_cache, index)

function MOI.delete(m::CachingOptimizer, index::MOI.Index)
    if m.state == ATTACHED_OPTIMIZER
        if !MOI.is_valid(m, index)
            # The index thrown by m.model_cache would be xored
            throw(MOI.InvalidIndex(index))
        end
        index_optimizer = m.model_to_optimizer_map[index]
        if m.mode == AUTOMATIC
            try
                MOI.delete(m.optimizer, index_optimizer)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.delete(m.optimizer, index_optimizer)
        end
    end
    # The state may have changed in AUTOMATIC mode since reset_optimizer is
    # called in case the deletion is not supported
    if m.state == ATTACHED_OPTIMIZER
        delete!(m.optimizer_to_model_map, m.model_to_optimizer_map[index])
        delete!(m.model_to_optimizer_map, index)
    end
    MOI.delete(m.model_cache, index)
end

# TODO: add_constraints, transform

## CachingOptimizer get and set attributes

# Attributes are mapped through attribute_value_map (defined in copy.jl) before
# they are sent to the optimizer and when they are returned from the optimizer.
# This map currently only translates indices on MOI.AbstractFunction objects
# between the optimizer indices and the (user-facing) model_cache indices. As a result,
# all MOI.AbstractFunctions must implement mapvariables. Other attributes that
# store indices need to be handled with care.

function MOI.set(m::CachingOptimizer, attr::MOI.AbstractModelAttribute, value)
    if m.state == ATTACHED_OPTIMIZER
        optimizer_value = attribute_value_map(m.model_to_optimizer_map, value)
        if m.mode == AUTOMATIC
            try
                MOI.set(m.optimizer, attr, optimizer_value)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.set(m.optimizer, attr, optimizer_value)
        end
    end
    MOI.set(m.model_cache, attr, value)
end

function MOI.set(m::CachingOptimizer, attr::Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}, index::MOI.Index, value)
    if m.state == ATTACHED_OPTIMIZER
        optimizer_index = m.model_to_optimizer_map[index]
        optimizer_value = attribute_value_map(m.model_to_optimizer_map, value)
        if m.mode == AUTOMATIC
            try
                MOI.set(m.optimizer, attr, optimizer_index, optimizer_value)
            catch err
                if err isa MOI.NotAllowedError
                    reset_optimizer(m)
                else
                    rethrow(err)
                end
            end
        else
            MOI.set(m.optimizer, attr, optimizer_index, optimizer_value)
        end
    end
    MOI.set(m.model_cache, attr, index, value)
end

function MOI.supports(m::CachingOptimizer,
                      attr::Union{MOI.AbstractVariableAttribute,
                                  MOI.AbstractConstraintAttribute},
                      IndexType::Type{<:MOI.Index})
    return MOI.supports(m.model_cache, attr, IndexType) &&
        (m.state == NO_OPTIMIZER || MOI.supports(m.optimizer, attr, IndexType))
end

function MOI.supports(m::CachingOptimizer, attr::MOI.AbstractModelAttribute)
    return MOI.supports(m.model_cache, attr) &&
        (m.state == NO_OPTIMIZER || MOI.supports(m.optimizer, attr))
end

function MOI.get(model::CachingOptimizer, attr::MOI.AbstractModelAttribute)
    if MOI.is_set_by_optimize(attr)
        if state(model) == NO_OPTIMIZER
            if attr == MOI.TerminationStatus()
                return MOI.OPTIMIZE_NOT_CALLED
            elseif attr == MOI.PrimalStatus()
                return MOI.NO_SOLUTION
            elseif attr == MOI.DualStatus()
                return MOI.NO_SOLUTION
            else
                error("Cannot query $(attr) from caching optimizer because no" *
                      " optimizer is attached.")
            end
        end
        return attribute_value_map(model.optimizer_to_model_map,
                                   MOI.get(model.optimizer, attr))
    else
        return MOI.get(model.model_cache, attr)
    end
end
function MOI.get(model::CachingOptimizer,
                 attr::Union{MOI.AbstractVariableAttribute,
                             MOI.AbstractConstraintAttribute},
                 index::MOI.Index)
    if MOI.is_set_by_optimize(attr)
        if state(model) == NO_OPTIMIZER
            error("Cannot query $(attr) from caching optimizer because no " *
                  "optimizer is attached.")
        end
        return attribute_value_map(model.optimizer_to_model_map,
                                   MOI.get(model.optimizer, attr,
                                           model.model_to_optimizer_map[index]))
    else
        return MOI.get(model.model_cache, attr, index)
    end
end
function MOI.get(model::CachingOptimizer,
                 attr::Union{MOI.AbstractVariableAttribute,
                             MOI.AbstractConstraintAttribute},
                 indices::Vector{<:MOI.Index})
    if MOI.is_set_by_optimize(attr)
        if state(model) == NO_OPTIMIZER
            error("Cannot query $(attr) from caching optimizer because no " *
                  "optimizer is attached.")
        end
        return attribute_value_map(model.optimizer_to_model_map,
                                   MOI.get(model.optimizer, attr,
                                           map(index -> model.model_to_optimizer_map[index],
                                               indices)))
  else
      return MOI.get(model.model_cache, attr, indices)
  end
end

#####
##### Names
#####

# Names are not copied, i.e. we use the option `copy_names=false` in
# `attachoptimizer`, so the caching optimizer can support names even if the
# optimizer does not.
function MOI.supports(model::CachingOptimizer,
                      attr::Union{MOI.VariableName,
                                  MOI.ConstraintName},
                      IndexType::Type{<:MOI.Index})
    return MOI.supports(model.model_cache, attr, IndexType)
end
function MOI.set(model::CachingOptimizer,
                 attr::Union{MOI.VariableName, MOI.ConstraintName},
                 index::MOI.Index, value)
    MOI.set(model.model_cache, attr, index, value)
end

function MOI.supports(m::CachingOptimizer,
                      attr::MOI.Name)
    return MOI.supports(m.model_cache, attr)
end
function MOI.set(model::CachingOptimizer, attr::MOI.Name, value)
    MOI.set(model.model_cache, attr, value)
end

function MOI.get(m::CachingOptimizer, IdxT::Type{<:MOI.Index}, name::String)
    return MOI.get(m.model_cache, IdxT, name)
end

# TODO: MOI.set for MOI.AbstractOptimizerAttribute.

function MOI.get(model::CachingOptimizer, attr::MOI.AbstractOptimizerAttribute)
    if state(model) == NO_OPTIMIZER
        error("Cannot query $(attr) from caching optimizer because no " *
              "optimizer is attached.")
    end
    # TODO: Copyable attributes (e.g., TimeLimit) could also be stored in the
    # cache. When MOI.set is implemented for MOI.AbstractOptimizerAttribute,
    # make sure this case is handled correctly.
    return attribute_value_map(model.optimizer_to_model_map,
                               MOI.get(model.optimizer, attr))
end

# Force users to specify whether the attribute should be queried from the
# model_cache or the optimizer. Maybe we could consider a small whitelist of
# attributes to handle automatically.

# These are expert methods to get or set attributes directly in the model_cache
# or optimizer.

struct AttributeFromModelCache{T <: MOI.AnyAttribute}
    attr::T
end

struct AttributeFromOptimizer{T <: MOI.AnyAttribute}
    attr::T
end

function MOI.get(m::CachingOptimizer, attr::AttributeFromModelCache{T}) where {T <: MOI.AbstractModelAttribute}
    return MOI.get(m.model_cache, attr.attr)
end

function MOI.get(m::CachingOptimizer, attr::AttributeFromModelCache{T}, idx) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    return MOI.get(m.model_cache, attr.attr, idx)
end

function MOI.get(m::CachingOptimizer, attr::AttributeFromOptimizer{T}) where {T <: MOI.AbstractModelAttribute}
    @assert m.state == ATTACHED_OPTIMIZER
    return attribute_value_map(m.optimizer_to_model_map,MOI.get(m.optimizer, attr.attr))
end

function MOI.get(m::CachingOptimizer, attr::AttributeFromOptimizer{T}, idx::MOI.Index) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    @assert m.state == ATTACHED_OPTIMIZER
    return attribute_value_map(m.optimizer_to_model_map,MOI.get(m.optimizer, attr.attr, m.model_to_optimizer_map[idx]))
end

function MOI.get(m::CachingOptimizer, attr::AttributeFromOptimizer{T}, idx::Vector{<:MOI.Index}) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    @assert m.state == ATTACHED_OPTIMIZER
    return attribute_value_map(m.optimizer_to_model_map,MOI.get(m.optimizer, attr.attr, getindex.(m.model_to_optimizer_map,idx)))
end

function MOI.set(m::CachingOptimizer, attr::AttributeFromModelCache{T}, v) where {T <: MOI.AbstractModelAttribute}
    return MOI.set(m.model_cache, attr.attr, v)
end

function MOI.set(m::CachingOptimizer, attr::AttributeFromModelCache{T}, idx, v) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    return MOI.set(m.model_cache, attr.attr, idx, v)
end

function MOI.set(m::CachingOptimizer, attr::AttributeFromOptimizer{T}, v) where {T <: MOI.AbstractModelAttribute}
    @assert m.state == ATTACHED_OPTIMIZER
    return MOI.set(m.optimizer, attr.attr, attribute_value_map(m.model_to_optimizer_map,v))
end

# Map vector of indices into vector of indices or one index into one index
map_indices_to_optimizer(m::CachingOptimizer, idx::MOI.Index) = m.model_to_optimizer_map[idx]
map_indices_to_optimizer(m::CachingOptimizer, indices::Vector{<:MOI.Index}) = getindex.(Ref(m.model_to_optimizer_map), indices)
function MOI.set(m::CachingOptimizer, attr::AttributeFromOptimizer{T}, idx, v) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    @assert m.state == ATTACHED_OPTIMIZER
    return MOI.set(m.optimizer, attr.attr, map_indices_to_optimizer(m, idx), attribute_value_map(m.model_to_optimizer_map,v))
end

function MOI.supports(m::CachingOptimizer, attr::AttributeFromModelCache{T}) where {T <: MOI.AbstractModelAttribute}
    return MOI.supports(m.model_cache, attr.attr)
end

function MOI.supports(m::CachingOptimizer, attr::AttributeFromModelCache{T}, idxtype::Type{<:MOI.Index}) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    return MOI.supports(m.model_cache, attr.attr, idxtype)
end

function MOI.supports(m::CachingOptimizer, attr::AttributeFromOptimizer{T}) where {T <: MOI.AbstractModelAttribute}
    @assert m.state == ATTACHED_OPTIMIZER
    return MOI.supports(m.optimizer, attr.attr)
end

function MOI.supports(m::CachingOptimizer, attr::AttributeFromOptimizer{T}, idxtype::Type{<:MOI.Index}) where {T <: Union{MOI.AbstractVariableAttribute,MOI.AbstractConstraintAttribute}}
    @assert m.state == ATTACHED_OPTIMIZER
    return MOI.supports(m.optimizer, attr.attr, idxtype)
end

# TODO: get and set methods to look up/set name strings


text:
"""
Class for Gaussian Processes models

```julia
GP(X::AbstractArray{T1,N1},y::AbstractArray{T2,N2},
   kernel::Union{Kernel,AbstractVector{<:Kernel}};
   noise::Real=1e-5, verbose::Integer=0,
   Autotuning::Bool=true,atfrequency::Integer=1,
   IndependentPriors::Bool=true,ArrayType::UnionAll=Vector)
```

Argument list :

**Mandatory arguments**

 - `X` : input features, should be a matrix N×D where N is the number of observation and D the number of dimension
 - `y` : input labels, can be either a vector of labels for multiclass and single output or a matrix for multi-outputs (note that only one likelihood can be applied)
 - `kernel` : covariance function, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models

**Keyword arguments**

 - `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
 - `Autotuning` : Flag for optimizing hyperparameters
 - `atfrequency` : Choose how many variational parameters iterations are between hyperparameters optimization
 - `IndependentPriors` : Flag for setting independent or shared parameters among latent GPs
 - `ArrayType` : Option for using different type of array for storage (allow for GPU usage)
"""
mutable struct GP{L<:Likelihood,I<:Inference,T<:Real,V<:AbstractVector{T}} <: AbstractGP{L,I,T,V}
    X::Matrix{T} #Feature vectors
    y::LatentArray #Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSample::Int64 # Number of data points
    nDim::Int64 # Number of covariates per data point
    nFeature::Int64 # Number of features of the GP (equal to number of points)
    nLatent::Int64 # Number pf latent GPs
    IndependentPriors::Bool # Use of separate priors for each latent GP
    nPrior::Int64 # Equal to 1 or nLatent given IndependentPriors
    μ₀::LatentArray{MeanPrior{T}}
    Knn::LatentArray{Symmetric{T,Matrix{T}}}
    invKnn::LatentArray{Symmetric{T,Matrix{T}}}
    kernel::LatentArray{Kernel{T}}
    likelihood::Likelihood{T}
    inference::Inference{T}
    verbose::Int64 #Level of printing information
    Autotuning::Bool
    atfrequency::Int64
    Trained::Bool
end

function GP(X::AbstractArray{T1,N1},y::AbstractArray{T2,N2},kernel::Union{Kernel,AbstractVector{<:Kernel}};  noise::Real=1e-5,
            verbose::Integer=0,Autotuning::Bool=true,atfrequency::Integer=1,mean::Union{<:Real,AbstractVector{<:Real},MeanPrior}=ZeroMean(),
            IndependentPriors::Bool=true,ArrayType::UnionAll=Vector) where {T1<:Real,T2,N1,N2}
            likelihood = GaussianLikelihood(noise)
            inference = Analytic()
            X,y,nLatent,likelihood = check_data!(X,y,likelihood)

            nPrior = IndependentPriors ? nLatent : 1
            nFeature = nSample = size(X,1); nDim = size(X,2);
            kernel = ArrayType([deepcopy(kernel) for _ in 1:nPrior])

            Knn = LatentArray([Symmetric(Matrix{T1}(I,nFeature,nFeature)) for _ in 1:nPrior]);
            invKnn = copy(Knn)
            μ₀ = []
            if typeof(mean) <: Real
                μ₀ = [ConstantMean(mean) for _ in 1:nPrior]
            elseif typeof(mean) <: AbstractVector{<:Real}
                μ₀ = [EmpiricalMean(mean) for _ in 1:nPrior]
            else
                μ₀ = [mean for _ in 1:nPrior]
            end
            likelihood = init_likelihood(likelihood,inference,nLatent,nSample)
            inference = init_inference(inference,nLatent,nSample,nSample,nSample)

            model = GP{GaussianLikelihood{T1},Analytic{T1},T1,ArrayType{T1}}(X,y,
                    nFeature, nDim, nFeature, nLatent,
                    IndependentPriors,nPrior,
                    μ₀,Knn,invKnn,kernel,likelihood,inference,
                    verbose,Autotuning,atfrequency,false)
            computeMatrices!(model)
            model.Trained = true
            return model
end

function Base.show(io::IO,model::GP{<:Likelihood,<:Inference,T}) where T
    print(io,"Gaussian Process with a $(model.likelihood) infered by $(model.inference) ")
end


text:
mutable struct IfBlock
    condition::String # condition for code to run
    true_code::String # code that runs if true
end

mutable struct IfElseBlock
    condition::String # condition for code to run
    true_code::String # code that runs if true
end

mutable struct ElseBlock
    true_code::String
end

mutable struct LogicBlock

    # condition and code for if statement
    ifBlock::IfBlock

    nestedifBlock::LogicBlock # for code in the if statement

    # array for if elses
    ifElseBlock::Vector{IfElseBlock}

    nestedIfElseBlock::Vector{LogicBlock} # for code in if else statements

    # code for the else block
    else_code::String

    nestedElseBlock::LogicBlock
end

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
# =================================================================================================
# =                                 High level interface                                          =
# =================================================================================================

# Testing higher level interface

#using Psychro

Tk = 183.3:10.0:463.3
Pa = [0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 49.0]*1e5

# Test volume - saturated moist air
xvs(Tk, P) = Psychro.efactor(Tk,P)*Psychro.Pws(Tk)/P
for p1 in Pa, t1 in Tk
    x = xvs(t1, p1)
    if x <= 1.0
        v1 = volmoist(t1, p1)
        v2 = Psychro.volume(Psychro.MoistAir, t1, Psychro.RelHum, 1.0, p1)
        @test v1 ≈ v2 rtol=1e-6
    end
end

# Test volume - dry air

for p1 in Pa, t1 in Tk
    v1 = Psychro.volumeair(t1, p1)
    v2 = Psychro.volume(Psychro.MoistAir, t1, Psychro.RelHum, 0.0, p1)
    v3 = Psychro.volume(Psychro.DryAir, t1, p1)
    @test v1 ≈ v2 rtol=1e-6
    @test v2 ≈ v3 rtol=1e-6
end

# Test enthalpy - saturated moist air
for p1 in Pa, t1 in Tk
    x = xvs(t1, p1)
    if x <= 1.0
        v1 = hmoist(t1, p1)
        v2 = Psychro.enthalpy(Psychro.MoistAir, t1, Psychro.RelHum, 1.0, p1)/1000
        @test v1 ≈ v2 rtol=1e-4
    end
end

# Test enthalpy - dry air

for p1 in Pa, t1 in Tk
    v1 = Psychro.enthalpyair(t1, p1)
    v2 = Psychro.enthalpy(Psychro.DryAir, t1, p1)
    v3 = Psychro.enthalpy(Psychro.MoistAir, t1, Psychro.RelHum, 0.0, p1)
    #println(p1, " - ", t1, " - ", v1, " - ", v2, " - ", v3)
    @test v1 ≈ v2 rtol=1e-4
    @test v2 ≈ v3 rtol=1e-4
end

# Test entropy - saturated moist air
for p1 in Pa, t1 in Tk
    x = xvs(t1, p1)
    if x <= 1.0
        v1 = smoist(t1, p1)
        v2 = Psychro.entropy(Psychro.MoistAir, t1, Psychro.RelHum, 1.0, p1)/1000
        @test v1 ≈ v2 rtol=1e-6
    end
end

# Test entropy - dry air

for p1 in Pa, t1 in Tk
    v1 = Psychro.entropyair(t1, p1)
    v2 = Psychro.entropy(Psychro.DryAir, t1, p1)
    v3 = Psychro.entropy(Psychro.MoistAir, t1, Psychro.RelHum, 0.0, p1)
    #println(p1, " - ", t1, " - ", v1, " - ", v2, " - ", v3)
    @test v1 ≈ v2 rtol=1e-4
    @test v2 ≈ v3 rtol=1e-4
end

# Calculates M* in table 13 of reference [5]
function mfun(Tbu, P=101325.0)
    hs = Psychro.enthalpy(Psychro.MoistAir, Tbu, Psychro.RelHum, 1.0, P)
    ws = Psychro.humrat(Psychro.MoistAir, Tbu,  Psychro.RelHum, 1.0, P)
    hw = Psychro.enthalpywi(Tbu)
    return (hs - ws*hw)/1000
end

@test 12.94 ≈ mfun(2.0+273.15) atol=0.015
@test 20.50 ≈ mfun(6.0+273.15) atol=0.015
@test 29.03 ≈ mfun(10.0+273.15) atol=0.015
@test 38.78 ≈ mfun(14.0+273.15) atol=0.015
@test 50.03 ≈ mfun(18.0+273.15) atol=0.015
@test 63.11 ≈ mfun(22.0+273.15) atol=0.015
@test 78.46 ≈ mfun(26.0+273.15) atol=0.015
@test 96.57 ≈ mfun(30.0+273.15) atol=0.015
@test 118.07 ≈ mfun(34.0+273.15) atol=0.015
@test 143.74 ≈ mfun(38.0+273.15) atol=0.015
@test 174.58 ≈ mfun(42.0+273.15) atol=0.015
@test 192.31 ≈ mfun(44.0+273.15) atol=0.015
@test 233.36 ≈ mfun(48.0+273.15) atol=0.015

function wfun(t, tbu, P=101325.0)

    ms = mfun(tbu, P)
    ha = Psychro.enthalpyair(t, P)
    hg = Psychro.enthalpyvapor(t)
    hw = Psychro.enthalpywi(tbu)

    return (ms - ha) / (hg - hw)
end

# Testing Table 20 from reference [5].

Tk = vcat(5*ones(4), 25*ones(4), 50*ones(4)) .+ 273.15
Tbu = [5.0, 2.0, -1.0, -3.0, 25.0, 20.0, 15.0, 10.0, 25.0, 22.0, 20.0, 19.0] .+ 273.15
Td = [5.0, -2.16, -11.92, -37.23, 25.0, 17.60, 7.73, -10.42, 13.47, 4.16, -5.48, -14.12]
w = [5.42, 3.16, 1.35, 0.11, 20.17, 12.66, 6.56, 1.55, 9.67, 5.11, 2.39, 1.11]*1e-3
rel = [100.0, 58.6, 25.1, 2.0, 100, 63.5, 33.2, 7.9, 12.5, 6.7, 3.1, 1.5]
vol = [0.794, 0.792, 0.789, 0.788, 0.872, 0.862, 0.853, 0.846, 0.930, 0.923, 0.919, 0.917]
h = [18.64, 12.97, 8.42, 5.30, 76.50, 57.38, 41.85, 29.09, 75.40, 63.58, 56.51, 53.20]
s = [0.0697, 0.0490, 0.0319, 0.0195, 0.2698, 0.2048, 0.1506, 0.1039, 0.2610, 0.2192, 0.1934, 0.1808]

P = 101325.0
Td1 = Psychro.dewpoint.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P) .- 273.15
w1 =  Psychro.humrat.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P)
rel1 =  Psychro.relhum.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P) * 100
vol1 =  Psychro.volume.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P)
h1 =  Psychro.enthalpy.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P)/1000
s1 =  Psychro.entropy.(Psychro.MoistAir, Tk, Psychro.WetBulb, Tbu, P)/1000

@test maximum(abs, Td1-Td) ≈ 0.0 atol=0.01
@test maximum(abs, w1-w) ≈ 0.0 atol=0.00001
@test maximum(abs, rel1-rel) ≈ 0.0 atol=0.1
@test maximum(abs, vol1-vol) ≈ 0.0 atol=0.001
@test maximum(abs, h1-h) ≈ 0.0 atol=0.01
@test maximum(abs, s1-s) ≈ 0.0 atol=0.0001

Tbu1 = Psychro.wetbulb.(Psychro.MoistAir, Tk, Psychro.DewPoint, Td1 .+ 273.15, P)
@test maximum(abs, Tbu1-Tbu) ≈ 0.0 atol=1e-4
Tbu2 = Psychro.wetbulb.(Psychro.MoistAir, Tk, Psychro.RelHum, rel1/100, P)
@test maximum(abs, Tbu2-Tbu) ≈ 0.0 atol=1e-4
Tbu3 = Psychro.wetbulb.(Psychro.MoistAir, Tk, Psychro.HumRat, w1/100, P)
@test maximum(abs, Tbu2-Tbu) ≈ 0.0 atol=1e-4

Td2 = Psychro.dewpoint.(Psychro.MoistAir, Tk, Psychro.HumRat, w1, P) .- 273.15
@test maximum(abs, Td2-Td1) ≈ 0.0 atol=1e-4
Td3 = Psychro.dewpoint.(Psychro.MoistAir, Tk, Psychro.RelHum, rel1/100, P) .- 273.15
@test maximum(abs, Td3-Td1) ≈ 0.0 atol=1e-4


text:
module FileCache

  export maybeload

  type CachedFile
    mtime::Float64
    content::Any
  end

  const cache = Dict{String, Any}()

  function maybeload(f::Function, path::AbstractString)
    p = realpath(path)
    if !(haskey(cache, p) && mtime(p) == cache[p].mtime)
      cache[p] = CachedFile(mtime(p), f(p))
    end
    cache[p].content
  end

end # module


text:
module Quadmath
using Requires

export Float128, ComplexF128

import Base: (*), +, -, /,  <, <=, ==, ^, convert,
          reinterpret, sign_mask, exponent_mask, exponent_one, exponent_half,
          significand_mask, exponent, significand,
          promote_rule, widen,
          string, print, show, parse,
          acos, acosh, asin, asinh, atan, atanh, cosh, cos,
          exp, expm1, log, log2, log10, log1p, sin, sinh, sqrt,
          tan, tanh,
          ceil, floor, trunc, round, fma,
          copysign, flipsign, max, min, hypot, abs,
          ldexp, frexp, modf, nextfloat, eps,
          isinf, isnan, isfinite, isinteger,
          floatmin, floatmax, precision, signbit,
          Int32, Int64, Float64, BigFloat, BigInt

if Sys.isapple()
    const quadoplib = "libquadmath.0"
    const libquadmath = "libquadmath.0"
elseif Sys.isunix()
    const quadoplib = "libgcc_s.so.1"
    const libquadmath = "libquadmath.so.0"
elseif Sys.iswindows()
    const quadoplib = "libgcc_s_seh-1.dll"
    const libquadmath = "libquadmath-0.dll"
end

@static if Sys.isunix()
    # we use this slightly cumbersome definition to ensure that the value is passed
    # on the xmm registers, matching the x86_64 ABI for __float128.
    const Cfloat128 = NTuple{2,VecElement{Float64}}

    struct Float128 <: AbstractFloat
        data::Cfloat128
    end
    convert(::Type{Float128}, x::Number) = Float128(x)

    const ComplexF128 = Complex{Float128}

    Base.cconvert(::Type{Cfloat128}, x::Float128) = x.data

    # reinterpret
    function reinterpret(::Type{UInt128}, x::Float128)
        hi = reinterpret(UInt64, x.data[2].value)
        lo = reinterpret(UInt64, x.data[1].value)
        UInt128(hi) << 64 | lo
    end
    function reinterpret(::Type{Float128}, x::UInt128)
        fhi = reinterpret(Float64, (x >> 64) % UInt64)
        flo = reinterpret(Float64, x % UInt64)
        Float128((VecElement(flo), VecElement(fhi)))
    end
    reinterpret(::Type{Unsigned}, x::Float128) = reinterpret(UInt128, x)
    reinterpret(::Type{Signed}, x::Float128) = reinterpret(Int128, x)

    reinterpret(::Type{Int128}, x::Float128) =
        reinterpret(Int128, reinterpret(UInt128, x))
    reinterpret(::Type{Float128}, x::Int128) =
        reinterpret(Float128, reinterpret(UInt128, x))

elseif Sys.iswindows()
    primitive type Float128 <: AbstractFloat 128
    end
    const Cfloat128 = Float128
end

function __init__()
    @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" begin
        import .SpecialFunctions

        SpecialFunctions.erf(x::Float128) = Float128(ccall((:erfq, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.erfc(x::Float128) = Float128(ccall((:erfcq, libquadmath), Cfloat128, (Cfloat128, ), x))

        SpecialFunctions.besselj0(x::Float128) = Float128(ccall((:j0q, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.besselj1(x::Float128) = Float128(ccall((:j1q, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.bessely0(x::Float128) = Float128(ccall((:y0q, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.bessely1(x::Float128) = Float128(ccall((:y1q, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.besselj(n::Cint, x::Float128) = Float128(ccall((:jnq, libquadmath), Cfloat128, (Cint, Cfloat128), n, x))
        SpecialFunctions.bessely(n::Cint, x::Float128) = Float128(ccall((:ynq, libquadmath), Cfloat128, (Cint, Cfloat128), n, x))

        SpecialFunctions.gamma(x::Float128) = Float128(ccall((:tgammaq, libquadmath), Cfloat128, (Cfloat128, ), x))
        SpecialFunctions.lgamma(x::Float128) = Float128(ccall((:lgammaq, libquadmath), Cfloat128, (Cfloat128, ), x))
    end
end

sign_mask(::Type{Float128}) =        0x8000_0000_0000_0000_0000_0000_0000_0000
exponent_mask(::Type{Float128}) =    0x7fff_0000_0000_0000_0000_0000_0000_0000
exponent_one(::Type{Float128}) =     0x3fff_0000_0000_0000_0000_0000_0000_0000
exponent_half(::Type{Float128}) =    0x3ffe_0000_0000_0000_0000_0000_0000_0000
significand_mask(::Type{Float128}) = 0x0000_ffff_ffff_ffff_ffff_ffff_ffff_ffff

fpinttype(::Type{Float128}) = UInt128

# conversion
Float128(x::Float128) = x

## Float64
Float128(x::Float64) =
    Float128(ccall((:__extenddftf2, quadoplib), Cfloat128, (Cdouble,), x))
Float64(x::Float128) =
    ccall((:__trunctfdf2, quadoplib), Cdouble, (Cfloat128,), x)

Int32(x::Float128) =
    ccall((:__fixtfsi, quadoplib), Int32, (Cfloat128,), x)
Float128(x::Int32) =
    Float128(ccall((:__floatsitf, quadoplib), Cfloat128, (Int32,), x))

Float128(x::UInt32) =
    Float128(ccall((:__floatunsitf, quadoplib), Cfloat128, (UInt32,), x))

Int64(x::Float128) =
    ccall((:__fixtfdi, quadoplib), Int64, (Cfloat128,), x)
Float128(x::Int64) =
    Float128(ccall((:__floatditf, quadoplib), Cfloat128, (Int64,), x))

Float128(x::Rational{T}) where T = Float128(numerator(x))/Float128(denominator(x))

# comparison

(==)(x::Float128, y::Float128) =
    ccall((:__eqtf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) == 0

(<)(x::Float128, y::Float128) =
    ccall((:__letf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) == -1

(<=)(x::Float128, y::Float128) =
    ccall((:__letf2,quadoplib), Cint, (Cfloat128,Cfloat128), x, y) <= 0

# arithmetic

(+)(x::Float128, y::Float128) =
    Float128(ccall((:__addtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(-)(x::Float128, y::Float128) =
    Float128(ccall((:__subtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(*)(x::Float128, y::Float128) =
    Float128(ccall((:__multf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(/)(x::Float128, y::Float128) =
    Float128(ccall((:__divtf3,quadoplib), Cfloat128, (Cfloat128,Cfloat128), x, y))
(-)(x::Float128) =
    Float128(ccall((:__negtf2,quadoplib), Cfloat128, (Cfloat128,), x))
(^)(x::Float128, y::Float128) =
    Float128(ccall((:powq, libquadmath), Cfloat128, (Cfloat128,Cfloat128), x, y))
# math

## one argument
for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :cosh, :cos,
          :erf, :erfc, :exp, :expm1, :log, :log2, :log10, :log1p,
          :sin, :sinh, :sqrt, :tan, :tanh,
          :ceil, :floor, :trunc, )
    @eval function $f(x::Float128)
        Float128(ccall(($(string(f,:q)), libquadmath), Cfloat128, (Cfloat128, ), x))
    end
end
for (f,fc) in (:abs => :fabs,
               :round => :rint,)
    @eval function $f(x::Float128)
        Float128(ccall(($(string(fc,:q)), libquadmath), Cfloat128, (Cfloat128, ), x))
    end
end

## two argument
for f in (:copysign, :hypot, )
    @eval function $f(x::Float128, y::Float128)
       Float128(ccall(($(string(f,:q)), libquadmath), Cfloat128, (Cfloat128, Cfloat128), x, y))
    end
end

flipsign(x::Float128, y::Float128) = signbit(y) ? -x : x

function atan(x::Float128, y::Float128)
    Float128(ccall((:atan2q, libquadmath), Cfloat128, (Cfloat128, Cfloat128), x, y))
end

## misc
fma(x::Float128, y::Float128, z::Float128) =
    Float128(ccall((:fmaq,libquadmath), Cfloat128, (Cfloat128, Cfloat128, Cfloat128), x, y, z))

isnan(x::Float128) =
    0 != ccall((:isnanq,libquadmath), Cint, (Cfloat128, ), x)
isinf(x::Float128) =
    0 != ccall((:isinfq,libquadmath), Cint, (Cfloat128, ), x)
isfinite(x::Float128) =
    0 != ccall((:finiteq,libquadmath), Cint, (Cfloat128, ), x)
isinteger(x::Float128) = isfinite(x) && x === trunc(x)

signbit(x::Float128) = signbit(reinterpret(Int128, x))
precision(::Type{Float128}) = 113

eps(::Type{Float128}) = reinterpret(Float128, 0x3f8f_0000_0000_0000_0000_0000_0000_0000)
floatmin(::Type{Float128}) = reinterpret(Float128, 0x0001_0000_0000_0000_0000_0000_0000_0000)
floatmax(::Type{Float128}) = reinterpret(Float128, 0x7ffe_ffff_ffff_ffff_ffff_ffff_ffff_ffff)

ldexp(x::Float128, n::Cint) =
    Float128(ccall((:ldexpq, libquadmath), Cfloat128, (Cfloat128, Cint), x, n))
ldexp(x::Float128, n::Integer) =
    ldexp(x, clamp(n, typemin(Cint), typemax(Cint)) % Cint)

function frexp(x::Float128)
    r = Ref{Cint}()
    y = Float128(ccall((:frexpq, libquadmath), Cfloat128, (Cfloat128, Ptr{Cint}), x, r))
    return y, Int(r[])
end

function modf(x::Float128)
    isinf(x) && return (zero(Float128), x)
    ipart = trunc(x)
    fpart = x - ipart
    return fpart, ipart
end

significand(x::Float128) = frexp(x)[1] * 2
function exponent(x::Float128)
     !isfinite(x) && throw(DomainError("Cannot be NaN or Inf."))
     abs(x) > 0 && return frexp(x)[2] - 1
     throw(DomainError("Cannot be subnormal converted to 0."))
end

function nextfloat(f::Float128, d::Integer)
    F = typeof(f)
    fumax = reinterpret(Unsigned, F(Inf))
    U = typeof(fumax)

    isnan(f) && return f
    fi = reinterpret(Signed, f)
    fneg = fi < 0
    fu = unsigned(fi & typemax(fi))

    dneg = d < 0
    da = Base.uabs(d)
    if da > typemax(U)
        fneg = dneg
        fu = fumax
    else
        du = da % U
        if fneg ⊻ dneg
            if du > fu
                fu = min(fumax, du - fu)
                fneg = !fneg
            else
                fu = fu - du
            end
        else
            if fumax - fu < du
                fu = fumax
            else
                fu = fu + du
            end
        end
    end
    if fneg
        fu |= sign_mask(F)
    end
    reinterpret(F, fu)
end

Float128(::Irrational{:π}) =  reinterpret(Float128, 0x4000921fb54442d18469898cc51701b8)
Float128(::Irrational{:e}) =  reinterpret(Float128, 0x40005bf0a8b1457695355fb8ac404e7a)

function BigFloat(x::Float128; precision=precision(BigFloat))
    if !isfinite(x) || iszero(x)
        @static if VERSION < v"1.1"
            return BigFloat(Float64(x), precision)
        else
            return BigFloat(Float64(x), precision=precision)
        end
    end

    @static if VERSION < v"1.1"
        b = setprecision(BigFloat, max(precision,113)) do
            BigFloat()
        end
    else
        b = BigFloat(precision=max(precision,113))
    end

    y, k = frexp(x)
    b.exp = Clong(k)
    b.sign = signbit(x) ? Cint(-1) : Cint(1)
    u = (reinterpret(UInt128, y) << 15) | 0x8000_0000_0000_0000_0000_0000_0000_0000
    i = cld(precision, sizeof(Culong)*8)
    while u != 0
        w = (u >> (128-sizeof(Culong)*8)) % Culong
        unsafe_store!(b.d, w, i)
        i -= 1
        u <<= sizeof(Culong)*8
    end
    # set remaining bits to zero
    while i > 0
        unsafe_store!(b.d, zero(Culong), i)
        i -= 1
    end

    if precision < 113
        @static if VERSION < v"1.1"
            b2 = setprecision(BigFloat, precision) do
                BigFloat()
            end
            ccall((:mpfr_set, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32),
                  b2, b, Base.MPFR.ROUNDING_MODE[])
            return b2
        else
            return BigFloat(b, precision=precision)
        end
    else
        return b
    end
end

function Float128(x::BigFloat)
    if !isfinite(x) || iszero(x)
        return Float128(Float64(x))
    end

    y,k = frexp(x)
    if k >= -16381
        prec = 113
    elseif k >= -16381-112
        prec = 113 + (k + 16381)
    elseif k == -16381-113 && abs(y) > 0.5
        z = reinterepret(Float128, UInt128(1))
    else
        z = reinterepret(Float128, UInt128(0))
    end

    @static if VERSION < v"1.1"
        y = setprecision(BigFloat, prec) do
            BigFloat()
        end
        ccall((:mpfr_set, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32),
                  y, x, Base.MPFR.ROUNDING_MODE[])
    else
        y = BigFloat(x, precision=prec)
    end

    u = zero(UInt128)
    i = cld(prec, sizeof(Culong)*8)
    j = 113
    while i > 0
        j -= sizeof(Culong)*8
        u |= (unsafe_load(y.d, i) % UInt128) << j
        i -= 1
    end
    u &= significand_mask(Float128)
    u |= exponent_half(Float128)
    z = ldexp(reinterpret(Float128, u), y.exp)
    return copysign(z,x)
end

function BigInt(x::Float128)
    !isinteger(x) && throw(InexactError(BigInt, x))
    BigInt(BigFloat(x, precision=precision(Float128)))
end
function Float128(x::BigInt)
    @static if VERSION < v"1.1"
        y = setprecision(BigFloat, precision(Float128)) do
            BigFloat(x)
        end
    else
        y = BigFloat(x, precision=precision(Float128))
    end
    Float128(y)
end

promote_rule(::Type{Float128}, ::Type{Float16}) = Float128
promote_rule(::Type{Float128}, ::Type{Float32}) = Float128
promote_rule(::Type{Float128}, ::Type{Float64}) = Float128
promote_rule(::Type{Float128}, ::Type{<:Integer}) = Float128

#widen(::Type{Float64}) = Float128
widen(::Type{Float128}) = BigFloat

# TODO: need to do this better
function parse(::Type{Float128}, s::AbstractString)
    Float128(ccall((:strtoflt128, libquadmath), Cfloat128, (Cstring, Ptr{Ptr{Cchar}}), s, C_NULL))
end

function string(x::Float128)
    lng = 64
    buf = Array{UInt8}(undef, lng + 1)
    lng = ccall((:quadmath_snprintf,libquadmath), Cint, (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Cfloat128...), buf, lng + 1, "%.35Qe", x)
    return String(resize!(buf, lng))
end

print(io::IO, b::Float128) = print(io, string(b))
show(io::IO, b::Float128) = print(io, string(b))

end # module Quadmath


text:
type MEngineError <: Exception
    message::ASCIIString
end


text:
# The quadratic segment: 𝗾(r) = r²𝘂 + r𝘃 + 𝘅₁
# 𝘂 = 2(𝘅₁ + 𝘅₂ - 2𝘅₃) and 𝘃 = -(3𝘅₁ + 𝘅₂ - 4𝘅₃)
# The line segment: 𝗹(s) = 𝘅₄ + s𝘄
# 𝘅₄ + s𝘄 = r²𝘂 + r𝘃 + 𝘅₁
# s𝘄 = r²𝘂 + r𝘃 + (𝘅₁ - 𝘅₄)
# 𝟬 = r²(𝘂 × 𝘄) + r(𝘃 × 𝘄) + (𝘅₁ - 𝘅₄) × 𝘄
# The cross product of two vectors in the plane is a vector of the form (0, 0, k).
# Let a = (𝘂 × 𝘄)ₖ, b = (𝘃 × 𝘄)ₖ, c = ([𝘅₁ - 𝘅₄] × 𝘄)ₖ
# 0 = ar² + br + c
# If a = 0
#   r = -c/b
# else
#   r = (-b ± √(b²-4ac))/2a
# We must also solve for s
# 𝘅₄ + s𝘄 = 𝗾(r)
# s𝘄 = 𝗾(r) - 𝘅₄
# s = ([𝗾(r) - 𝘅₄] ⋅𝘄 )/(𝘄 ⋅ 𝘄)
#
# r is invalid if:
#   1) b² < 4ac
#   2) r ∉ [0, 1]   (Curve intersects, segment doesn't)
# s is invalid if:
#   1) s ∉ [0, 1]   (Line intersects, segment doesn't)
function Base.intersect(l::LineSegment{Point{2,T}},
                        q::QuadraticSegment{Point{2,T}}) where {T}
    pmiss = Point(T(1e6),T(1e6))
    # Check if the segment is effectively straight.
    # Project P₃ onto the line from P₁ to P₂, call it P₄
    𝘃₁₃ = q[3] - q[1]
    𝘃₁₂ = q[2] - q[1]
    v₁₂ = norm²(𝘃₁₂)
    𝘃₁₄ = (𝘃₁₃ ⋅ 𝘃₁₂)*inv(v₁₂)*𝘃₁₂
    # Determine the distance from P₃ to P₄ (P₄ - P₃ = P₁ + 𝘃₁₄ - P₃ = 𝘃₁₄ - 𝘃₁₃)
    d² = norm²(𝘃₁₄ - 𝘃₁₃)
    if d² < ϵ_Point^2 # Use line segment intersection, segment is effectively straight
        # Line segment intersection looks like the following.
        # We want to reuse quantities we have already computed
        # Here l₁ = l, l₂ = LineSegment(q[1], q[2])
        #    𝘄 = l₂[1] - l₁[1]
        #    𝘂₁= l₁[2] - l₁[1]
        #    𝘂₂= l₂[2] - l₂[1]
        #    z = 𝘂₁ × 𝘂₂
        #    r = (𝘄 × 𝘂₂)/z
        #    s = (𝘄 × 𝘂₁)/z
        #    valid = 0 ≤ r && r ≤ 1 && 0 ≤ s && s ≤ 1
        𝘄 = q[1] - l[1]
        𝘂₁= l[2] - l[1]
        # 𝘂₂= 𝘃₁₂
        z = 𝘂₁ × 𝘃₁₂
        r = (𝘄 × 𝘃₁₂)/z
        s = (𝘄 × 𝘂₁)/z
        valid = 0 ≤ r && r ≤ 1 && 0 ≤ s && s ≤ 1
        return valid ? Vec(l(r), pmiss) : Vec(pmiss, pmiss)
    else
        𝘂 = 2𝘃₁₂ - 4𝘃₁₃
        𝘃 = 4𝘃₁₃ - 𝘃₁₂
        𝘄 = l[2] - l[1]
        a = 𝘂 × 𝘄
        b = 𝘃 × 𝘄
        c = (q[1] - l[1]) × 𝘄
        w² = 𝘄  ⋅ 𝘄  # 0 ≤ w²
        if a == 0
            r = -c/b
            0 ≤ r ≤ 1 || return Vec(pmiss, pmiss)
            p = q(r)
            s = (p - l[1]) ⋅ 𝘄
            # Since 0 ≤ w², we may test 0 ≤ s ≤ w², and avoid a division by
            # w² in computing s
            return 0 ≤ s && s ≤ w² ? Vec(p, pmiss) : Vec(pmiss, pmiss)
        elseif b^2 ≥ 4a*c
            r₁ = (-b - √(b^2 - 4a*c))/2a
            r₂ = (-b + √(b^2 - 4a*c))/2a
            p₁ = pmiss
            p₂ = pmiss
            if 0 ≤ r₁ ≤ 1
                x = q(r₁)
                s₁ = (x - l[1])⋅𝘄
                if 0 ≤ s₁ && s₁ ≤ w²
                    p₁ = x
                end
            end
            if 0 ≤ r₂ ≤ 1
                y = q(r₂)
                s₂ = (y - l[1])⋅𝘄
                if 0 ≤ s₂ && s₂ ≤ w²
                    p₂ = y
                end
            end
            return Vec(p₁, p₂)
        else
            return Vec(pmiss, pmiss)
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Debugger

println("hello")
function hhh(uio)
    println("hhh in ")

    println(uio + 5)
    retval = 4
    return retval
end

# this line will run debugger in console
# @enter hhh(4)

# this line  should be called from REPL
# Juno.@enter hhh(4)


text:
@testset "ModeEstimation.jl" begin
    @testset "gdemo" begin
        @testset "MLE" begin
            Random.seed!(222)
            true_value = [0.0625, 1.75]

            f1 = optim_function(gdemo_default, MLE();constrained=false)
            p1 = GalacticOptim.OptimizationProblem(f1.func, f1.init(true_value), nothing)

            p2 = optim_objective(gdemo_default, MLE();constrained=false)

            p3 = optim_problem(gdemo_default, MLE();constrained=false, init_theta=true_value)

            m1 = solve(p1, NelderMead())
            m2 = solve(p1, LBFGS())
            m3 = solve(p1, BFGS())
            m4 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), NelderMead())
            m5 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), LBFGS())
            m6 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), BFGS())
            m7 = solve(p3.prob, NelderMead())
            m8 = solve(p3.prob, LBFGS())
            m9 = solve(p3.prob, BFGS())

            @test all(isapprox.(f1.transform(m1.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m2.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m3.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m4.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m5.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m6.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m7.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m8.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m9.minimizer) - true_value, 0.0, atol=0.01))
        end

        @testset "MAP" begin
            Random.seed!(222)
            true_value = [49 / 54, 7 / 6]

            f1 = optim_function(gdemo_default, MAP();constrained=false)
            p1 = GalacticOptim.OptimizationProblem(f1.func, f1.init(true_value), nothing)

            p2 = optim_objective(gdemo_default, MAP();constrained=false)

            p3 = optim_problem(gdemo_default, MAP();constrained=false,init_theta=true_value)

            m1 = solve(p1, NelderMead())
            m2 = solve(p1, LBFGS())
            m3 = solve(p1, BFGS())
            m4 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), NelderMead())
            m5 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), LBFGS())
            m6 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), p2.init(true_value), BFGS())
            m7 = solve(p3.prob, NelderMead())
            m8 = solve(p3.prob, LBFGS())
            m9 = solve(p3.prob, BFGS())

            @test all(isapprox.(f1.transform(m1.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m2.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m3.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m4.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m5.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m6.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m7.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m8.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m9.minimizer) - true_value, 0.0, atol=0.01))
        end

        @testset "MLE constrained" begin
            Random.seed!(222)
            true_value = [0.0625, 1.75]
            lb = [0.0, 0.0]
            ub = [2.0, 2.0]

            f1 = optim_function(gdemo_default, MLE();constrained=true)
            p1 = GalacticOptim.OptimizationProblem(f1.func, f1.init(true_value), nothing; lb=lb, ub=ub)

            p2 = optim_objective(gdemo_default, MLE();constrained=true)

            p3 = optim_problem(gdemo_default, MLE();constrained=true, init_theta=true_value, lb=lb, ub=ub)

            m1 = solve(p1, Fminbox(LBFGS()))
            m2 = solve(p1, Fminbox(BFGS()))
            m3 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), lb, ub, p2.init(true_value), Fminbox(LBFGS()))
            m4 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), lb, ub, p2.init(true_value), Fminbox(BFGS()))
            m5 = solve(p3.prob, Fminbox(LBFGS()))
            m6 = solve(p3.prob, Fminbox(BFGS()))

            @test all(isapprox.(f1.transform(m1.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m2.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m3.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m4.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m5.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m6.minimizer) - true_value, 0.0, atol=0.01))
        end

        @testset "MAP constrained" begin
            Random.seed!(222)
            true_value = [49 / 54, 7 / 6]
            lb = [0.0, 0.0]
            ub = [2.0, 2.0]

            f1 = optim_function(gdemo_default, MAP();constrained=true)
            p1 = GalacticOptim.OptimizationProblem(f1.func, f1.init(true_value), nothing; lb=lb, ub=ub)

            p2 = optim_objective(gdemo_default, MAP();constrained=true)

            p3 = optim_problem(gdemo_default, MAP();constrained=true, init_theta=true_value, lb=lb, ub=ub)

            m1 = solve(p1, Fminbox(LBFGS()))
            m2 = solve(p1, Fminbox(BFGS()))
            m3 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), lb, ub, p2.init(true_value), Fminbox(LBFGS()))
            m4 = optimize(p2.obj, (G,z) -> p2.obj(nothing,G,nothing,z), lb, ub, p2.init(true_value), Fminbox(BFGS()))
            m5 = solve(p3.prob, Fminbox(LBFGS()))
            m6 = solve(p3.prob, Fminbox(BFGS()))

            @test all(isapprox.(f1.transform(m1.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(f1.transform(m2.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m3.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p2.transform(m4.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m5.minimizer) - true_value, 0.0, atol=0.01))
            @test all(isapprox.(p3.transform(m6.minimizer) - true_value, 0.0, atol=0.01))
        end
    end
end


text:
function FieldKeplerianOrbit(arg0::FieldOrbit)
    return FieldKeplerianOrbit((FieldOrbit,), arg0)
end

function FieldKeplerianOrbit(arg0::FieldPVCoordinates, arg1::Frame, arg2::FieldAbsoluteDate, arg3::RealFieldElement)
    return FieldKeplerianOrbit((FieldPVCoordinates, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3)
end

function FieldKeplerianOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg10::RealFieldElement, arg11::RealFieldElement, arg12::PositionAngle, arg13::Frame, arg14::FieldAbsoluteDate, arg15::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::RealFieldElement, arg7::RealFieldElement, arg8::RealFieldElement, arg9::RealFieldElement)
    return FieldKeplerianOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldKeplerianOrbit(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::RealFieldElement, arg6::PositionAngle, arg7::Frame, arg8::FieldAbsoluteDate, arg9::RealFieldElement)
    return FieldKeplerianOrbit((RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, PositionAngle, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function FieldKeplerianOrbit(arg0::TimeStampedFieldPVCoordinates, arg1::Frame, arg2::RealFieldElement)
    return FieldKeplerianOrbit((TimeStampedFieldPVCoordinates, Frame, RealFieldElement), arg0, arg1, arg2)
end

function add_kepler_contribution(obj::FieldKeplerianOrbit, arg0::PositionAngle, arg1::RealFieldElement, arg2::Vector{RealFieldElement})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, RealFieldElement, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function elliptic_eccentric_to_mean(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "ellipticEccentricToMean", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function elliptic_eccentric_to_true(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "ellipticEccentricToTrue", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function get_a(obj::FieldKeplerianOrbit)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_a_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getADot", RealFieldElement, ())
end

function get_anomaly(obj::FieldKeplerianOrbit, arg0::PositionAngle)
    return jcall(obj, "getAnomaly", RealFieldElement, (PositionAngle,), arg0)
end

function get_anomaly_dot(obj::FieldKeplerianOrbit, arg0::PositionAngle)
    return jcall(obj, "getAnomalyDot", RealFieldElement, (PositionAngle,), arg0)
end

function get_e(obj::FieldKeplerianOrbit)
    return jcall(obj, "getE", RealFieldElement, ())
end

function get_e_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEDot", RealFieldElement, ())
end

function get_eccentric_anomaly(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEccentricAnomaly", RealFieldElement, ())
end

function get_eccentric_anomaly_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEccentricAnomalyDot", RealFieldElement, ())
end

function get_equinoctial_ex(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEquinoctialEx", RealFieldElement, ())
end

function get_equinoctial_ex_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEquinoctialExDot", RealFieldElement, ())
end

function get_equinoctial_ey(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEquinoctialEy", RealFieldElement, ())
end

function get_equinoctial_ey_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getEquinoctialEyDot", RealFieldElement, ())
end

function get_hx(obj::FieldKeplerianOrbit)
    return jcall(obj, "getHx", RealFieldElement, ())
end

function get_hx_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getHxDot", RealFieldElement, ())
end

function get_hy(obj::FieldKeplerianOrbit)
    return jcall(obj, "getHy", RealFieldElement, ())
end

function get_hy_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getHyDot", RealFieldElement, ())
end

function get_i(obj::FieldKeplerianOrbit)
    return jcall(obj, "getI", RealFieldElement, ())
end

function get_i_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getIDot", RealFieldElement, ())
end

function get_le(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLE", RealFieldElement, ())
end

function get_le_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLEDot", RealFieldElement, ())
end

function get_lm(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lm_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLMDot", RealFieldElement, ())
end

function get_lv(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_lv_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getLvDot", RealFieldElement, ())
end

function get_mean_anomaly(obj::FieldKeplerianOrbit)
    return jcall(obj, "getMeanAnomaly", RealFieldElement, ())
end

function get_mean_anomaly_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getMeanAnomalyDot", RealFieldElement, ())
end

function get_perigee_argument(obj::FieldKeplerianOrbit)
    return jcall(obj, "getPerigeeArgument", RealFieldElement, ())
end

function get_perigee_argument_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getPerigeeArgumentDot", RealFieldElement, ())
end

function get_right_ascension_of_ascending_node(obj::FieldKeplerianOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNode", RealFieldElement, ())
end

function get_right_ascension_of_ascending_node_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getRightAscensionOfAscendingNodeDot", RealFieldElement, ())
end

function get_true_anomaly(obj::FieldKeplerianOrbit)
    return jcall(obj, "getTrueAnomaly", RealFieldElement, ())
end

function get_true_anomaly_dot(obj::FieldKeplerianOrbit)
    return jcall(obj, "getTrueAnomalyDot", RealFieldElement, ())
end

function get_type(obj::FieldKeplerianOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::FieldKeplerianOrbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function hyperbolic_eccentric_to_mean(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "hyperbolicEccentricToMean", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function hyperbolic_eccentric_to_true(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "hyperbolicEccentricToTrue", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function interpolate(obj::FieldKeplerianOrbit, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldKeplerianOrbit, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function mean_to_elliptic_eccentric(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "meanToEllipticEccentric", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function mean_to_hyperbolic_eccentric(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "meanToHyperbolicEccentric", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function shifted_by(obj::FieldKeplerianOrbit, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldKeplerianOrbit, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldKeplerianOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldKeplerianOrbit, (jdouble,), arg0)
end

function to_orbit(obj::FieldKeplerianOrbit)
    return jcall(obj, "toOrbit", KeplerianOrbit, ())
end

function to_string(obj::FieldKeplerianOrbit)
    return jcall(obj, "toString", JString, ())
end

function true_to_elliptic_eccentric(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "trueToEllipticEccentric", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function true_to_hyperbolic_eccentric(::Type{FieldKeplerianOrbit}, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(FieldKeplerianOrbit, "trueToHyperbolicEccentric", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end


text:
using Exercise1
using Test

@testset "Exercise1.jl" begin
    # Write your tests here.
    n = 10
    A = randn(n,n)
    B = randn(n,n)
    C = A*B
    Ctest = similar(C)
    rowmatmul!(Ctest, A, B)
    @test Ctest ≈ C
    colmatmul!(Ctest, A, B)
    @test Ctest ≈ C
end

@testset "dot versions" begin
    n = 10
    A = randn(n,n)
    B = randn(n,n)
    C = A*B
    Ctest = similar(C)
    rowmatmuldot!(Ctest, A, B)
    @test Ctest ≈ C
end


text:
"""
    Poisson

* Parametrization
parametrized in \\lambda

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
Poisson

function score!(score_til::Matrix{T}, y::Int, ::Type{Poisson}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = y/param[t, 1] - 1
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{Poisson}, param::Matrix{T}, t::Int) where T
    aux.fisher[1, 1] = 1/param[1]
    return
end

function log_likelihood(::Type{Poisson}, y::Vector{Int}, param::Matrix{T}, n::Int) where T
    loglik = zero(T)
    for t in 1:n
        loglik += y[t] * log(param[t, 1]) - param[t, 1] - logfactorial(y[t])
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{Poisson}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{Poisson}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Poisson}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
    return
end

# utils
function update_dist(::Type{Poisson}, param::Matrix{T}, t::Int) where T
    return Poisson(param[t, 1])
end

function params_sdm(d::Poisson)
    return Distributions.params(d)
end

function num_params(::Type{Poisson})
    return 1
end


text:
module cordoba
# Write your package code here.
import JuMP
import JSON
import CSV
import Memento
import PowerModels
import PowerModelsACDC
import FlexPlan
const _PM = PowerModels
const _PMACDC = PowerModelsACDC
const _FP = FlexPlan
import InfrastructureModels
const _IM = InfrastructureModels
const _MOI = _IM._MOI # MathOptInterface

using JuMP: with_optimizer, optimizer_with_attributes
export with_optimizer, optimizer_with_attributes

include("prob/cordoba_storage_tnep.jl")
include("prob/common.jl")
include("core/objective.jl")
include("core/storage.jl")
include("io/profile_data.jl")
include("post_process/printing.jl")
include("clustering/bins.jl")
include("clustering/common.jl")
include("economics/main.jl")
end


text:
# predicates

isnan(a::Float120)        = isnan(a.hi)
isinf(a::Float120)        = isinf(a.hi)
isfinite(a::Float120)     = isfinite(a.hi)
issubnormal(a::Float120)  = issubnormal(a.hi)

iszero(a::Float120)       = (a.hi == zero(typeof(a.hi)))
ispos(a::Float120)        = (a.hi >= zero(typeof(a.hi)))
ispoz(a::Float120)        = !signbit(a.hi)
isneg(a::Float120)        = signbit(a.hi)


text:
module FermionicOperatorStringAnalyzer
    using Reexport
    @reexport using SymEngine

    include("FO.jl")
    include("FOS.jl")
    include("vanisher.jl")
    include("simplify.jl")
end


text:
struct BamRegion
	bamFile::BamFile
	refID::Int # 0-based index, since that is what each read has stored
	region::UnitRange{Int}
end
#BamRegion{T<:Integer}(f::BamFile, refID::Int, region::UnitRange{T})

BamRegion(f::BamFile, refID::Int) = BamRegion(f,refID,1:sequences(f)[refID+1][2])
BamRegion(f::BamFile, refName::AbstractString) = BamRegion(f,refID(f,refName))
BamRegion(f::BamFile, refName::AbstractString, region::UnitRange) = BamRegion(f,refID(f,refName), region)

region(f::BamFile, args...) = BamRegion(f, args...)

function ChunkQueue(reg::BamRegion)
	@assert hasindex(reg.bamFile) "Missing bai index for $(filename(reg.bamFile))."
	ChunkQueue(index(reg.bamFile), reg.refID, reg.region)
end

# TODO: type or immutable?
mutable struct BamRegionState
	refID::Int
	region::UnitRange{Int}
	cq::ChunkQueue
	currChunk::Chunk
	beginOffs::VirtualOffset # linear index offset of the 16kbp window containing start of region
	s::BGZFBlockStream
	currRead::BamRead
	nextRead::BamRead
	isdone::Bool # useful since there are many indicators that we are done
end

function BamRegionState(reg::BamRegion)
	linearIndex = index(reg.bamFile).refs[reg.refID+1].linearIndex
	window = ((reg.region[1]-1)>>14)+1
	#beginOffs = linearIndex[window]
	beginOffs = window<=length(linearIndex) ? linearIndex[window] : VirtualOffset(0xffffffffffffffff)

	s = BamRegionState(reg.refID,
	                   reg.region,
	                   ChunkQueue(reg),
	                   Chunk(VirtualOffset(0),VirtualOffset(0)),
	                   beginOffs,
	                   BGZFBlockStream(filename(reg.bamFile)),
	                   BamRead(reg.bamFile),
	                   BamRead(reg.bamFile),
	                   false)
	# put in valid state such that we can call next() and done()

	nextimpl(reg.bamFile, s)
	s

	# s.done = isempty(s.cq)
	# s.done && return s # no chunks - done
	# nextimpl(s)
	# s
end

# start(reg::BamRegion) = BamRegionState(reg)

function nextimpl(bamFile::BamFile, s::BamRegionState)

	# println("--- nextimpl ---")

	while true

		# println("ftell: ", position(s.s))
		# println("currChunk end: ", s.currChunk.e)

		# move to a new chunk if needed
		if position(s.s) >= s.currChunk.e # after current chunk?
			while true
				isempty(s.cq) && (s.isdone=true; return)
				s.currChunk = shift!(s.cq) # get next chunk
				position(s.s) >= s.currChunk.e && continue # skip if it doesn't end after last chunk?

				# println("currChunk: ", s.currChunk)
				# println("linear index offset: ", s.beginOffs)

				#if true # isinteresting(currChunk)
				if s.currChunk.e > s.beginOffs # use linear index to check if chunk could be of interest
					# println("Using chunk: ", s.currChunk)

					# seek if next chunk starts after current
					position(s.s) < s.currChunk.s && seek(s.s, s.currChunk.s)
					break
				else
					# println("Skipping chunk: ", s.currChunk)
				end
			end
		end

		# println("block size: ", length(s.s.block))
		eof(s.s) && (isdone=true; return)

		# we are now in a chunk with interesting reads
		# println("nextread!")
		#nextread!(s.s,s.nextRead)
		trynextread!(s.s,s.nextRead) || (s.isdone=true; return)
		# println("after nextread!")

		# read starts beyond region?
		# println("pnext: ", pos(s.nextRead), ", region end: ", s.region[end])
		pos(s.nextRead) > s.region[end] && (s.isdone=true; return)

		# # read ends before region?
		# println("pendnext: ", endpos(s.nextRead), ", region start: ", s.region[1])
		# endpos(s.nextRead) < s.region[1] && continue

		# read ends before region?
		# endpos is much more expensive to compute, only do it for reads where it is needed!
		pos(s.nextRead) < s.region[1] && endpos(s.nextRead) < s.region[1] && continue

		break
	end
end

# function next(reg::BamRegion, s::BamRegionState)
# 	s.currRead,s.nextRead = s.nextRead,s.currRead # swap
# 	nextimpl(reg.bamFile, s)
# 	(s.currRead, s)
# end

# done(::BamRegion, s::BamRegionState) = s.isdone

function iterate(reg::BamRegion, s::BamRegionState=BamRegionState(reg))
	s.isdone && return nothing

	s.currRead,s.nextRead = s.nextRead,s.currRead # swap
	nextimpl(reg.bamFile, s)
	(s.currRead, s)
end

eltype(::Type{BamRegion}) = BamRead

# immutable BamState
# 	io::IO
# 	zio::IO
# 	read::BamRead
# end

# BamState(f::BamFile, io::IO, zio::IO) = BamState(io,zio,BamRead(f))
# function BamState(f::BamFile, io::IO)
# 	zio = ZlibInflateInputStream(io,reset_on_end=true)
# 	skip_bam_header(zio) # just skip header for now
# 	BamState(f, io, zio)
# end
# BamState(f::BamFile) = BamState(f,open(f.filename))

# Base.start(f::BamFile) = BamState(f)
# Base.next(::BamFile, s::BamState) = (nextread!(s.zio, s.read), s)
# Base.done(::BamFile, s::BamState) = eof(s.zio)
# Base.eltype(::Type{BamFile}) = BamRead