text:
module EmpiricalDist
include("utils.jl")
import .Utils: checkpkgs
checkpkgs("PrettyTables", "Statistics", "StatsBase", "Distributions", "RecipesBase", "DataFrames", "Interpolations", "KernelDensity")

import PrettyTables, Statistics, StatsBase, RecipesBase, DataFrames
import Random, Interpolations, KernelDensity, Plots, Distributions

abstract type AbstractSeries{T <: Real, S} <: AbstractVector{T} end
Base.eltype(::AbstractSeries{T}) where T = T
Base.keytype(::AbstractSeries{T,S}) where {T,S} = S
Base.size(x::AbstractSeries) = size(values(x))
Base.findfirst(p::Function, x::AbstractSeries) = findfirst(p, keys(x))
Base.argmax(x::AbstractSeries) = keyat(x, argmax(values(x)))
valarraytype(x::AbstractSeries) = typeof(values(x))
keyat(x::AbstractSeries, i) = keys(x)[i]
valueat(x::AbstractSeries, i) = values(x)[i]
valueat!(x::AbstractSeries, i, v) = values(x)[i] = v
Base.iterate(x::AbstractSeries, i=1) = (i - 1) < length(x) ? (valueat(x, i), i + 1) : nothing
Base.pairs(x::AbstractSeries) = map(x -> Pair(x[1], x[2]),zip(keys(x), values(x)))

_check_matching_lengths(v, k) = length(v) == length(k) || error("keys and values must have same length")

function makeseries(AS::Type{<:AbstractSeries}, p, q; normalize = false, sort = false, name = nothing)
    _check_matching_lengths(p, q)
    k, v = q, p
    if normalize
        v = p ./ sum(p)
    end
    if sort
        perm = sortperm(k)
        k, v = k[perm], v[perm]
    end
    AS(v, k, name)
end

makeseries(AS::Type{<:AbstractSeries}, p::Real, q; kwargs...) = makeseries(AS, fill(p, length(q)), q; kwargs...)
makeseries(AS::Type{<:AbstractSeries{T,S}}; name = nothing) where {T,S} = AS(Vector{T}(), Vector{S}(), name)

function makeseries(AS::Type{<:AbstractSeries}, p::AbstractVector, q::AbstractRange; kwargs...)
    :sort in keys(kwargs) && kwargs[:sort] && error("sort not supported for type $(typeof(q))")
    makeseries(AS, p, collect(q); kwargs...)
end

function makeseries(AS::Type{<:AbstractSeries}, x::AbstractDict; kwargs...)
    p = collect(x)
    k = map(x->x.first, p)
    v = map(x->x.second, p)
    makeseries(AS, v, k; kwargs...)
end

function makeseries(AS::Type{<:AbstractSeries}, a::AbstractVector; kwargs...)
    # if index appears in kwargs, `a` is *values*, otherwise, it's keys
    if :index in keys(kwargs)
        k = kwargs[:index]
        v = a
        kwargs = filter(x -> x.first != :index, kwargs) # remove index
    else
        c = collect(StatsBase.countmap(a))
        k = map(x->x.first, c)
        v = map(x->x.second, c)
        v = float(v) # force float
    end
    makeseries(AS, v, k; kwargs...)
end

function Base.first(x::T, n::Integer=1) where{T<:AbstractSeries}
    upperbound = n < length(x) ? n : :end
    makeseries(T, valueat(x, 1:upperbound), keyat(x, 1:upperbound); sort=false)
end

function Base.push!(x::AbstractSeries{T,S}, p::Pair{S, T}) where {T,S}
    push!(keys(x), p.first)
    push!(values(x), p.second)
    return x
end

Base.sort(x::AS) where {AS<:AbstractSeries} = makeseries(AS, valus(x), keys(x), sort=true)

function _check_matching_keys(series)::Bool
    if length(series) > 1
        k = keys(series[1])
        return all(s -> (keys(s) == k), series[2:end])
    end
    true
end

"""
creates an uninitialzied series of the same type as input series
if `empty` is not set, keys are copied to the new series
Otherwise, an empty series is returned

`name` is passed to the constructor as the new series name
"""
function Base.similar(x::T; empty::Bool = false, name=nothing) where{T<: AbstractSeries}
    if !empty
        v = valarraytype(x)(undef, size(x))
        k = copy(keys(x))
        T(v, k, name)
    else
        makeseries(T; name=name)
    end
end

findloc(x::AbstractSeries, i) = findfirst(==(i), keys(x)) # XXX pretty bad

function Base.push!(x::AbstractSeries, p::Pair{S, T}) where {T,S}
    push!(keys(x), p.first)
    push!(values(x), p.second)
    return x
end

function Base.getindex(x::AbstractSeries, i)
    if (loc = findloc(x, i)) !== nothing
        return valueat(x,loc)
    else
        throw(KeyError(i))
    end
end

function Base.getindex(x::AbstractSeries, i::BitVector)
    length(x) == length(i) || error("mask must have same length as series")
    makeseries(typeof(x), valueat(x, i), keyat(x, i))
end

(x::AbstractSeries{T})(i) where T = (loc = findloc(x, i); isnothing(loc) ? T(0) : valueat(x,loc))
(x::AbstractSeries{T})(a::AbstractVector) where T = [x(i) for i in a]

function Base.setindex!(x::AbstractSeries, v, i)
    if (loc = findloc(x, i)) !== nothing
        valueat!(x, loc, v)
    else
        push!(x, (i => v))
        return v
    end
end

colheader(x::AbstractSeries) = isnothing(name(x)) ? "values" : name(x)
function _show(io::IO, x::AbstractSeries; backend = Val(:text))
    PrettyTables.pretty_table(io, values(x); row_names=keys(x), header=[colheader(x)],
                    crop = :both, compact_printing = true, backend = backend)
end
Base.show(io::IO, m::MIME{Symbol("text/plain")}, x::AbstractSeries) = _show(io, x, backend = Val(:text))
Base.show(io::IO, m::MIME{Symbol("text/html")}, x::AbstractSeries) = _show(io, x, backend = Val(:html))
Base.show(io::IO, x::AbstractSeries) = _show(io, x)

struct SeriesStyle <: Base.Broadcast.BroadcastStyle end
Base.Broadcast.BroadcastStyle(::Type{<:AbstractSeries}) = SeriesStyle()
Base.Broadcast.BroadcastStyle(::SeriesStyle, ::SeriesStyle) = SeriesStyle()
Base.Broadcast.BroadcastStyle(::SeriesStyle, ::Base.Broadcast.BroadcastStyle) = SeriesStyle()
Base.Broadcast.BroadcastStyle(::Base.Broadcast.BroadcastStyle, ::SeriesStyle) = SeriesStyle()

function mergeseries(::typeof(+), x::T, y::T) where {T<:AbstractSeries}
    # XXX merged series is sorted. Original order is lost
    # XXX only defined for addition
    xdict = Dict(zip(keys(x), values(x)))
    ydict = Dict(zip(keys(y), values(y)))

    merged = merge(+, xdict, ydict)
    return makeseries(T, merged; sort=true)
end

for op in (*, +, -, /, ^)
    @eval function Base.broadcasted(::typeof($op), x::T, y) where {T<:AbstractSeries}
        return T($op.(values(x), y), keys(x), name(x))
    end
    # XXX doesn't support adding AbstractSeries{Int} + AbstractSeries{Float}
    @eval function Base.broadcasted(::typeof($op), x::T, y::T) where {T<:AbstractSeries}
        if _check_matching_keys((x, y))
            T($op.(values(x), values(y)), keys(x), name(x))
        else
            mergeseries($op, x, y)
        end
    end
    @eval function Base.broadcasted(::typeof($op), x, y::T) where {T<:AbstractSeries}
        return T($op.(x, values(y)), keys(y), name(y))
    end
end

function Base.broadcasted(::typeof(Base.literal_pow), ^, x::T, ::Val{y}) where {T<:AbstractSeries, y}
    return T(values(x) .^ y, keys(x), name(x))
end

for op in (>, <, ==, >=, <=)
    # boolean operators: return a bit vector
    @eval Base.broadcasted(::typeof($op), x::AbstractSeries, y) = Base.broadcasted($op, values(x), y)
    @eval Base.broadcasted(::typeof($op), x, y::AbstractSeries) = Base.broadcasted($op, x, values(y))
    @eval function Base.broadcasted(::typeof($op), x::AbstractSeries, y::AbstractSeries)
        if !_check_matching_keys((x, y))
            if ==($op,==) # yep, this is valid Julia code
                return false
            else
                error("$op is undefined for series with different keys")
            end
        end
        Base.broadcasted($op, values(x), values(y))
    end
end

Base.isapprox(p1::AbstractSeries, p2::AbstractSeries) = _check_matching_keys((p1, p2)) && all(Base.isapprox.(values(p1), values(p2)))

RecipesBase.@recipe function f(p::AbstractSeries; fill_below = nothing, fill_above=nothing)
    n = name(p)
    if !isnothing(n) && length(n) > 0
        label --> n
    end

    if get(plotattributes, :label, nothing) !== nothing
        legend --> true
    end

    if !isnothing(fill_below)
        RecipesBase.@series begin
            fill --> (0,0.5,:grey)
            primary := false
            p[keys(p) .<= fill_below]
        end
    end
    if !isnothing(fill_above)
        RecipesBase.@series begin
            fill --> (0,0.5,:grey)
            primary := false
            p[keys(p) .>= fill_above]
        end
    end

    keys(p), values(p)
end

Base.copy(src::T) where {T<:AbstractSeries} = makeseries(T, copy(values(src)), copy(keys(src)), name=name(src))

Base.copyto!(dest::AbstractVector, src::AbstractSeries) = copyto!(dest, values(src))
function Base.copyto!(d::T, bc::Base.Broadcast.Broadcasted) where {T<:AbstractSeries}
    bcf = Base.Broadcast.flatten(bc)
    data = Base.copyto!(values(d), bcf)
    T(data, keys(d), name(d))
end

function normalize!(x::AbstractSeries)
    s = sum(x)
    if !iszero(s) # XXX what's the correct behavior?
        copyto!(values(x), values(x) ./ s)
    end
    return s
end

DataFrames.transform(f::Function, x::AS) where{AS<:AbstractSeries} = makeseries(AS, copy(values(x)), f.(keys(x)))

"""
Add a series to a DataFrame as a new column.
The series name is used as the new column name if `colname` is not given
if `overwrite` is `true`, the new column will replace an existing column with same name, otherwise
(the default), it will raise an exception

The series keys must match the `index` column (if `index` doesn't exist, it will be added)
"""
function Base.insert!(df::DataFrames.AbstractDataFrame, x::AbstractSeries; colname = nothing, overwrite = false)
    n = isnothing(colname) ? name(x) : Symbol(colname)
    if !isnothing(n)
        if (string(n) in names(df)) && !overwrite
            # XXX more meaningful exception?
            error("column $n already exists in dataframe")
        end
    end

    k = keys(x)
    v = values(x)

    # check index column
    if "index" in names(df)
        if k != df[!, :index]
            error("series keys do not match index column")
        end
    else
        df[!, "index"] = k
    end

    if !isnothing(n)
        df[!, n] = v
    else
        # use default column name
        df[!, end+1] = v
    end
    nothing
end

function make_dataframe(x::T) where {T<:AbstractVector{<:AbstractSeries}}
    df = DataFrames.DataFrame()
    for s in x
        append!(df, Dict(zip(Symbol.(keys(s)), values(s))), cols = :union)
    end

    return df
end

make_dataframe(x::Tuple{Vararg{AbstractSeries}}) = make_dataframe(collect(x))

struct Series{T,S} <: AbstractSeries{T,S}
    values::Vector{T}
    index::Vector{S}
    name::Union{AbstractString, Nothing}
end
Base.values(x::Series) = x.values
Base.keys(x::Series) = x.index
name(x::Series) = x.name

struct Pmf{T <: AbstractFloat, S} <: AbstractSeries{T,S}
    ps::Vector{T}
    qs::Vector{S}
    name::Union{AbstractString, Nothing}
end

Pmf(a::AbstractSeries) = Pmf(values(a), keys(a), name(a))

# extract/deduce T, S
const _deffloattype = typeof(float(0))
Pmf(; kwargs...) = makeseries(Pmf{_deffloattype, Any}; kwargs...) # default floating point
Pmf{T,S}(; kwargs...) where {T<:Real, S} = makeseries(Pmf{T, S}; kwargs...)

Pmf(x::AbstractDict{K,V}; sort=true, kwargs...) where {K,V}= makeseries(Pmf{V,K}, x; sort=sort, kwargs...)
Pmf(a::AbstractVector; kwargs...) where S = makeseries(Pmf, a; kwargs...)
Pmf(p, q; kwargs...) where T = makeseries(Pmf, float.(p), q; kwargs...) # force floating point
Pmf(n::Real, s::AbstractSet; kwargs...) = makeseries(Pmf, n, collect(s); kwargs...)

function makeuniform(q; kwargs...)
    pmf = Pmf(1.0, q; kwargs...)
    normalize!(pmf)
    return pmf
end

colheader(x::Pmf) = isnothing(name(x)) ? "probs" : name(x)

Base.values(x::Pmf) = x.ps
Base.keys(x::Pmf) = x.qs
name(x::Pmf) = x.name

function Base.sort!(x::Pmf)
    p = sortperm(x.qs)
    x.qs .= x.qs[p]
    x.ps .= x.ps[p]
    x
end

"""Probability of quantities greater than threshold."""
prob_ge(p::Pmf{<: Real}, threshold) = sum(p[keys(p) .>= threshold])
prob_gt(p::Pmf{<: Real}, threshold) = sum(p[keys(p) .> threshold])

"""Probability of quantities less than threshold."""
prob_le(p::Pmf{<: Real}, threshold) = sum(p[keys(p) .<= threshold])
prob_lt(p::Pmf{<: Real}, threshold) = sum(p[keys(p) .< threshold])

prob_eq(p::Pmf{<: Real}, threshold) = sum(p[keys(p) .== threshold])

function pmf_outer(op::Function, p1::Pmf{<:Real, <:Real}, p2::Pmf{<:Real, <:Real})
    qm = op.(keys(p1), keys(p2)');
    pm = values(p1) .* values(p2)';
    return sum(qm .* pm)
end

prob_ge(p1::Pmf{<: Real, <: Real}, p2::Pmf{<: Real, <: Real}) = pmf_outer(>=, p1, p2)
prob_gt(p1::Pmf{<: Real, <: Real}, p2::Pmf{<: Real, <: Real}) = pmf_outer(>, p1, p2)
prob_le(p1::Pmf{<: Real, <: Real}, p2::Pmf{<: Real, <: Real}) = pmf_outer(<=, p1, p2)
prob_lt(p1::Pmf{<: Real, <: Real}, p2::Pmf{<: Real, <: Real}) = pmf_outer(<, p1, p2)
prob_eq(p1::Pmf{<: Real, <: Real}, p2::Pmf{<: Real, <: Real}) = pmf_outer(==, p1, p2)

ge_dist = prob_ge
gt_dist = prob_gt
le_dist = prob_le
lt_dist = prob_lt
eq_dist = prob_eq

maxprob(x::Pmf) = argmax(x)
pmffromseq(seq::AbstractVector; sort=true, normalize=true, name=nothing) = Pmf(seq; sort, normalize, name)
pmffromseq(seq::AbstractArray; sort=true, normalize=true, name=nothing)  = pmffromseq(vec(seq); sort, normalize, name)

"""Make a discrete approximation.

dist: distribution object
qs: quantities

returns: Pmf
"""
function pmffromdist(dist, qs)
    ps = Distributions.pdf.(dist, qs)
    pmf = Pmf(ps, qs)
    normalize!(pmf)
    return pmf
end

Statistics.mean(x::Pmf) = sum(values(x) .* keys(x))
Statistics.var(x::Pmf) = sum((keys(x) .- Statistics.mean(x)) .^ 2 .* values(x))
Statistics.std(x::Pmf) = sqrt(Statistics.var(x))

function Statistics.quantile(x::Pmf, a::AbstractVector)
    # XXX we could just rely on Cdf.
    # E.g.:
    #    quantile(makecdf(x), a)
    k = keys(x)
    v = cumsum(values(x))
    [k[searchsortedfirst(v, i)] for i in a]
end

Statistics.quantile(x::Pmf, i::T) where{T<:AbstractFloat} = Statistics.quantile(x,[i])[1]
Statistics.median(x::Pmf)  = Statistics.quantile(x,[0.5])[1]

StatsBase.sample(p::Pmf) = StatsBase.sample(p.qs, StatsBase.pweights(p.ps))
StatsBase.sample(p::Pmf, n::Integer) = StatsBase.sample(p.qs, StatsBase.pweights(p.ps), n)
Random.rand(p::Pmf) = StatsBase.sample(p)
Random.rand(p::Pmf, n::Integer) = StatsBase.sample(p, n)

credibleinterval(x::Pmf, p::Real) = Statistics.quantile(x, [(1-p)/2, 1-(1-p)/2])

function convolvedist(ufunc::Function, xs::Vararg{Pmf{T,S}}) where {T<:Real, S<:Real}
    length(xs) > 1 || error("at least 2 series are needed for convolvedist")
    qs = keys(xs[1])
    ps = values(xs[1])
    for i in 2:length(xs)
        q, p = keys(xs[i]), values(xs[i])
        qs = vec(ufunc.(qs,q'))
        ps = vec(ps .* p')
    end
    d = StatsBase.addcounts!(Dict{S, T}(), qs, ps)
    pmf = Pmf(d, sort=true)
end

adddist(xs::Vararg{Pmf}) = convolvedist(+, xs...)
adddist(xs::AbstractVector{T}) where {T<:Pmf} = convolvedist(+, xs...)
adddist(x::Pmf, s::T) where{T<:Real} = Pmf(values(x), keys(x) .+ s)

subdist(xs::Vararg{Pmf}) = convolvedist(-, xs...)
subdist(xs::AbstractVector{T}) where {T<:Pmf} = convolvedist(-, xs...)
subdist(x::Pmf, s::T) where{T<:Real} = Pmf(values(x), keys(x) .- s)

muldist(xs::Vararg{Pmf}) = convolvedist(*, xs...)
muldist(xs::AbstractVector{T}) where {T<:Pmf} = convolvedist(*, xs...)

divdist(xs::Vararg{Pmf}) = convolvedist(/, xs...)
divdist(xs::AbstractVector{T}) where {T<:Pmf} = convolvedist(/, xs...)

function expand_seq(sp::AbstractVector{<:Pmf{T}}) where {T<:Real}
    # quite inefficient; assumes keys are sorted
    ks = sort(union([keys(p) for p in sp]...))
    k2l = Dict([x[2] =>x[1] for x in enumerate(ks)])
    m = zeros(T, (length(sp), length(ks)))
    for (i,s) in enumerate(sp)
        for (q,p) in pairs(s)
            m[i, k2l[q]] = p
        end
    end
    return (m, ks)
end

function makemixture(p::Pmf, s::AbstractVector{<:Pmf}; name = nothing)
    # XXX assumes Pmfs have same "base" (or can the extended to one)
    m, k = expand_seq(s)
    v = vec(sum(values(p) .* m, dims=1))
    Pmf(v, k, name)
end

"""Make a kernel density estimate from a sample."""
function kde_from_sample(sample::AbstractVector{<: Real}, qs::AbstractVector{<:Real})
    ps = Distributions.pdf(KernelDensity.InterpKDE(KernelDensity.kde(sample)), qs)
    pmf = Pmf(ps, qs)
    normalize!(pmf)
    return pmf
end

KernelDensity.kde(p::Pmf) = KernelDensity.kde(p.qs; weights = StatsBase.pweights(p.ps))

function marginal(p::Pmf{<:Real, <:Tuple}, d::Integer)
    _extract_tuple_ntype(::Type{<:Tuple{T}}, ::Val{1}) where T = T
    _extract_tuple_ntype(::Type{<:Tuple{T, <:Any}}, ::Val{1}) where T = T
    _extract_tuple_ntype(::Type{<:Tuple{<:Any, T}}, ::Val{2}) where T = T
    _extract_tuple_ntype(::Type{<:Tuple{T, <:Any, <:Any}}, ::Val{1}) where T = T
    _extract_tuple_ntype(::Type{<:Tuple{<:Any, T, <:Any}}, ::Val{2}) where T = T
    _extract_tuple_ntype(::Type{<:Tuple{<:Any, <:Any, T}}, ::Val{3}) where T = T

    T = eltype(keys(p))
    V = eltype(values(p))
    K = _extract_tuple_ntype(T, Val{d}())
    dict = Dict{K, V}()
    for i in 1:length(p)
        k = keyat(p, i)[d]
        (k in keys(dict)) || (dict[k] = 0)
        dict[k] += valueat(p, i)
    end
    Pmf(dict)
end

struct Cdf{T <: AbstractFloat, S} <: AbstractSeries{T,S}
    cps::Vector{T}
    qs::Vector{S}
    name::Union{AbstractString, Nothing}
end

Cdf(; kwargs...) = makeseries(Cdf{_deffloattype, Any}; kwargs...)
Cdf{T,S}(; kwargs...) where {T<:Real, S} = makeseries(Cdf{T, S}; kwargs...)
Cdf(p, q; kwargs...) where T = makeseries(Cdf, float.(p), q; kwargs...) # force floating point

Base.values(x::Cdf) = x.cps
Base.keys(x::Cdf) = x.qs
name(x::Cdf) = x.name

Base.cumsum(x::AbstractSeries) = Cdf(cumsum(values(x)), collect(keys(x)), nothing)

function makecdf(x::Pmf; normalize = true)
    cdf = cumsum(x)
    if normalize
        normalize!(cdf)
    end
    return cdf
end

function cdffromseq(x::AbstractVector; normalize = true, name = nothing)
    pmf = makeseries(Pmf, x; normalize=false, sort = true, name = name)
    makecdf(pmf; normalize)
end

# TODO: memoize/cache the result (or store in Cdf struct)
interpolation_func(s, d) = Interpolations.LinearInterpolation(s, d)
forward_inerpolate(x::Cdf, s) = interpolation_func(keys(x), values(x))(s)
(x::Cdf)(i) = forward_inerpolate(x, i)

# see https://discourse.julialang.org/t/problem-with-interpolations-did-it-change/60685
_prepare_for_intepolation(x, y) = (idx = unique(i -> x[i], 1:length(x)); (x[idx], y[idx]))
backward_interpolate(x::Cdf, s) = interpolation_func(_prepare_for_intepolation(values(x), keys(x))...)(s)
Statistics.quantile(x::Cdf, s) = backward_interpolate(x, s)
credibleinterval(x::Cdf, p::Real) = Statistics.quantile(x, [(1-p)/2, 1-(1-p)/2])

colheader(x::Cdf) = isnothing(name(x)) ? "cumulative" : name(x)
normalize!(x::Cdf) = copyto!(values(x), values(x) ./ values(x)[end])

function makepmf(x::Cdf{T, S}) where {T<: AbstractFloat, S}
    v = diff(vcat(zeros(T, 1), values(x))) # need to "restore" first element
    return Pmf(v, keys(x), name(x))
end

Statistics.mean(x::Cdf) = Statistics.mean(makepmf(x))
Statistics.var(x::Cdf) = Statistics.var(makepmf(x))
Statistics.std(x::Cdf) = Statistics.std(makepmf(x))

prob_ge(p::Cdf, threshold) = prob_ge(makepmf(p), threshold)
prob_gt(p::Cdf, threshold) = prob_gt(makepmf(p), threshold)
prob_le(p::Cdf, threshold) = prob_le(makepmf(p), threshold)
prob_lt(p::Cdf, threshold) = prob_lt(makepmf(p), threshold)
prob_eq(p::Cdf, threshold) = prob_lt(makepmf(p), threshold)

"""Make a binomial Pmf."""
makebinomial(n, p; kwargs...) = makeseries(Pmf, Distributions.pdf.(Distributions.Binomial(n,p), 0:n), 0:n; kwargs...)

maxdist(x::Cdf, n::Real) = Cdf(values(x) .^ n, keys(x), nothing)
mindist(x::Cdf, n::Real) = Cdf(1 .- (1 .- values(x)) .^ n, keys(x), nothing)

abstract type AbstractJointDistribution end

Base.size(j::AbstractJointDistribution) = size(distribution(j))
Base.maximum(j::AbstractJointDistribution) = maximum(distribution(j))
Base.sum(j::AbstractJointDistribution) = sum(distribution(j))
Base.length(j::AbstractJointDistribution) = length(rows(j))
Base.first(j::JD, n::Integer=1) where{JD<:AbstractJointDistribution} = JD(distribution(j)[1:n, :], rows(j)[1:n], columns(j))
Base.transpose(j::JD)  where{JD<:AbstractJointDistribution} = JD(distribution(j)', columns(j), rows(j))
Base.adjoint(j::AbstractJointDistribution) = transpose(j)

function DataFrames.stack(j::AbstractJointDistribution)
    m, x, y = distribution(j), columns(j), rows(j)
    q = [(y[u]..., x[v]...) for u in 1:length(y) for v in 1:length(x)]
    p = [m[u, v] for u in 1:length(y) for v in 1:length(x)]
    return Pmf(p, q)
end

function marginal(j::AbstractJointDistribution, dim)
    if (dim == :x) || (dim == 1)
        return Pmf(vec(sum(distribution(j), dims=1)), columns(j))
    elseif (dim == :y) || (dim == 2)
        return Pmf(vec(sum(distribution(j), dims=2)), rows(j))
    else
        error("unknown dimension")
    end
end

Plots.contour(j::AbstractJointDistribution, args...; levels = 4, kwargs...) =
                Plots.contour(columns(j), rows(j), distribution(j), args...; levels=levels, kwargs... )
Plots.contour!(j::AbstractJointDistribution, args...; levels = 4, kwargs...) =
                Plots.contour!(columns(j), rows(j), distribution(j), args...; levels=levels, kwargs... )

function _show(io::IO, j::AbstractJointDistribution; backend)
    PrettyTables.pretty_table(io, distribution(j); row_names=rows(j), header=columns(j),
                    crop = :both, compact_printing = true, backend = backend)
end

Base.show(io::IO, m::MIME{Symbol("text/plain")}, x::AbstractJointDistribution) = _show(io, x, backend = Val(:text))
Base.show(io::IO, m::MIME{Symbol("text/html")}, x::AbstractJointDistribution) = _show(io, x, backend = Val(:html))

function Base.argmax(j::AbstractJointDistribution)
    m, x, y = distribution(j), columns(j), rows(j)
    p = argmax(m)
    return y[p[1]], x[p[2]]
end

struct JDStyle <: Base.Broadcast.BroadcastStyle end
Base.Broadcast.BroadcastStyle(::Type{AbstractJointDistribution}) = JDStyle()
Base.Broadcast.BroadcastStyle(::JDStyle, ::JDStyle) = JDStyle()
Base.Broadcast.BroadcastStyle(::JDStyle, ::Base.Broadcast.BroadcastStyle) = JDStyle()
Base.Broadcast.BroadcastStyle(::Base.Broadcast.BroadcastStyle, ::JDStyle) = JDStyle()

function Base.broadcasted(::typeof(*), j::J, k::J) where{J<:AbstractJointDistribution}
    m1, x1, y1 = distribution(j), columns(j), rows(j)
    m2, x2, y2 = distribution(k), columns(k), rows(k)

    (x1 == x2 && y1 == y2) || error("dimensions must be the same")
    return J(m1 .* m2, y1, x1)
end

function Base.broadcasted(::typeof(*), j::J, k) where{J<:AbstractJointDistribution}
    m, x, y = distribution(j), columns(j), rows(j)
    return J(m .* k, y, x)
end

Base.copy(j::J) where{J<:AbstractJointDistribution} = J(copy(distribution(j)), copy(rows(j)), copy(columns(j)))

function Base.similar(j::J) where{J<: AbstractJointDistribution}
    m = similar(distribution(j))
    J(m, rows(j), columns(j))
end

# TODO: support naming axes
struct JointDistribution{M<:AbstractMatrix{<:Real}, V<:AbstractVector, U<:AbstractVector} <: AbstractJointDistribution
    m::M
    y::V # ↓
    x::U # →
end

distribution(j::JointDistribution) = j.m
columns(j::JointDistribution) = j.x
rows(j::JointDistribution) = j.y
index(j::JointDistribution) = rows(j)

function DataFrames.unstack(p::Pmf{<:Real, <:Tuple{<:Any, <:Any}})
    x = unique([k[2] for k in keys(p)])
    y = unique([k[1] for k in keys(p)])
    lenx = length(x)
    leny = length(y)

    m = [valueat(p, (i-1)*lenx+j) for i in 1:leny, j in 1:lenx]
    return JointDistribution(m, y, x)
end

DataFrames.unstack(p::Pmf{<:Real, <:Tuple{<:Any}}) = Pmf(values(p), [k[1] for k in keys(p)], name(p))

function DataFrames.unstack(p::Pmf{<:Real, <:Tuple})
    x = unique([k[end] for k in keys(p)])
    y = unique([k[1:end-1] for k in keys(p)])
    lenx = length(x)
    leny = length(y)

    m = [valueat(p, (i-1)*lenx+j) for i in 1:leny, j in 1:lenx]
    return JointDistribution(m, y, x)
end

normalize!(j::JointDistribution) = (s = sum(j.m); j.m ./= sum(j.m); s)
makejoint(a::Pmf, b::Pmf) = JointDistribution(values(a)' .* values(b), keys(b), keys(a))

Statistics.cov(df::DataFrames.AbstractDataFrame) = Statistics.cov(Matrix(df))

StatsBase.mean_and_cov(df::DataFrames.AbstractDataFrame) = StatsBase.mean_and_cov(Matrix(df))

"""
calculates a covariance matrix from a dataframe and returns a dataframe with an index column
"""
covdf(df::DataFrames.AbstractDataFrame) =
        DataFrames.DataFrame("index" => names(df),
                        map(i -> (names(df)[i] => Statistics.cov(df)[:, i]), 1:DataFrames.ncol(df))...)

"""Print the mean and CI of a distribution.

posterior: Pmf
digits: number of digits to round to
prob: probability in the CI
"""
function summarize(posterior::Pmf, digits::Integer=3, prob=0.9)
    μ = round(Statistics.mean(posterior), digits=3)
    ci = credibleinterval(posterior, prob)
    println("$μ, $ci")
end

export Pmf, Cdf, JointDistribution
export makepmf, normalize!, maxprob, pmffromseq, prob_ge, prob_le, prob_gt, prob_lt, maxprob, credibleinterval
export adddist, subdist, muldist, divdist, makebinomial, makecdf, cdffromseq, make_dataframe, maxdist, mindist
export prob_gt, prob_eq, lt_dist, gt_dist, le_dist, ge_dist, makemixture, makejoint, makeuniform, distribution
export summarize, covdf, marginal, rows, columns, kde_from_sample, pmffromdist, transform

end; # module


text:
const _default_maxbuf_size = 2048

"""
Documentation
"""
:readarrays, :readarrays!

@generated function readarrays!(toks::SpaceTokenizer, as...)
  N = length(as)
  Ts = map(eltype, as)

  # the basic loop is simple, we just keep reading!
  expr = quote
    while true
    end
    return as
  end

  for i=1:length(as)
    expr_read = quote
        curlen = step!(toks)
        if curlen <= 0
          break
        end
        push!(as[$i], myparse($(Ts[i]), toks.buf, 1, curlen))
    end
    push!(expr.args[2].args[2].args, expr_read)
  end
  return expr
end

function readarrays!(io, as...; maxbuf=_default_maxbuf_size)
  toks = SpaceTokenizer(io, maxbuf)
  return readarrays!(toks, as...)
end

function readarrays!(filename::AbstractString, as...; kwargs...)
  open(filename, "r") do fh
    return readarrays!(fh, as...; kwargs...)
  end
end

readarrays(io, Ts...; kwargs...) =
  readarrays!(io, map(T -> zeros(T,0), Ts)...; kwargs...)

function readarrays(filename::AbstractString, as...; kwargs...)
  open(filename, "r") do fh
    return readarrays(fh, as...; kwargs...)
  end
end

"""
"""
:readarray, :readarray!

readarray!(io, a; kwargs...) = readarrays!(io, a; kwargs...)[1]

readarray(::Type{T}, io; kwargs...) where {T <: ParsableNumbers} =
  readarray!(io, zeros(T,0); kwargs...)
readarray(io; kwargs...) = readarray(Float64, io; kwargs...)

#=
# Given a filename, refer to the type with an IO
function readarray!(filename::AbstractString, a; kwargs...)
  open(filename, "r") do fh
    return readarray!(fh, a; kwargs...)
  end
end

readarray(io; kwargs...) = readarray(Float64, io; kwargs...)
readarray(::Type{T}, io; kwargs...) where {T <: ParsableNumbers} =
  readarray!(io, zeros(T,0); kwargs...)

# Everything with filenames refers to the types of io
readarray(filename::AbstractString; kwargs...) =
  readarray(Float64, filename; kwargs...)

function readarray(::Type{T}, filename::AbstractString; kwargs...) where {T <: ParsableNumbers}
  open(filename, "r") do fh
    return readarray(T, fh; kwargs...)
  end
end
=#

"""
Documentation
"""
:readmatrix

function readmatrix(::Type{T}, io; transpose::Bool=true, kwargs...) where {T <: ParsableNumbers}
  firstline = readline(io) # read the first line
  a = zeros(T, 0)
  readarray!(IOBuffer(firstline), a; kwargs...)
  ncols = length(a)
  readarray!(io, a; kwargs...)
  nrows, nrem = divrem(length(a),ncols)
  if nrem != 0
    throw(ArgumentError("first line had $(ncols) entries," *
      " but overall there were $(length(a)) entries for " *
      "$(nrows) rows and $(nrem) left"))
  else
    if transpose
      return Base.transpose(reshape(a, ncols, nrows))
    else
      return reshape(a, ncols, nrows)
    end
  end
end
readmatrix(io; kwargs...) = readmatrix(Float64, io; kwargs...)

function readmatrix(::Type{T}, filename::AbstractString; kwargs...) where {T <: ParsableNumbers}
  open(filename, "r") do fh
    return readmatrix(T, fh; kwargs...)
  end
end
readmatrix(filename::AbstractString; kwargs...) = readmatrix(Float64, filename; kwargs...)

ext: jl


text:
import RLCore
import Flux

mutable struct GVFNetwork{M, H<:RLCore.AbstractHorde, N}
    gvf_model::M
    horde::H
    dqn_model::N
end

(m::GVFNetwork)(x) = m.dqn_model(Flux.data(m.gvf_model(x)))


text:
function LText(parent::Scene, text; kwargs...)
    LText(parent; text = text, kwargs...)
end

function LText(parent::Scene; bbox = nothing, kwargs...)

    default_attrs = default_attributes(LText, parent).attributes
    theme_attrs = subtheme(parent, :LText)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (text, textsize, font, color, visible, halign, valign,
        rotation, padding)

    layoutobservables = LayoutObservables{LText}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox)

    textpos = Node(Point3f0(0, 0, 0))

    # this is just a hack until boundingboxes in abstractplotting are perfect
    alignnode = lift(halign, rotation) do h, rot
        # left align the text if it's not rotated and left aligned
        if rot == 0 && (h == :left || h == 0.0)
            (:left, :center)
        else
            (:center, :center)
        end
    end

    t = text!(parent, text, position = textpos, textsize = textsize, font = font, color = color,
        visible = visible, align = alignnode, rotation = rotation, raw = true)[end]

    textbb = Ref(BBox(0, 1, 0, 1))

    onany(text, textsize, font, rotation, padding) do text, textsize, font, rotation, padding
        textbb[] = FRect2D(boundingbox(t))
        autowidth = width(textbb[]) + padding[1] + padding[2]
        autoheight = height(textbb[]) + padding[3] + padding[4]
        layoutobservables.autosize[] = (autowidth, autoheight)
    end

    onany(layoutobservables.computedbbox, padding) do bbox, padding

        tw = width(textbb[])
        th = height(textbb[])

        box = bbox.origin[1]
        boy = bbox.origin[2]

        # this is also part of the hack to improve left alignment until
        # boundingboxes are perfect
        tx = if rotation[] == 0 && (halign[] == :left || halign[] == 0.0)
            box + padding[1]
        else
            box + padding[1] + 0.5 * tw
        end
        ty = boy + padding[3] + 0.5 * th

        textpos[] = Point3f0(tx, ty, 0)
    end

    # trigger first update, otherwise bounds are wrong somehow
    text[] = text[]
    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    lt = LText(parent, layoutobservables, t, attrs)

    lt
end

function Base.delete!(lt::LText)
    GridLayoutBase.disconnect_layoutobservables!(lt.layoutobservables.gridcontent)
    GridLayoutBase.remove_from_gridlayout!(lt.layoutobservables.gridcontent)
    empty!(lt.layoutobservables.suggestedbbox.listeners)
    empty!(lt.layoutobservables.computedbbox.listeners)
    empty!(lt.layoutobservables.reportedsize.listeners)
    empty!(lt.layoutobservables.autosize.listeners)
    empty!(lt.layoutobservables.protrusions.listeners)

    # remove the plot object from the scene
    delete!(lt.parent, lt.textobject)
end


text:
module JinGen

export JinGenParams, JinGenBuilder, FlawedJinGen, generate

using RouletteWheels
using Compat

immutable JinGenParams
    r_0::Float64
    r_1::Float64
    γ::Float64
    z_star::Int
    barrier_iters::Int
    n::Int
end

JinGenParams(n) = JinGenParams(
    0.0005, 2.0, 0.005, 5, 1, n
)

immutable JinGenBuilder
    params::JinGenParams

    edges::Set{Tuple{Int,Int}}
    adj_list::Dict{Int,Vector{Int}}
    degrees::RouletteWheel
    mutual_meet_factors::RouletteWheel
    num_pairs::Int
    n_meetings::Int

    function JinGenBuilder(params)
        freqs = zeros(Int, params.n)
        num_pairs = @compat(round(Int, params.n * (params.n - 1) / 2))
        new(
            params,
            @compat(Set{Tuple{Int,Int}}()),
            @compat((Int => Vector{Int})[i => Vector{Int}() for i in 1:params.n]),
            StochasticAcceptance(freqs),
            StochasticAcceptance(freqs),
            num_pairs,
            @compat(round(Int, num_pairs * params.r_0))
        )
    end
end

"""
    :return: a random number, i, in [0, n), except for i == excluding.
"""
function random_idx(n, excluding)
    selected = excluding
    r = 1:n
    while selected == excluding
        selected = rand(r)
    end
    return selected
end

"""
    :return: a tuple of (i, j) representing an edge in a network of [0, n)
             nodes. This tuple has the property of i < j, forall edges.
"""
function random_edge_as_indices(n)
    a = rand(1:n)
    b = random_idx(n, a)

    # The following constraint is useful. Consider a set of edges in a
    # undirected graph. In this case, (1, 2) is semantically equivalent to
    # (2, 1). But, from the perspective of a set, these tuples are unique.
    # Using the canonical form of smaller number first enforces the correct
    # semantics.
    return a < b ? (a, b) : (b, a)
end

"""
:return: A unique pair of randomly selected nodes in edge canonical form.
"""
function random_pair_from_list(nodes)
    i, j = random_edge_as_indices(length(nodes))
    a, b = nodes[i], nodes[j]

    return a < b ? (a, b) : (b, a)  # Canonical form.
end

############################################################################
# All add and remove operations on edges must use add_edge! and remove_edge!
# respectively. They maintain the set of edges, as expected. They
# also maintain the degrees and mutual meet factors (i.e. k * (k-1)) for
# each node.
#
# Both these functions require edges in canonical form. That is, the smaller
# index must come first in the pair of indices representing the edge.
############################################################################

function inc_edge!(jingen::JinGenBuilder, a, b)
    push!(jingen.adj_list[a], b)
    jingen.degrees[a] += 1
    k = jingen.degrees[a]
    jingen.mutual_meet_factors[a] = k * (k - 1)
end

function add_edge!(jingen, edge)
    assert(edge[1] < edge[2])

    push!(jingen.edges, edge)
    inc_edge!(jingen, edge[1], edge[2])
    inc_edge!(jingen, edge[2], edge[1])
end

function dec_edge!(jingen::JinGenBuilder, a, b)
    deleteat!(
        jingen.adj_list[a],
        findfirst(jingen.adj_list[a], b) # FIX: XXX: High z* problem
    )
    jingen.degrees[a] -= 1
    k = jingen.degrees[a]
    jingen.mutual_meet_factors[a] = k * (k - 1)
end

function remove_edge!(jingen, edge)
    assert(edge[1] < edge[2])

    delete!(jingen.edges, edge)
    dec_edge!(jingen, edge[1], edge[2])
    dec_edge!(jingen, edge[2], edge[1])
end

# From Jin et al (2001):
#   If a pair meet and there is not already a connection between them, a new
#   connection is established unless one of them already has z* connections,
#   in which case nothing happens. In other words, z* forms a hard upper
#   limit on the degree z of any vertex, beyond which no more edges can
#   be added.
function connection_can_be_made(jingen, edge)
    if edge in jingen.edges
        return false
    end

    a, b = edge
    return (
        jingen.degrees[a] < jingen.params.z_star &&
        jingen.degrees[b] < jingen.params.z_star
    )
end

# From Jin et al (2001):
#   At each time-step, we choose $n_p * r_0$ pairs of vertices
#   uniformly at random from the network to meet. If a pair
#   meet who do not have a pre-existing connection, and if neither
#   of them has the maximum number of z* connections then a new
#   connection is established between them.
function meet_randomly(jingen)
    for _ in 1:jingen.n_meetings
        edge = random_edge_as_indices(jingen.params.n)
        if connection_can_be_made(jingen, edge)
            add_edge!(jingen, edge)
        end
    end
end

# From Jin et al (2001):
#   At each time-step, we choose $n_m * r_1$ vertices at random, with
#   probabilities proportional to $z_i(z_i-1)$. For each vertex chosen
#   we randomly choose one pair of its neighbors to meet, and establish
#   a new connection between them if they do not have a pre-existing
#   connection and if neither of them already has the maximum number
#   z* of connections.
function introduce_friends(jingen)
    num_mutual = @compat round(Int, 0.5 * sum(jingen.mutual_meet_factors))
    n_introductions = @compat round(Int, num_mutual * jingen.params.r_1)

    for _ in 1:n_introductions
        ego = rand(jingen.mutual_meet_factors)
        alters = jingen.adj_list[ego]

        # You can't introduce your friends if you only have 1 or none.
        # But, this property is true by the random proportional selection
        # algorithm. k * (k - 1) = 0 for k = 0 and k = 1. I don't have
        # to check it in production code.
        assert(length(alters) >= 2)

        edge = random_pair_from_list(alters)
        if connection_can_be_made(jingen, edge)
            add_edge!(jingen, edge)
        end
    end
end

function lose_touch(jingen)
    n_lost = @compat round(Int, (length(jingen.edges) * jingen.params.γ))
    for _ in 1:n_lost
        ego = rand(jingen.degrees)
        alter = rand(jingen.adj_list[ego])
        edge = ego < alter ? (ego, alter) : (alter, ego)
        remove_edge!(jingen, edge)
    end
end

function noop(i, jingen)
end

function generate(jingen::JinGenBuilder, iterations=30000, callback=noop)
    for i in 1:iterations
        meet_randomly(jingen)
        introduce_friends(jingen)
        if i > jingen.params.barrier_iters
            lose_touch(jingen)
        end

        callback(i, jingen)
    end

    return jingen.edges
end

# See: https://github.com/jbn/agents/blob/6eef634f80fa048e4c2fc5c62d85b384545b0e02/Repast%20J/demos/src/uchicago/src/repastdemos/jinGirNew/JinGirNewModelII.java
immutable FlawedJinGen
    builder::JinGenBuilder
end

function lose_touch(flawed::FlawedJinGen)
    jingen = flawed.builder
    n_lost = @compat round(Int, (length(jingen.edges) * jingen.params.γ))
    for _ in 1:n_lost
        ego = rand(jingen.mutual_meet_factors)  # This is wrong
        alter = rand(jingen.adj_list[ego])
        edge = ego < alter ? (ego, alter) : (alter, ego)
        remove_edge!(jingen, edge)
    end
end

function generate(flawed::FlawedJinGen, iterations=30000, callback=noop)
    jingen = flawed.builder
    for i in 1:iterations
        meet_randomly(jingen)
        introduce_friends(jingen)
        if i > jingen.params.barrier_iters
            lose_touch(flawed)
        end

        callback(i, jingen)
    end

    return jingen.edges
end

end


text:
# The List implementation is inspired by https://github.com/JuliaCollections/DataStructures.jl
export List, Nil, Cons, list, cons, nil, head, tail
abstract type List{T} end
struct Nil{T} <: List{T} end
struct Cons{T} <: List{T}
  head::T
  tail::List{T}
end

nil(T) = Nil{T}()
nil() = nil(Any)
cons(h, t::List{T}) where {T} = Cons{T}(h, t)
list() = nil()
Base.isempty(lst::Nil) = true
Base.isempty(lst::Cons) = false
Base.getindex(lst::Nil, i) = throw(BoundsError(lst, i))
Base.getindex(lst::Cons, i) = i == 1 ? lst.head : getindex(lst.tail, i-1)

list(elts...) =
  let l = nil()
    for i = length(elts):-1:1
      l = cons(elts[i], l)
    end
    l
  end

list(elts::T...) where {T} =
  let l = nil(T)
    for i = length(elts):-1:1
      l = cons(elts[i], l)
    end
    l
  end
head(x::Cons) = x.head
tail(x::Cons) = x.tail

import Base.==
==(x::Nil, y::Nil) = true
==(x::Cons, y::Cons) = (x.head == y.head) && (x.tail == y.tail)

show(io::IO, lst::Nil) = print(io, "list()")
show(io::IO, lst::Cons) =
  begin
    print(io, "list(")
    show(io, head(lst))
    for e in tail(lst)
      print(io, ", ")
      show(io, e)
    end
    print(io, ")")
  end

Base.length(l::Nil) = 0
Base.length(l::Cons) =
  let n = 0
    for i in l
      n += 1
    end
    n
  end

Base.map(f::Base.Callable, lst::List) = list((f(e) for e in lst)...)
Base.filter(f::Function, lst::List) = list((e for e in lst if f(e))...)

Base.cat() = list()
Base.cat(lst::List, lsts::List...) =
  let T = typeof(lst).parameters[1]
    n = length(lst)
    for l in lsts
      T2 = typeof(l).parameters[1]
      T = typejoin(T, T2)
      n += length(l)
    end
    elems = Vector{T}(undef, n)
    i = 1
    for e in lst
      elems[i] = e
      i += 1
    end
    for lst in lsts
      for e in lst
      elems[i] = e
      i += 1
      end
    end
    let l = nil(T)
      for i = i-1:-1:1
        l = cons(elems[i], l)
      end
      l
    end
  end

Base.iterate(l::List, ::Nil) = nothing
Base.iterate(l::List, state::Cons = l) = state.head, state.tail

# Lists can be converted to Arrays

Base.convert(::Type{Array{S,1}}, l::List{T}) where {S, T <: S} = collect(T, l)

export random, random_range, set_random_seed,
       division, map_division

#previous_random::Int = 12345
previous_random = 12345

set_random_seed(v::Int) =
  global previous_random = v

next_random(previous_random::Int) =
  let test = 16807*rem(previous_random,127773) - 2836*div(previous_random,127773)
    if test > 0
      if test > 2147483647
        test - 2147483647
      else
        test
      end
    else
      test + 2147483647
    end
  end

next_random!() =
  begin
    global previous_random = next_random(previous_random)
    previous_random
  end

random(x::Int) = rem(next_random!(), x)

random(x::Real) = x*next_random!()/2147483647.0

random_range(x0, x1) =
  if x0 == x1
    x0
  else
    x0 + random(x1 - x0)
  end

export RGB, rgb, rgb_radiance
const rgb = RGB

rgb_radiance(c::RGB) = 0.265*red(c)+0.67*green(c)+0.065*blue(c)

required() = error("Required parameter")

division(t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    collect(include_last ? iter : take(iter, n))
  end

# Generic objects are processed using map_division
division(obj::Any, n::Real) = map_division(identity, obj, n)

map_division(f, t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    map(f, include_last ? iter : take(iter, n))
  end

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv),
               u0, u1, nu, include_last_u)

map_division(f, u0, u1, nu::Real, v0, v1, nv::Real, include_last_v::Bool=true) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu)

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real, include_last_v::Bool) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu, include_last_u)

# Grasshopper compatibility

export series, crossref, remap, cull

series(start::Real, step::Real, count::Int) =
  range(start, step=step, length=count)

export crossref_holistic
crossref_holistic(arr1, arr2) =
  vcat([arr1[i] for i in range(1, stop=length(arr1)) for j in arr2]...),
  vcat([arr2 for i in arr1]...)

crossref(as, bs) = [(a, b) for a in as, b in bs]

remap(in, (min_in, max_in), (min_out, max_out)) =
  min_out + (max_out-min_out)/(max_in-min_in)*(in-min_in)

cull(template, as) =
  [a for (a, t) in zip(as, cycle(template)) if t]

# To create paths from paths

path_replace_suffix(path::String, suffix::String) =
  let (base, old_suffix) = splitext(path)
    base * suffix
  end

#

replace_in(expr::Expr, replacements) =
    if expr.head == :.
        Expr(expr.head,
             replace_in(expr.args[1], replacements), expr.args[2])
    elseif expr.head == :quote
        expr
    else
        Expr(expr.head,
             map(arg -> replace_in(arg, replacements), expr.args) ...)
    end
replace_in(expr::Symbol, replacements) =
    get(replacements, expr, esc(expr))
replace_in(expr::Any, replacements) =
    expr

showit(s, a) = begin
    print(s)
    print(":")
    println(a)
    a
end

#
function process_named_params(def)
  call, body = def.args[1], def.args[2]
  name, params = call.args[1], call.args[2:end]
  idx = findfirst(p->p.head==:(::), params)
  mand, opts = isnothing(idx) ? ([], params) : (esc.(params[1:idx]), params[idx+1:end])
  opt_names = map(opt -> opt.args[1].args[1], opts)
  opt_types = map(opt -> esc(opt.args[1].args[2]), opts)
  opt_inits = map(opt -> opt.args[2], opts)
  opt_renames = map(Symbol ∘ string, opt_names)
  opt_replacements = Dict(zip(opt_names, opt_renames))
  mk_param(name,typ,init) = Expr(:kw, Expr(:(::), name, typ), init)
  opt_params = map(mk_param, opt_renames, opt_types, map(init -> replace_in(init, opt_replacements), opt_inits))
  key_params = map(mk_param, opt_names, opt_types, opt_renames)
  func_name = esc(name)
  :($(func_name)($(mand...), $(opt_params...); $(key_params...)) = $(esc(body)))
end

macro named_params(def)
  process_named_params(def)
end


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

"""
res=metabolicEP(S,b,lb,ub,...)

The output in res is of type `EPOut`: there are several fields:
-   `μ::Vector`: A parameter linked to the mean of the posterior probability
-   `σ::Vector`: A parameter linked to the std  of the posterior probability
-   `av::Vector`: The mean posterior probability
-   `va::Vector`: The variance of the posterior probability
-   `sol::EPFields`: The internal field status. From this value we can restart the sampling from a specific state.
-   `status::Symbol`: either ``:converged`` or ``:unconverged``.

Input (required)
----
- `S`: MxN matrix (either sparse or dense) please note that if you input a dense version, the algorithm is slighlty more efficient. Dense matrices can be create from sparse ones with `Matrix(S)`.
- `b`: a vector of M intakes/uptakes
- `lb`: a vector of length N of lower bounds.
- `ub`: a vector of length N of upper bounds.

Input (optional arguments).
----
- `alpha` (inverse temperature::``Real``): default 10^7;  the zero temperature algorithm is run setting ``beta=Inf``.
- `beta_vec` (MaxEnt distribution inverse temperature vector:: ``Vector{T}``): default ``T[]``
- `verbose` (``true`` or ``false``): default ``true``
- `damp` (∈ (0,1) newfield = damp * oldfield + (1-damp)* newfield): default 0.9
- `epsconv` (convergence criterion): default 1e-6
- `maxiter` (maximum number of iterations): default 2000
- `maxvar`  (threshold on maximum variance): default 1e50
- `minvar`  (threshold on minimum variance): default 1e-50
- `solution` (start from solution. Is of type ``EPOut``): default: ``nothing``
- `expval` (fix to posterior probability of mean and/or variance to values): default ``nothing``. expval can be either at ``Tuple{Float64,Float64,Int}`` or a ``Vector{Tuple{Float64,Float64,Int}}``. Values can be fixed as``expval=(0.2,0.4,4)`` meaning that for flux index 4 the mean is set to 0.2 and the variance to 0.4. Fixing more values ``expval=[(0.2, 0.3, 4), (0.4, nothing, 5)]``: in this case, we fix the posterior of flux 4 to 0.2 (mean) and 0.3 (variance), while for flux 5 we fix the mean to 0.4 and we keep the variance free.
"""

function maxent_ep(
        S::AbstractArray{T,2}, b::Array{T,1}, lb::Array{T,1}, ub::Array{T,1};
        alpha::Real=Inf,                             # inverse temperature
        beta_vec::AbstractVector{T} = T[],           # maxent inverse temperature vector
        verbose::Bool=true,                          # output verbosity
        damp::Real=0.9,                              # damp ∈ (0,1) newfield = damp * oldfield + (1-damp)* newfield
        epsconv::Real=1e-6,                          # convergence criterion
        maxiter::Int=2000,                           # maximum iteration count
        maxvar::Real=1e50,                           # maximum numerical variance
        minvar::Real=1e-50,                          # minimum numerical variance
        solution::Union{EPOut{T}, Nothing}=nothing,  # start from a solution
        expval=nothing,                              # fix posterior probability experimental values for std and mean
        iter0 = 0,                                   # the started iteration count
        drop_epfields = false,                       # if the final EPOut object will export the epfields
        # callbacks
        oniter::Function = (it, epmodel) -> (false, nothing)
    ) where T<:Real

    epmodel = EPModel(S, b, lb, ub; alpha, beta_vec, solution, expval)
    epout = converge_ep!(epmodel; verbose, damp, epsconv,
        maxiter, maxvar, minvar, iter0, drop_epfields, oniter
    )
end

maxent_ep(metnet::MetNet; kwargs...) =
    maxent_ep(metnet.S, metnet.b, metnet.lb, metnet.ub; kwargs...)


text:
@testset "Basis" begin
    @variables u[1:3]
    @parameters w[1:2] t
    h = [u[1]; u[2]; cos(w[1]*u[2]+w[2]*u[3]); u[3]+u[2]]
    h_not_unique = [u[1]; u[1]; u[1]^1; h]
    basis = Basis(h_not_unique, u, parameters = w, iv = t)

    @test isequal(variables(basis), u)
    @test isequal(parameters(basis), w)
    @test isequal(independent_variable(basis), t)
    @test free_parameters(basis) == 6
    @test free_parameters(basis, operations = [+, cos]) == 7
    @test DataDrivenDiffEq.count_operation((ModelingToolkit.Constant(1) + cos(u[2])*sin(u[1]))^3, [+, cos, ^, *]) == 4

    basis_2 = unique(basis)
    @test isequal(basis, basis_2)
    @test size(basis) == size(h)
    @test basis([1.0; 2.0; π], [0. 1.]) ≈ [1.0; 2.0; -1.0; π+2.0]
    @test size(basis) == size(basis_2)
    push!(basis_2, sin(u[2]))
    @test size(basis_2)[1] == length(h)+1

    basis_3 = merge(basis, basis_2)
    @test size(basis_3) == size(basis_2)
    @test isequal(variables(basis_3), variables(basis_2))
    @test isequal(parameters(basis_3), parameters(basis_2))

    merge!(basis_3, basis)
    @test basis_3 == basis_2

    push!(basis, u[3]+u[2])
    unique!(basis)
    @test size(basis) == size(h)

    @variables a
    g = [u[1]; u[3]; a]
    basis = Basis(g, [u; a])
    @test basis([1; 2; 3; 4]) == [1; 3; 4]
    g = [u[1]; u[3]; u[2]]
    basis = Basis(g, u, parameters = [])
    X = ones(Float64, 3, 10)
    X[1, :] .= 3*X[1, :]
    X[3, :] .= 5*X[3, :]
    # Check the array evaluation
    @test basis(X) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0] * X
    Y = similar(X)
    basis(Y, X)
    @test Y ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0] * X
    f = jacobian(basis)
    @test f([1;1;1], [0.0 0.0], 0.0) ≈ [1.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0]
    @test_nowarn sys = ODESystem(basis)
    @test_nowarn [xi for xi in basis]
    @test_nowarn basis[2:end]; basis[2]; first(basis); last(basis); basis[:]

    @variables u[1:3] t

    g = [u[2]; -sin(u[1])*exp(-t); u[2]+u[3]]
    basis = Basis(g, u, iv = t)
    @test_nowarn ODESystem(basis)

    f(u, p, t) = [u[3]; u[2]*u[1]; p[1]*sin(u[1])*u[2]; p[2]*t]
    b = Basis(f, u, parameters = w, iv = t)
    @test f([1; 2; 3], [2; 0], 3.0) ≈ b([1; 2; 3], [2; 0], 3.0)
    @test_throws AssertionError ODESystem(b)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/multiplication.jl
# Defines multiplication operator.
#############################################################################

# Nonsmooth multiplication kernel definition
@inline function mul1_u1pos_u2pos(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
    xLc = z.lo
    xUc = z.hi
  	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
  	if cv1 > cv2
    	cv = cv1
    	cv_grad = x2.Intv.hi*x1.cv_grad
  	else
    	cv = cv2
    	cv_grad = x2.Intv.lo*x1.cv_grad
  	end

  	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
  	if cc1 < cc2
    	cc = cc1
    	cc_grad = x2.Intv.lo*x1.cc_grad
  	else
    	cc = cc2
    	cc_grad = x2.Intv.hi*x1.cc_grad
  	end
  	return MCNoGrad(cv, cc, z, cnst)
end
@inline function mul1_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
    if cv1 > cv2
        cv = cv1
        cv_grad = x2.Intv.hi*x1.cv_grad
    else
        cv = cv2
        cv_grad = x2.Intv.lo*x1.cc_grad
    end

    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
    if cc1 < cc2
        cc = cc1
        cc_grad = x2.Intv.lo*x1.cv_grad
    else
        cc = cc2
        cc_grad = x2.Intv.hi*x1.cc_grad
    end
    return MCNoGrad(cv, cc, z, cnst)
end
@inline function mul1_u1mix_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
    cv = cv1
    if cv1 > cv2
        cv = cv1
        cv_grad = x2.Intv.hi*x1.cv_grad
    else
        cv = cv2
        cv_grad = x2.Intv.lo*x1.cc_grad
    end
    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
    cc = cc1
    if cc1 < cc2
        cc = cc1
        cc_grad = x2.Intv.lo*x1.cv_grad
    else
        cc = cc2
        cc_grad = x2.Intv.hi*x1.cc_grad
    end
    return MCNoGrad(cv, cc, z, cnst)
end
@inline function mul2_u1pos_u2pos(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64})
	xLc = z.lo
	xUc = z.hi
	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if cv1 > cv2
		cv = cv1
		cv_grad = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.lo*x2.cv_grad
	end
	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if cc1 < cc2
		cc = cc1
		cc_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end
	return MCNoGrad(cv, cc, z, cnst)
end
@inline function mul2_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
    if cv1 > cv2
        cv = cv1
		cvt2 = x1.Intv.hi
		cv2_is_cvg = 1
    else
        cv = cv2
		cvt2 = x1.Intv.lo
		cv2_is_cvg = 0
    end
	cvt1 = 0.0
	cv1_is_cvg = 2

    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
    if cc1 < cc2
        cc = cc1
		cct2 = x1.Intv.hi
    else
        cc = cc2
		cct2 = x1.Intv.lo
    end
	cct1 = 0.0
	cc1_is_ccg = 2
	cc2_is_ccg = 1

	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
end
@inline function mul2_u1mix_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64})
	xLc = z.lo
	xUc = z.hi
  	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
	if cv1 > cv2
		cv = cv1
		cvt1 = x2.Intv.hi
		cvt2 = x1.Intv.hi
		cv1_is_cvg = 1
		cv2_is_cvg = 1
	else
		cv = cv2
		cvt1 = x2.Intv.lo
		cvt2 = x1.Intv.lo
		cv1_is_cvg = 0
		cv2_is_cvg = 0
	end

	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
	if cc1 < cc2
		cc = cc1
		cct1 = x2.Intv.lo
		cct2 = x1.Intv.hi
		cc1_is_ccg = 0
		cc2_is_ccg = 1
	else
		cc = cc2
		cct1 = x2.Intv.hi
		cct2 = x1.Intv.lo
		cc1_is_ccg = 1
		cc2_is_ccg = 0
	end

	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
end
@inline function mul3_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64})
	xLc = z.lo
	xUc = z.hi
    cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if cv1 > cv2
		cv = cv1
		cvt1 = x2.Intv.hi
		cvt2 = x1.Intv.hi
		cv1_is_cvg = 1
	else
		cv = cv2
		cvt1 = x2.Intv.lo
		cvt2 = x1.Intv.lo
		cv1_is_cvg = 0
	end
	cv2_is_cvg = 1

	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if cc1 < cc2
		cc = cc1
		cct1 = x2.Intv.lo
		cct2 = x1.Intv.hi
		cc1_is_ccg = 0
	else
		cc = cc2
		cct1 = x2.Intv.hi
		cct2 = x1.Intv.lo
		cc1_is_ccg = 1
	end
	cc2_is_ccg = 1
	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
end

@inline function multiply_STD_NS(t::NS, x1::MCNoGrad, x2::MCNoGrad, y::Interval{Float64})
	if x2.Intv.lo >= 0.0
    	x2.cnst && (return mul1_u1pos_u2pos(t, x1, x2, y, x1.cnst))
    	x1.cnst && (return mul1_u1pos_u2pos(t, x2, x1, y, x2.cnst))
  		return mul2_u1pos_u2pos(t, x1, x2, y)
	elseif x2.Intv.hi <= 0.0
		return -mult_kernel(t, x1, -x2, -y)
	else
    	x2.cnst && (return mul1_u1pos_u2mix(t, x1, x2, y, x1.cnst))
    	x1.cnst && (return mul2_u1pos_u2mix(t, x1, x2, y, x2.cnst))
	end
	mul3_u1pos_u2mix(t, x1, x2, y)
end

@inline function mult_kernel(t::NS, x1::MCNoGrad, x2::MCNoGrad, y::Interval{Float64})
	isone(x1) && (return x2)
	isone(x2) && (return x1)
	if x1.Intv.lo >= 0.0
		return multiply_STD_NS(t, x1, x2, y)
	elseif x1.Intv.hi <= 0.0
		(x2.Intv.lo >= 0.0) && (return -mult_kernel(t, -x1, x2, -y))
	    (x2.Intv.hi <= 0.0) && (return mult_kernel(t, -x1, -x2, y))
		return -mult_kernel(t, x2, -x1, -y)
	elseif x2.Intv.lo >= 0.0
		return mult_kernel(t, x2, x1, y)
	elseif x2.Intv.hi <= 0.0
		return -mult_kernel(t, -x2, x1, -y)
	else
    	x2.cnst && (return mul1_u1mix_u2mix(t, x1, x2, y, x1.cnst))
    	x1.cnst && (return mul1_u1mix_u2mix(t, x2, x1, y, x2.cnst))
	end
	return mul2_u1mix_u2mix(t, x1, x2, y)
end

@inline function *(t::NS, x1::MCNoGrad, x2::MCNoGrad)
	mult_kernel(t, x1, x2, x1.Intv*x2.Intv)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#!/usr/bin/env julia

using DifferentialEquations, Test, SafeTestsets
@time begin
@time @safetestset "Default Discrete Algorithm" begin include("default_discrete_alg_test.jl") end
@time @safetestset "Default ODE Algorithm" begin include("default_ode_alg_test.jl") end
@time @safetestset "Default Steady State Algorithm" begin include("default_steady_state_alg_test.jl") end
@time @safetestset "Default SDE Algorithm" begin include("default_sde_alg_test.jl") end
@time @safetestset "Default RODE Algorithm" begin include("default_rode_alg_test.jl") end
@time @safetestset "Default DDE Algorithm" begin include("default_dde_alg_test.jl") end
@time @safetestset "Default DAE Algorithm" begin include("default_dae_alg_test.jl") end
@time @safetestset "Default BVP Algorithm" begin include("default_bvp_alg_test.jl") end
end

max_stars_repo_name: ChrisRackauckas/DiffEq

max_issues_repo_name: ChrisRackauckas/DiffEq


text:
"""
This `AeolianSources` module is to create aeolian sources for use in AIBECS-generated models.
"""
module AeolianSources

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using NCDatasets
using Unitful               # for units
using Reexport
using Interpolations
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

# Aeolian source URLs
url() = "http://www.geo.cornell.edu/eas/PeoplePlaces/Faculty/mahowald/dust/Chienetal2016/post.aerosols.2x2.seasonal.nc"

# Hashes
sha() = "33b0036bcca87019e0636ec2d59e910be989693819a2d57f2e4ab053a564101f" # On 2 June 2020

# Create registry entry for the aeolian source in BSON format
function register_AeolianSource()
    register(
        DataDep(
            "AIBECS-Chienetal2016",
            """
            References:
            - $(citation())
            """,
            url(),
            sha(),
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load

Returns the 2D deposition fields.
"""
function load()
    register_AeolianSource()
    nc_file = @datadep_str string("AIBECS-Chienetal2016/", "post.aerosols.2x2.seasonal.nc")
    s_A_2D = Dataset(nc_file, "r") do ds
        Dict(
            :lat => ds["lat"][:],
            :lon => ds["lon"][:],
            :Fires                    => ds["dep"][:,:,:,1],
            :Biofuels                 => ds["dep"][:,:,:,2],
            :FossilFuel               => ds["dep"][:,:,:,3],
            :Dust                     => ds["dep"][:,:,:,4],
            :Seasalts                 => ds["dep"][:,:,:,5],
            :PrimaryBiogenicParticles => ds["dep"][:,:,:,6],
            :Volcanoes                => ds["dep"][:,:,:,7]
        )
    end # ds is closed
    @info """You are about to use the Chien et al. (2016) data for aeolian deposition.
          If you use it for research, please cite:

          - $(citation())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "Chien_etal_2016" key.)
          """
    return s_A_2D
end

citation() = "Chien, C.-T., K. R. M. Mackey, S. Dutkiewicz, N. M. Mahowald, J. M. Prospero, and A. Paytan (2016), Effects of African dust deposition on phytoplankton in the western tropical Atlantic Ocean off Barbados, Global Biogeochem. Cycles, 30, doi:10.1002/2015GB005334."

end # end module

export AeolianSources


text:
# Package objects provide a _slightly_ higher-level API for dealing with
# installing/uninstalling tarballs to a Prefix.
export BinaryPackage, install, uninstall, satisfied

"""
A `BinaryPackage` collects all the information needed to download and install a
tarball containing binary objects; it has the `url` to download form, the
`hash` to verify package integrity, a list of `products` to check for proper
functioning after installation, and a list of `dependencies` that must also be
installed before this package can be used.

There exist `install()`, `uninstall()` and `satisfied()` methods for
`BinaryPackage` objects, similar to the lower-level versions that take direct
`url` and `hash` arguments.
"""
struct BinaryPackage
    url::String
    hash::String
    platform::Platform
    products::Vector{Product}
    dependencies::Vector{BinaryPackage}

    function BinaryPackage(url::AbstractString,
                           hash::AbstractString,
                           platform::Platform,
                           products::Vector{Product}=Product[],
                           dependencies::Vector{BinaryPackage}=BinaryPackage[])
        return new(url, hash, platform, products, dependencies)
    end
end

function pkg_name(pkg::BinaryPackage)
    name = basename(pkg.url)
    if endswith(name, ".tar.gz")
        return name[1:end-7]
    end
    return name
end

"""
`install(pkg::BinaryPackage; verbose::Bool = false, kwargs...)`

A thin wrapper over the main `install(url, path; ...)` method.  Installs all
of `pkg`'s dependencies, then installs `pkg`, but only of it is not already
satisfied.
"""
function install(pkg::BinaryPackage; verbose::Bool = false, kwargs...)
    name = pkg_name(pkg)

    # If we are already satisfied, don't do nuthin'
    if satisfied(pkg)
        if verbose
            info("Not installing $(name), as it is already satisfied")
        end
        return true
    end

    # Begin by installing all the dependencies if they are not already
    for dep in pkg.dependencies
        # TODO: We may want to handle this through `Pkg3` operations
        install(dep; verbose=verbose, kwargs...)
    end

    # Finally, install ourselves
    install(pkg.url, pkg.hash; verbose=verbose, kwargs...)

    # Check to see if we are actually satisfied
    if !satisfied(pkg; verbose=verbose)
        warn("$(name) did not satisfy itself after installation!")
        return false
    end

    return true
end

"""
manifest_path(pkg::BinaryPackage; prefix::Prefix = global_prefix(),
                                  verbose::Bool = false)

Discovers the manifest path for the given `BinaryPackage` within the given
`Prefix`.  First attempts to guess from the `url` what the manifest file would
have been named, if that doesn't work, will search for manifests that contain
any of the `products` that are within `pkg`.  If neither approach works, throws
an error.
"""
function manifest_path(pkg::BinaryPackage; prefix::Prefix = global_prefix(),
                                           verbose::Bool = false)
    name = pkg_name(pkg)
    # First, see if we can auto-guess the manifest file path:
    manifest_path = manifest_from_url(pkg.url, prefix=prefix)
    if isfile(manifest_path)
        if verbose
            info("Correctly auto-guessed manifest path $(manifest_path)")
        end
        return manifest_path
    end

    if verbose
        info("Could not auto-guess manifest path for $(name)")
    end

    # Otherwise, let's try to guess from our products
    if isempty(pkg.products)
        msg = """
        Cannot find manifest path for package $(name) with unguessable manifest
        file and no products.
        """
        error(replace(strip(msg),"\n", " "))
    end

    for product in pkg.products
        product_path = locate(product, platform=pkg.platform)
        if product_path != nothing
            try
                manifest_path = manifest_for_file(product_path; prefix=prefix)
                relmani = relpath(manifest_path, prefix.path)
                relprod = relpath(product_path, prefix.path)
                info("Found $(relmani) for product $(relprod)")
                return manifest_path
            end
        end
    end

    error("Cannot find manifest path for package $(name)")
end

"""
uninstall(pkg::BinaryPackage; prefix::Prefix = global_prefix,
                              verbose::Bool = false)

Uninstall `pkg` from the given `prefix` by automatically determining the
manifest path created when the package was installed.  Throws an error if
this `pkg` was not installed in the first place.
"""
function uninstall(pkg::BinaryPackage; prefix::Prefix = global_prefix,
                                       verbose::Bool = false)
    # Find the manifest path for this pkg, then uninstall it
    manipath = manifest_path(pkg; prefix=prefix, verbose=verbose)
    uninstall(manipath; verbose=verbose)
end

"""
`satisfied(pkg::BinaryPackage; verbose::Bool = false)`

Returns `true` if all products defined within `pkg` are satisfied.
"""
function satisfied(pkg::BinaryPackage; verbose::Bool = false)
    s = p -> satisfied(p; platform=pkg.platform, verbose=verbose)
    return all(s(p) for p in pkg.products)
end


text:
struct NormPSU{T,D} <: ObjectiveFunction
    dims::Dims{D}
    δt::Float64 # timestep (fixed)
    Ut::Matrix{ComplexF64} # target unitary
    Hd::Matrix{T} # drift Hamiltonian
    Hc::Vector{Matrix{T}} # control Hamiltonian(s)
    u_last::Vector{Float64} # last control amplitudes, ordered by control, and then time
    # Propagator and eigendecompotion storage
    U::Vector{Matrix{ComplexF64}} # individual propagators
    X::Vector{Matrix{ComplexF64}} # forward cumulative propagators
    P::Vector{Matrix{ComplexF64}} # backwards cumulative propagators
    D::Vector{Vector{Float64}} # eigenvalues associated to each propagator
    V::Vector{Matrix{T}} # eigenvector matrix associated to each propagator
    # Workspace variables
    H::Hermitian{T,Matrix{T}} # temporary storage for full Hamiltonian
    A::Matrix{ComplexF64} # temporary
    Jkj::Matrix{ComplexF64} # temporary
    cisDj::Vector{ComplexF64} # temporary
end

NormPSU(dims,δt,Ut,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj) = NormPSU{eltype(H),length(dims)}(dims,δt,Ut,Hd,Hc,u_last,U,X,P,D,V,H,A,Jkj,cisDj)

function NormPSU(Ut::Operator,Hd::Operator,Hc::Vector{<:Operator},t::Real,n::Integer)
    # TODO: Prob should do some error checking
    N = prod(dims(Hd))
    m = length(Hc)
    # Make sure we pass dense operators
    Ut_d = Array(Ut)
    Hd_d = Array(Hd)
    Hc_d = Array.(Hc)
    # Storage for last control ampitudes, NaN for first run
    u_last = fill(NaN64,m*n)
    # Generate cache for various objects
    H = Hermitian(zeros(promote_eltype(typeof(Hd_d),eltype(Hc_d)),N,N))
    A = Matrix{ComplexF64}(undef,N,N)
    U = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    X = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    P = [Matrix{ComplexF64}(undef,N,N) for i=1:n]
    copyto!(P[end],Ut_d)
    # For the exact derivative we need to store eigenvectors and eigenvalues
    D = [Vector{Float64}(undef,N) for i=1:n] # eigenvalues are always real
    V = [similar(H.data) for i=1:n]
    # More temporary storage
    Jkj = Matrix{ComplexF64}(undef,N,N)
    cisDj = Vector{ComplexF64}(undef,N)
    return NormPSU(dims(Hd),t/n,Ut_d,Hd_d,Hc_d,u_last,U,X,P,D,V,H,A,Jkj,cisDj)
end

function objective(O::NormPSU,u)
    N = size(O.Ut,1); Uf = O.X[end]
    # Calculate forward propagators
    calc_fprops!(O,u)
    # Calculate PSU (projective special unitary) norm infidelity:
    # fₑ = 1 - |⟨Ut,Uf⟩|/N
    return 1 - abs(inner(O.Ut,Uf))/N
end

function gradient!(O::NormPSU,fp,u)
    n = length(O.U); m = length(O.Hc); N = size(O.Ut,1); Uf = O.X[end]
    # Calculate forward and backward propagators
    calc_fprops!(O,u)
    calc_bprops!(O)
    # Calculate exact derivative of fidelity error function:
    # ∂fₑ/∂uₖⱼ = -1/2*⟨Pⱼ,∂Uⱼ/∂uₖⱼ*Xⱼ₋₁⟩*cis(∠⟨Xⱼ,Pⱼ⟩)/N + c.c.
    #          = -Re(⟨Pⱼ,Jₖⱼ*Xⱼ₋₁⟩*cis(∠⟨Uf,Ut⟩))/N where Jₖⱼ = ∂Uⱼ/∂uₖⱼ
    a = normalize(inner(Uf,O.Ut))
    for j = 1:n
        O.cisDj .= cis.(O.D[j])
        for k = 1:m
            Jmat!(O,k,j)
            j==1 ? copyto!(O.A,O.Jkj) : mul!(O.A,O.Jkj,O.X[j-1])
            fp[(k-1)*n+j] = -real(inner(O.P[j],O.A)*a)/N
        end
    end
    return fp
end


text:
using DelaySSAToolkit
using Catalyst
using Statistics
using Test
using BenchmarkTools
using DataFrames
using CSV
reltol = 0.1
@parameters t
@variables X(t)
burst_sup = 30
a, b = [0.0282, 3.46]
rxs = [Reaction(a*b^i/(1+b)^(i+1),nothing,[X],nothing,[i]) for i in 1:burst_sup]
rxs = vcat(rxs)
@named rs = ReactionSystem(rxs,t,[X],[])

# convert the ReactionSystem to a JumpSystem
jumpsys = convert(JumpSystem, rs, combinatoric_ratelaws=false)

u0 = [0]
de_chan0 = [[]]
tf = 200.
tspan = (0,tf)
dprob = DiscreteProblem(jumpsys,u0,tspan)
τ = 130.

delay_trigger = Dict([Pair(i, [1=>fill(τ, i)]) for i in 1:burst_sup])
delay_complete = Dict(1=>[1=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

algo_list = [DelayDirect(),DelayRejection(),DelayMNRM(),DelayDirectCR()]

jprob = DelayJumpProblem(jumpsys, dprob, algo_list[1], delayjumpset, de_chan0, save_positions=(false, false))
a=@benchmark solve(jprob, SSAStepper())

jprob = DelayJumpProblem(jumpsys, dprob, algo_list[2], delayjumpset, de_chan0, save_positions=(false, false))
a=@benchmark solve(jprob, SSAStepper())

jprob = DelayJumpProblem(jumpsys, dprob, algo_list[3], delayjumpset, de_chan0, save_positions=(false, false))
a=@benchmark solve(jprob, SSAStepper())

jprob = DelayJumpProblem(jumpsys, dprob,algo_list[4] , delayjumpset, de_chan0, save_positions=(false, false))
a=@benchmark solve(jprob, SSAStepper())

meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]
function testalgo(algo_list)
    for algo in algo_list
        #de_chan0 = [[]]
        jprob = DelayJumpProblem(jumpsys, dprob, algo, delayjumpset, de_chan0, save_positions=(false, false))
        a=@benchmark solve(jprob, SSAStepper())
        push!(meanlist,mean(a).time/1e9)
        push!(medianlist,median(a).time/1e9)
        push!(minlist,minimum(a).time/1e9)
        push!(maxlist,maximum(a).time/1e9)
        push!(stdlist,std(a).time/1e9)
    end
end

@time testalgo(algo_list)

df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("test/zxytest/results/bursty_save_positions=f.csv",df)

df=CSV.read("C:/Users/86158/Desktop/algotest/bursty_save_positions=f.csv",DataFrame)

medianlist=df.median

algo_name = ["Direct","Rejection","MNRM","DirectCR"]

meanvalue=[string(round(mt,digits=8),"s") for mt in meanlist]
medianvalue=[string(round(mt,digits=8),"s") for mt in medianlist]

using Plots
bar(algo_name,meanlist,legend=:false,title="bursty",ylabel="meantime")
scatter!(algo_name, 0.3e-5 .+ meanlist , markeralpha=0, series_annotations=meanvalue)

bar(algo_name,medianlist,legend=:false,title="bursty",ylabel="mediantime")
scatter!(algo_name, 0.3e-5 .+ medianlist , markeralpha=0, series_annotations=medianvalue)


text:
print("Enter latitude       => ")
lat = parse(Float64, readline(STDIN))
print("Enter longitude      => ")
lng = parse(Float64, readline(STDIN))
print("Enter legal meridian => ")
ref = parse(Float64, readline(STDIN))
println()

slat = sin(deg2rad(lat))
@printf "    sine of latitude:   %.3f\n" slat
@printf "    diff longitude:     %.3f\n" (lng - ref)

println("\nHour, sun hour angle, dial hour line angle from 6am to 6pm\n")

for h in -6:6
  hra = 15 * h
  hra -= lng - ref
  hla = rad2deg(atan(slat * tan(deg2rad(hra))))
  @printf "HR = %3d; HRA = %7.3f; HLA = %7.3f\n" h hra hla
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
 Tests for  https://link.springer.com/article/10.1007/s10543-007-0130-3 with test problems as in the paper.
 RS1, RS2
 and for https://www.sciencedirect.com/science/article/pii/S0377042706003906
 NON, COM, NON2
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
using DiffEqDevTools
#using DiffEqGPU

seed = 137475

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test non-commutative noise SDEs (iip)
"""

@info "Non-commutative noise"

u₀ = [0.1,0.1]
function f2!(du,u,p,t)
  du[1] = 5//4*u[2]-5//4*u[1]
  du[2] = 1//4*u[1]-1//4*u[2]
end
function g2!(du,u,p,t)
  du[1,1] = sqrt(3)/2*(u[1]-u[2])
  du[1,2] = 1//2*(u[1]+u[2])
  #du[2,1] = 0
  du[2,2] = u[1]
end
dts = 1 .//2 .^(4:-1:1)
tspan = (0.0,1.0)

h2(z) = z^2 # E(x_i) = 1/10 exp(1/2t) or E(x_1* x_2) = 1/100 exp(2t)

prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2))
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h2(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(3e7)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,RS1(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(2)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RS1:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,RS2(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(2)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("RS2:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,NON(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(2)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("NON:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,NON2(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(2)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("NON2:", sim.𝒪est[:weak_final])

numtraj = Int(2e7)
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,COM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(2)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("COM:", sim.𝒪est[:weak_final])
# COM tests are passing; problem might be not hard enough..

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c07e59de-7de1-11eb-31ba-9febfeacb93f
using Distributions

# ╔═╡ 73290c6a-7dd7-11eb-363c-1ff5a2aad1f5
using Plots

# ╔═╡ c4145b48-7de1-11eb-059d-870e7d3004c5
begin
	μ = @bind mu html"<input type='range' min='-10.0' max='10.0' step='0.1' value='0.0'>"
	σ = @bind sigma html"<input type='range' min='1' max='10.0' step='0.1' value='1.0'>"
	md"""**distribution parameters:**
	$$\mu$$ $μ
	$$\sigma$$ $σ
	"""
end

# ╔═╡ b0fa3884-7de1-11eb-11e1-7dcc4792c795
x = [mu - 5 * sigma:0.1:mu + 5 * sigma;]

# ╔═╡ 8ad224ba-7e03-11eb-3390-fd30f91b85b6
md"""
 $$\mu$$ = $mu

 $$\sigma$$ = $sigma
"""

# ╔═╡ 019515ce-7de3-11eb-39ac-f5a83fb11fb7
md"""
Mean: $$\mu$$ = $mu

Variance: $$\sigma^2$$ = $(sigma^2)
"""

# ╔═╡ b940194a-7de2-11eb-2f76-ad1abefdb11d
N = Normal(mu, sigma^2);

# ╔═╡ b15f8620-7de2-11eb-323f-2f98407f01f3
plot(x, pdf.(N, x), label = "Normal( $mu, $(round(sigma^2, digits = 2))",
	linewidth = 2, 	ylims = (0, 0.45 / sigma^0.5))

# ╔═╡ Cell order:
# ╟─b0fa3884-7de1-11eb-11e1-7dcc4792c795
# ╟─c4145b48-7de1-11eb-059d-870e7d3004c5
# ╟─8ad224ba-7e03-11eb-3390-fd30f91b85b6
# ╟─019515ce-7de3-11eb-39ac-f5a83fb11fb7
# ╠═b940194a-7de2-11eb-2f76-ad1abefdb11d
# ╟─b15f8620-7de2-11eb-323f-2f98407f01f3
# ╟─c07e59de-7de1-11eb-31ba-9febfeacb93f
# ╟─73290c6a-7dd7-11eb-363c-1ff5a2aad1f5

max_stars_repo_name: roualdes/Normal

max_issues_repo_name: roualdes/Normal

max_forks_repo_name: roualdes/Normal


text:
using Test
using GeoPhyInv
using LinearMaps
using Misfits
using LinearAlgebra

# ### SeisInvExpt
pa=SeisInvExpt(FdtdBorn(), LS(), :pizza, rfields=[:P, :Vz]); # use linearized modeling attribute
F=LinearMap(pa); # generate LinearMap

# ### Linearity test
x1=randn(size(F,2)); # random input 1
x2=randn(size(F,2)); # random input 2
x12=x1.+x2; # sum

#-

d12=F*x12;
d1=F*x1; #
d2=F*x2;
d12new=d1.+d2;

#-

@test Misfits.error_squared_euclidean!(nothing, d12, d12new, nothing, norm_flag=true)<1e-25

# ### Adjoint test

x=randn(size(F,2)); # random x
y=randn(size(F,1)); # random y
a=LinearAlgebra.dot(y,F*x);
b=LinearAlgebra.dot(x,adjoint(F)*y);
c=LinearAlgebra.dot(x, transpose(F)*F*x);

#-
println("adjoint test: ", a, "\t", b)
@test isapprox(a,b,rtol=1e-5)

#-

println("x*transpose(F)*F*x be positive: ", c)
@test c>0.0


text:
isdefined(Turf, :BBox2D) || const BBox2D = Vector{Float64}(undef, 4)
isdefined(Turf, :BBox3D) || const BBox3D = Vector{Float64}(undef, 6)

"""
    bbox(geojson::T) where {T<:AbstractFeatureCollection}

Take a set of features, calculate the bbox of all input features, and returns a bounding box.
"""
function bbox(geojson::T) where {T<:AbstractFeatureCollection}
    result = [Inf, Inf, -Inf, -Inf]

    for feat in geojson.features
        geom = feat.geometry
        coords = geom.coordinates

        if geotype(geom) === :LineString || geotype(geom) === :MultiPoint
            for i in eachindex(coords)
                result[1] > coords[i][1] && (result[1] = coords[i][1])

                result[2] > coords[i][2] && (result[2] = coords[i][2])

                result[3] < coords[i][1] && (result[3] = coords[i][1])

                result[4] < coords[i][2] && (result[4] = coords[i][2])
            end
        elseif geotype(geom) === :Polygon || geotype(geom) === :MultiLineString
            for j in eachindex(coords)
                for i in eachindex(coords[j])
                    result[1] > coords[j][i][1] && (result[1] = coords[j][i][1])

                    result[2] > coords[j][i][2] && (result[2] = coords[j][i][2])

                    result[3] < coords[j][i][1] && (result[3] = coords[j][i][1])

                    result[4] < coords[j][i][2] && (result[4] = coords[j][i][2])
                end
            end
        elseif geotype(geom) === :Point
                result[1] > coords[1] && (result[1] = coords[1])

                result[2] > coords[2] && (result[2] = coords[2])

                result[3] < coords[1] && (result[3] = coords[1])

                result[4] < coords[2] && (result[4] = coords[2])
        elseif geotype(geom) === :MultiPolygon
            for i in eachindex(coords)
                for j in eachindex(coords[i])
                    for k in eachindex(coords[i][j])
                        result[1] > coords[i][j][k][1] && (result[1] = coords[i][j][k][1])

                        result[2] > coords[i][j][k][2] && (result[2] = coords[i][j][k][2])

                        result[3] < coords[i][j][k][1] && (result[3] = coords[i][j][k][1])

                        result[4] < coords[i][j][k][2] && (result[4] = coords[i][j][k][2])
                    end
                end
            end
        end
    end
    return result
end

"""
    bbox(geojson::T) where {T <: AbstractGeometry}

Take a GeoJSON Geometry and calculate its bounding box.
"""
function bbox(geojson::T) where {T <: AbstractGeometry}
    result = [Inf, Inf, -Inf, -Inf]

    coords = geojson.coordinates

    if geotype(geojson) === :LineString || geotype(geojson) === :MultiPoint
        for i in eachindex(coords)
            result[1] > coords[i][1] && (result[1] = coords[i][1])

            result[2] > coords[i][2] && (result[2] = coords[i][2])

            result[3] < coords[i][1] && (result[3] = coords[i][1])

            result[4] < coords[i][2] && (result[4] = coords[i][2])
        end
    elseif geotype(geojson) === :Polygon || geotype(geojson) === :MultiLineString
        for j in eachindex(coords)
            for i in eachindex(coords[j])
                result[1] > coords[j][i][1] && (result[1] = coords[j][i][1])

                result[2] > coords[j][i][2] && (result[2] = coords[j][i][2])

                result[3] < coords[j][i][1] && (result[3] = coords[j][i][1])

                result[4] < coords[j][i][2] && (result[4] = coords[j][i][2])
            end
        end
    elseif geotype(geojson) === :Point
            result[1] > coords[1] && (result[1] = coords[1])

            result[2] > coords[2] && (result[2] = coords[2])

            result[3] < coords[1] && (result[3] = coords[1])

            result[4] < coords[2] && (result[4] = coords[2])
    elseif geotype(geojson) === :MultiPolygon
        for i in eachindex(coords)
            for j in eachindex(coords[i])
                for k in eachindex(coords[i][j])
                    result[1] > coords[i][j][k][1] && (result[1] = coords[i][j][k][1])

                    result[2] > coords[i][j][k][2] && (result[2] = coords[i][j][k][2])

                    result[3] < coords[i][j][k][1] && (result[3] = coords[i][j][k][1])

                    result[4] < coords[i][j][k][2] && (result[4] = coords[i][j][k][2])
                end
            end
        end
    end

    return result
end

"""
     bbox(geojson::T) where {T<: AbstractFeature}

Take a Feature and return a bounding box around its geometry.
"""
bbox(geojson::T) where {T<: AbstractFeature} = bbox(geojson.geometry)

"""
    bbox_polygon(bbox::Vector{T}) where {T <: Real}

Take a bbox and return an equivalent Polygon.
"""
function bbox_polygon(bbox::Vector{T})::Polygon where {T <: Real}
    west = bbox[1]
    south = bbox[2]
    east = bbox[3]
    north = bbox[4]

    length(bbox) === 6 && throw(error("BBoxes with 6 positions are not supported."))

    lowLeft = [west, south]
    topLeft = [west, north]
    topRight = [east, north]
    lowRight = [east, south]

    return Polygon([[lowLeft, lowRight, topRight, topLeft, lowLeft]])
end


text:
# Return a Dict with plot settings. Easier to modify than a bunch of constants.
# - tickSize governs sizes of various fonts
function plot_settings()
    return Dict([
        :barLineWidth => 0,  # does not affect error bars
        :lineAlpha => nothing,  # for bar graphs
        :subPlotMargin => 3.0mm,
        :tickSize => 6,
        :markerStrokeWidth => 2,  # affects error bars
        :showScalarDev => true  # Show scalar deviations in plots?
    ])
end

bar_defaults() = Dict([
    :leg => :none,
    :xlabel => "",
    :ylabel => "",
    :margin => plot_setting(:subPlotMargin),
    :lw => 0
]);

grouped_bar_defaults() = Dict([
    :bar_position => :dodge,
    :leg => :bottomright,
    :labels => data_model_labels(),
    :markerstrokewidth => plot_setting(:markerStrokeWidth),
    :markerstrokecolor => :black,
    :linecolor => :black,
    :lw => 0
]);

"""
	$(SIGNATURES)

Retrieve a named plot setting.

# Example
```
plot_setting(:tickSize) == 7
```
"""
function plot_setting(sName :: Symbol)
    ps = plot_settings();
    @assert haskey(ps, sName)  "$sName not found."
    return ps[sName]
end

function plot_defaults()
    ps = plot_settings();
    # Backend
    gr();
    # theme(:wong);
    theme(:sand);
    # guidefont: x and y axis labels
    fontColor = :match;
    tickSize = ps[:tickSize];
    tickFont = (tickSize, fontColor);
    default(linewidth = 2, xtickfont = tickFont, ytickfont = tickFont,
        titlefontsize = tickSize + 2, legendfont = (tickSize, fontColor),
        guidefont = (tickSize, fontColor));
end

# ----------