text:
"""
    bea_query(url, querydict)

Internal method for handling requests to the BEA data API
"""
function bea_query(url, querydict)
    # Catch "broken pipe" IOError
    # https://github.com/stephenbnicar/BeaData.jl/issues/10
    # https://github.com/JuliaWeb/HTTP.jl/issues/382
    response = @repeat 3 try
        HTTP.request("GET", url, query = querydict)
    catch e
        @delay_retry if isa(e, Base.IOError) end
    end

    response_body = String(response.body)
    response_json = JSON.parse(response_body)

    # Check for bad requests
    if !haskey(response_json["BEAAPI"], "Results")
        if haskey(response_json["BEAAPI"], "Error")
            error(string(response_json["BEAAPI"]["Error"]["APIErrorDescription"]),
                " ", response_json["BEAAPI"]["Error"]["ErrorDetail"]["Description"])
        end
    elseif haskey(response_json["BEAAPI"]["Results"], "Error")
        error(response_json["BEAAPI"]["Results"]["Error"]["APIErrorDescription"])
    end

    return response_json["BEAAPI"]["Results"]
end

"""
    bea_datasets(;user_id::String = USER_ID) -> DataFrame

Return a `DataFrame` of names and descriptions for datasets accessible through the BEA data API.

Example:
```julia
# User ID stored in ~/.beadatarc
julia> bea_datsets();

julia> show(ans[1:3, :])
3×2 DataFrames.DataFrame
│ Row │ DatasetName        │ Description                          │
│     │ String             │ String                               │
├─────┼────────────────────┼──────────────────────────────────────┤
│ 1   │ NIPA               │ Standard NIPA tables                 │
│ 2   │ NIUnderlyingDetail │ Standard NI underlying detail tables │
│ 3   │ MNE                │ Multinational Enterprises            │
```
"""
function bea_datasets(;user_id = USER_ID)
    # These datasets are deprecated by BEA, though still show up in requests
    excluded_sets = ["RegionalData", "RegionalIncome", "RegionalProduct",
        "APIDatasetMetaData"]

    api_method = "GetDatasetList"

    querydict = Dict("UserID" => user_id,
                     "Method" => api_method,
                     "ResultFormat" => "JSON")

    r_json = bea_query(API_URL, querydict)
    r_dict = r_json["Dataset"]

    ids = [d["DatasetName"] for d in r_dict if d["DatasetName"] ∉ excluded_sets]
    descriptions = [d["DatasetDescription"] for d in r_dict if d["DatasetName"] ∉ excluded_sets]

    return DataFrame(DatasetName = ids, Description = descriptions)
end

"""
    bea_parameterlist(dataset::String; user_id = USER_ID) -> DataFrame

Return a `DataFrame` of parameter names and attributes for `dataset`.

Example:
```julia
# User ID stored in ~/.beadatarc
julia> bea_parameterlist("NIPA")
5×3 DataFrames.DataFrame. Omitted printing of 1 columns
│ Row │ ParameterName │ ParameterDescription                                           │
│     │ String        │ String                                                         │
├─────┼───────────────┼────────────────────────────────────────────────────────────────┤
│ 1   │ Frequency     │ A - Annual, Q-Quarterly, M-Monthly                             │
│ 2   │ ShowMillions  │ A flag indicating that million-dollar data should be returned. │
│ 3   │ TableID       │ The standard NIPA table identifier                             │
│ 4   │ TableName     │ The new NIPA table identifier                                  │
│ 5   │ Year          │ List of year(s) of data to retrieve (X for All)                │
```
"""
function bea_parameterlist(dataset::String; user_id = USER_ID)

    api_method = "GetParameterList"

    querydict = Dict("UserID" => user_id,
                     "Method" => api_method,
                     "DatasetName" => dataset,
                     "ResultFormat" => "JSON")

    r_json = bea_query(API_URL, querydict)
    r_dict = r_json["Parameter"]

    nparams = length(r_dict)
    parameter_id = Array{String}(undef, nparams)
    parameter_desc = Array{String}(undef, nparams)
    parameter_req = Array{String}(undef, nparams)
    for i = 1:nparams
        parameter_id[i] = r_dict[i]["ParameterName"]
        parameter_desc[i] = r_dict[i]["ParameterDescription"]
        parameter_req[i] = r_dict[i]["ParameterIsRequiredFlag"] == "1" ? "Yes" : "No"
    end

    return DataFrame(ParameterName = parameter_id, ParameterDescription = parameter_desc,
        ParamaterRequired = parameter_req)
end

"""
    bea_parametervalues(dataset::String, param_name:: String;
        user_id = USER_ID) -> DataFrame

Return a `DataFrame` of permissible values for `param_name`, with descriptions.

Example:
```julia
# User ID stored in ~/.beadatarc
julia> bea_parametervalues("NIPA", "TableName");

julia> show(ans[1:3, :])
3×2 DataFrame
│ Row │ Value  │ Description                                                                              │
│     │ String │ String                                                                                │
├─────┼────────┼──────────────────────────────────────────────────────────────────────────────────────────┤
│ 1   │ T10101 │ Table 1.1.1. Percent Change From Preceding Period in Real Gross Domestic Product (A) (Q) │
│ 2   │ T10102 │ Table 1.1.2. Contributions to Percent Change in Real Gross Domestic Product (A) (Q)      │
│ 3   │ T10103 │ Table 1.1.3. Real Gross Domestic Product, Quantity Indexes (A) (Q)                       │

```
"""
function bea_parametervalues(dataset::String, param_name:: String; user_id = USER_ID)

    no_year = ["NIPA", "NIUnderlyingDetail", "FixedAssets", "GDPbyIndustry"]
    no_table_id = ["NIPA", "NIUnderlyingDetail"]
    no_frequency = ["UnderlyingGDPbyIndustry"]

    if  (param_name == "Year" && dataset ∈ no_year) ||
        (param_name == "TableID" && dataset ∈ no_table_id) ||
        (param_name == "Frequency" && dataset ∈ no_frequency)
        error("Parameter $param_name in dataset $dataset not supported.")
    end

    api_method = "GetParameterValues"

    querydict = Dict("UserID" => user_id,
                     "Method" => api_method,
                     "DatasetName" => dataset,
                     "ParameterName" => param_name,
                     "ResultFormat" => "JSON")

    r_json = bea_query(API_URL, querydict)
    r_dict = r_json["ParamValue"]

    nvalues = length(r_dict)
    param_vals = Array{String}(undef, nvalues)
    param_desc = Array{String}(undef, nvalues)

    k = if in(param_name, keys(r_dict[1]))
        param_name
    elseif in("Key", keys(r_dict[1]))
        "Key"
    elseif in("key", keys(r_dict[1]))
        "key"
    else
        param_name * "ID"
    end

    desc = if in("Description", keys(r_dict[1]))
        "Description"
    elseif in("Desc", keys(r_dict[1]))
        "Desc"
    elseif in("Descr", keys(r_dict[1]))
        "Descr"
    else
        "desc"
    end

    for i in 1:nvalues
        param_vals[i] = r_dict[i][k]
        param_desc[i] = r_dict[i][desc]
    end

    return DataFrame(Value = param_vals, Description = param_desc)
end


text:
function download_and_check_md5sum(url::String, filename::String; md5_goal::String = "" )
    need_to_download = true
    if isfile(download_filename)
        if Sys.which("md5sum") == nothing
            println("# System doesn't have md5sum so forcing redownload of " * filename )
        else
            md5_here = ""
            try
                md5sum_output = read(`md5sum $filename`, String)
                md5_here = split(md5sum_output)[1]
                if length(md5_goal) == 0 || md5_here == md5_goal
                    println("# Found ", filename, " so not downloading.")
                    need_to_download = false
                else
                    println("# Found ", filename, " but md5sum doesn't match, so should redownload.")
                end
            catch
                println("# Warning couldn't compute md5sum to verify download of ", filename)
                need_to_download = false
            end
        end
    end
    if need_to_download
        println("# Downloading ", filename)
        download(url, filename)
        if Sys.which("md5sum") == nothing
            @warn "System doesn't have md5sum to check that download of " * filename * " was successful."
        else
            md5sum_output = read(`md5sum $filename`, String)
            md5_here = split(md5sum_output)[1]
            if md5_here != md5_goal
                println("# The md5sum didn't match expectations, make sure the download was successfuly completed.")
            end
        end
    end
    return
end

ext: jl

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# Type definition
mutable struct TaylorSeries{T, S<:Real} <: FunctionSeries{T}
    coeffs::Vector{T}
    offset::S
    function TaylorSeries(coeffs::Vector{T}, offset = 0.) where {T}
        @assert length(coeffs) >= 1
        S = typeof(offset)
        new{T,S}(coeffs, offset)
    end
end

# Basic properties
offset(f::TaylorSeries) = f.offset
domain(f::TaylorSeries) = (-Inf,+Inf)

degree(f::TaylorSeries) = length(f.coeffs) - 1
coefficients(f::TaylorSeries) = f.coeffs

function shift!(f::TaylorSeries{T}, b) where T
    a = offset(f)
    a == b && return f
    for k = 0:degree(f)
        for l = degree(f):-1:k+1
            if T <: AbstractArray
                axpy!((b-a)^(l-k)*binomial(l,k), f[l], f[k])
            else
                f[k] += f[l]*(b-a)^(l-k)*binomial(l,k)
            end
        end
    end
    f.offset = b
    return f
end
shift(f::TaylorSeries, b) = b == offset(f) ? f : shift!(copy(f), b)

function Base.:(==)(f1::TaylorSeries, f2::TaylorSeries)
    offset(f1) == offset(f2) || return (f1 == shift(f2, offset(f1)))
    K = max(degree(f1), degree(f1))
    for k = 0:K
        f1[k] == f2[k] || return false
    end
    return true
end

# Indexing, getting and setting coefficients
Base.eachindex(f::TaylorSeries) = 0:degree(f)
function Base.getindex(f::TaylorSeries, k)
    k < 0 && throw(BoundsError(f, k))
    @inbounds if k > degree(f)
        return zero(f.coeffs[1])
    else
        return f.coeffs[k+1]
    end
end
function Base.setindex!(f::TaylorSeries, v, k)
    k < 0 && return BoundsError(f, k)
    @inbounds if k > degree(f)
        f0 = f[0]
        while degree(f) < k-1
            push!(f.coeffs, zero(f0))
        end
        push!(f.coeffs, v)
        return v
    else
        setindex!(f.coeffs, v, k+1)
        return v
    end
end

# Use as function
function (f::TaylorSeries)(x::Real)
    k = degree(f)
    xa = x - offset(f)
    (xa == zero(xa) || k == 0) && return copy(f[0])
    v = copy(f[k])
    if v isa AbstractArray
        while k > 0
            v = LinearAlgebra.axpby!(one(xa), f[k-1], xa, v)
            k -= 1
        end
    else
        while k > 0
            v = muladd(xa, v, f[k-1])
            k -= 1
        end
    end
    return v
end

# Change number of coefficients
function truncate!(f::TaylorSeries; Kmax::Integer = degree(f), tol::Real = 0, dx::Real = 1.)
    ktol = findlast(k->(norm(f[k])*abs(dx)^k >= tol), eachindex(f))
    Ktol = ktol === nothing ? 0 : ktol-1
    Kmax = min(Kmax, Ktol)
    if Kmax < degree(f)
        resize!(f.coeffs, Kmax + 1)
    end
    return f
end

function setdegree!(f::TaylorSeries, K::Int)
    while degree(f) < K
        push!(f.coeffs, zero(f[0]))
    end
    if degree(f) > K
        resize!(f.coeffs, K+1)
    end
    return f
end

# special purpose constructor
Base.similar(f::TaylorSeries, ::Type{T} = scalartype(f)) where {T} =
    TaylorSeries([zero(T)*f[0]], offset(f))

Base.zero(f::TaylorSeries) = TaylorSeries([zero(f[0])], offset(f))
Base.one(f::TaylorSeries) = TaylorSeries([one(f[0])], offset(f))

# Arithmetic (out of place)
Base.copy(f::TaylorSeries) = TaylorSeries(copy.(coefficients(f)), offset(f))

Base.:-(f::TaylorSeries) = TaylorSeries(.-coefficients(f), offset(f))
Base.:+(f::TaylorSeries) = TaylorSeries(.+coefficients(f), offset(f))

Base.:*(f::TaylorSeries, a::Const) =
    TaylorSeries([c*a for c in coefficients(f)], offset(f))
Base.:*(a::Const, f::TaylorSeries) =
    TaylorSeries([a*c for c in coefficients(f)], offset(f))
Base.:/(f::TaylorSeries, a::Const) =
    TaylorSeries([c/a for c in coefficients(f)], offset(f))
Base.:\(a::Const, f::TaylorSeries) =
    TaylorSeries([a\c for c in coefficients(f)], offset(f))

function Base.:+(f1::TaylorSeries, f2::TaylorSeries)
    offset(f1) == offset(f2) || return (f1 + shift(f2, offset(f1)))
    K = min(degree(f1), degree(f2))
    coeffs = [f1[j]+f2[j] for j = 0:K]
    for j = K+1:degree(f1)
        push!(coeffs, copy(f1[j]))
    end
    for j = K+1:degree(f2)
        push!(coeffs, copy(f2[j]))
    end
    return TaylorSeries(coeffs, offset(f1))
end

function Base.:-(f1::TaylorSeries, f2::TaylorSeries)
    offset(f1) == offset(f2) || return (f1 - shift(f2, offset(f1)))
    K = min(degree(f1), degree(f2))
    coeffs = [f1[j]-f2[j] for j = 0:K]
    for j = K+1:degree(f1)
        push!(coeffs, copy(f1[j]))
    end
    for j = K+1:degree(f2)
        push!(coeffs, -(f2[j]))
    end
    return TaylorSeries(coeffs, offset(f1))
end

Base.:*(f1::TaylorSeries, f2::TaylorSeries) = truncmul(f1, f2)

function truncmul(f1::TaylorSeries, f2::TaylorSeries; kwargs...)
    if offset(f1) != offset(f2)
        _truncmul(f1, shift(f2, offset(f1)); kwargs...)
    else
        _truncmul(f1, f2; kwargs...)
    end
end

function _truncmul(f1::TaylorSeries, f2::TaylorSeries;
                    Kmax = degree(f1) + degree(f2), tol::Real = 0, dx = 1)
    K = min(Kmax, degree(f1) + degree(f2))
    f = TaylorSeries(sizehint!([zero(f1[0])*zero(f2[0])], K+1), offset(f1))
    return truncmul!(f, f1, f2, true, true; Kmax = Kmax, tol = tol, dx = dx)
end

Base.conj(f::TaylorSeries) = TaylorSeries(map(conj, coefficients(f)), offset(f))
Base.adjoint(f::TaylorSeries) = TaylorSeries(map(adjoint, coefficients(f)), offset(f))
Base.transpose(f::TaylorSeries) = TaylorSeries(map(transpose, coefficients(f)), offset(f))

LinearAlgebra.tr(f::TaylorSeries) = TaylorSeries(map(tr, coefficients(f)), offset(f))

Base.real(f::TaylorSeries) = TaylorSeries(map(real, coefficients(f)), offset(f))
Base.imag(f::TaylorSeries) = TaylorSeries(map(imag, coefficients(f)), offset(f))

# Arithmetic (in place / mutating methods)
function Base.copy!(fdst::TaylorSeries, fsrc::TaylorSeries)
    fdst.offset = offset(fsrc)
    setdegree!(fdst, degree(fsrc))
    for j = 0:degree(fsrc)
        fdst[j] = copy(fsrc[j])
    end
    return fdst
end
function LinearAlgebra.rmul!(f::TaylorSeries, α::Number)
    if eltype(f) <: Number
        rmul!(coefficients(f), α)
    else
        for k in eachindex(f)
            rmul!(f[k], α)
        end
    end
    return f
end
function LinearAlgebra.lmul!(α::Number, f::TaylorSeries)
    if eltype(f) <: Number
        lmul!(α, coefficients(f))
    else
        for k in eachindex(f)
            lmul!(α, f[k])
        end
    end
    return f
end
LinearAlgebra.axpy!(α::Number, fx::TaylorSeries, fy::TaylorSeries) =
    truncadd!(fy, fx, α)
LinearAlgebra.axpby!(α::Number, fx::TaylorSeries, β::Number, fy::TaylorSeries) =
    truncadd!(fy, fx, α, β)
LinearAlgebra.mul!(fy::TaylorSeries, s::Number, fx::TaylorSeries, α = true, β = false) =
    truncmul!(fy, s, fx, α, β)
LinearAlgebra.mul!(fy::TaylorSeries, fx::TaylorSeries, s::Number, α = true, β = false) =
    truncmul!(fy, fx, s, α, β)
LinearAlgebra.mul!(f::TaylorSeries, f1::TaylorSeries, f2::TaylorSeries,
                    α = true, β = false) = truncmul!(f, f1, f2, α, β)

function truncadd!(fy::TaylorSeries, fx::TaylorSeries, α = true, β = true;
                    Kmax::Integer = max(iszero(β) ? 0 : degree(fy), degree(fx)),
                    tol::Real = 0, dx = 1)
    offset(fx) == offset(fy) || throw(DomainMismatch())
    Kx = degree(fx)
    Ky = min(Kmax, iszero(β) ? Kx : max(Kx, degree(fy)))
    setdegree!(fy, Ky)
    # if offset(fx) == offset(fy)
        if eltype(fy) <: Number
            if Kx > Ky
                axpby!(α, view(coefficients(fx, 1:Ky+1)), β, coefficients(fy))
            else
                axpby!(α, coefficients(fx), β, view(coefficients(fy), 1:Kx+1))
                lmul!(β, view(coefficients(fy), (Kx+2):(Ky+1)))
            end
        else
            for k in 0:Kx
                axpby!(α, fx[k], β, fy[k])
            end
            for k in Kx+1:Ky
                lmul!(β, fy[k])
            end
        end
    # else
    #     dx = offset(fy) - offset(fx)
    #     lmul!(β, fy)
    #     for k = 0:Kx
    #         for l = 0:min(k, Ky)
    #             if eltype(fy) <: Number
    #                 fy[l] += α*fx[k]*dx^(l-k)*binomial(l,k)
    #             else
    #                 axpy!(α*dx^(l-k)*binomial(l,k), fx[k], fy[l])
    #             end
    #         end
    #     end
    # end
    return truncate!(fy; tol = tol, dx = dx)
end

truncmul!(fdst::TaylorSeries, α₁::Number, fsrc::TaylorSeries, α₂ = true, β = false;
            kwargs...) = truncadd!(fdst, fsrc, α₁ * α₂, β; kwargs...)
truncmul!(fdst::TaylorSeries, fsrc::TaylorSeries, α₁::Number, α₂ = true, β = false;
            kwargs...) = truncadd!(fdst, fsrc, α₁ * α₂, β; kwargs...)

function truncmul!(F::TaylorSeries, F1::TaylorSeries, F2::TaylorSeries,
                    α = true, β = false;
                    Kmax::Integer = max(iszero(β) ? 0 : degree(F), degree(F1)+degree(F2)),
                    tol::Real = 0,
                    dx = 1)
    offset(F) == offset(F1) == offset(F2) || throw(DomainMismatch())
    K1 = degree(F1)
    K2 = degree(F2)
    K = min(Kmax, max(iszero(β) ? 0 : degree(F), K1+K2))
    setdegree!(F, K)
    Threads.@threads for k = 0:K
        fk = F[k]
        if fk isa AbstractArray
            T = eltype(fk)
            if iszero(β)
                fill!(fk, T(β))
            elseif !isone(β)
                rmul!(fk, T(β))
            end
            for k1 = max(0, k-K2):min(K1, k)
                k2 = k - k1
                fk1 = F1[k1]
                fk2 = F2[k2]
                if fk1 isa AbstractArray && fk2 isa AbstractArray
                    mul!(fk, fk1, fk2, T(α), T(true))
                elseif fk1 isa AbstractArray
                    axpy!(T(fk2*α), fk1, fk)
                elseif fk2 isa AbstractArray
                    axpy!(T(fk1*α), fk2, fk)
                else
                    @warn "unexpected branch"
                    fk .= fk + (fk1*fk2*α)
                end
            end
        else
            if iszero(β)
                fk = zero(fk)
            elseif !isone(β)
                fk *= β
            end
            for k1 = max(0, k-K2):min(K1, k)
                k2 = k - k1
                fk1 = F1[k1]
                fk2 = F2[k2]
                fk += fk1*fk2*α
            end
            F[k] = fk
        end
    end
    if tol != 0
        truncate!(F; tol = tol, dx = dx)
    end
    return F
end

# Inner product and norm
function localdot(f1::TaylorSeries, f2::TaylorSeries)
    offset(f1) == offset(f2) || localdot(f1, shift(f2, offset(f1)))
    K1 = degree(f1)
    K2 = degree(f2)
    coeffs = let K = K1 + K2
        [sum(dot(f1[k1], f2[k-k1]) for k1 = max(0, k-K2):min(k, K1)) for k=0:K]
    end
    return TaylorSeries(coeffs, offset(f1))
end

function differentiate(f::TaylorSeries)
    if degree(f) >= 1
        return TaylorSeries([k*f[k] for k = 1:degree(f)], offset(f))
    else
        return zero(f)
    end
end
function integrate(f::TaylorSeries, (a,b)::Tuple{Real,Real})
    s = f[0]*(b-a)
    c = offset(f)
    db = b - c
    da = a - c
    for k = 1:degree(f)
        s += f[k]*(db^(k+1) - da^(k+1))/(k+1)
    end
    return s
end

# Simple Taylor series least squares fitting:
# typically not needed for large number of points,
# but should work with matrix valued functions etc
# => generality over efficiency
function fit(f, ::Type{TaylorSeries}, domain::Tuple{Real,Real};
                Kmax = 10, numpoints = Kmax+1, offset = (domain[1]+domain[2])/2)

    a, b = domain
    x = range(a, stop = b, length = numpoints)
    dx = x .- offset
    coeffs = similar(x, (numpoints, Kmax+1))
    coeffs[:, 1] .= 1
    for k = 1:Kmax
        coeffs[:, k+1] .= view(coeffs, :, k) .* dx
    end
    TaylorSeries(pinv(coeffs) * f.(x), offset)
end

# Inverse (and square root)
# Warning: experimental
function Base.inv(t::TaylorSeries; tol = 1e-10, Kmax = 10*degree(t))
    t0inv = inv(t[0])
    coeffs = [t0inv]
    K = degree(t)
    for k = 1:Kmax
        invtk = -t0inv*sum(coeffs[i+1]*t[k-i] for i = max(0,k-K):(k-1))
        if norm(invtk) < tol && norm(coeffs[end]) < tol
            break
        end
        push!(coeffs, invtk)
    end
    return TaylorSeries(coeffs, offset(t))
end


text:
using DiffEqSensitivity, OrdinaryDiffEq, LinearAlgebra, SparseArrays, Zygote, LinearSolve
using AlgebraicMultigrid: AlgebraicMultigrid
using Test

foop(u, p, t) = jac(u, p, t) * u
jac(u, p, t) = spdiagm(0=>p)
paramjac(u, p, t) = SparseArrays.spdiagm(0=>u)
@Zygote.adjoint foop(u, p, t) = foop(u, p, t), delta->(jac(u, p, t)' * delta, paramjac(u, p, t)' * delta, zeros(length(u)))

n = 2
p = collect(1.0:n)
u0 = ones(n)
tspan = [0.0, 1]
odef = ODEFunction(foop; jac=jac, jac_prototype=jac(u0, p, 0.0), paramjac=paramjac)
function g_helper(p; alg=Rosenbrock23(linsolve=LUFactorization()))
    prob = ODEProblem(odef, u0, tspan, p)
    soln = Array(solve(prob, alg; u0=prob.u0, p=prob.p, abstol=1e-4, reltol=1e-4, sensealg=InterpolatingAdjoint()))[:, end]
    return soln
end
function g(p; kwargs...)
    soln = g_helper(p; kwargs...)
    return sum(soln)
end
@test isapprox(exp.(p), g_helper(p); atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), Zygote.gradient(g, p)[1]; atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), g_helper(p; alg=Rosenbrock23(linsolve=KLUFactorization())); atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), Zygote.gradient(p->g(p; alg=Rosenbrock23(linsolve=KLUFactorization())), p)[1]; atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=LUFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), Zygote.gradient(p->g(p; alg=ImplicitEuler(linsolve=LUFactorization())), p)[1]; atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), Zygote.gradient(p->g(p; alg=ImplicitEuler(linsolve=UMFPACKFactorization())), p)[1]; atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), Zygote.gradient(p->g(p; alg=ImplicitEuler(linsolve=KrylovJL_GMRES())), p)[1]; atol=1e-1, rtol=1e-1)


text:
@testset "rsoc problems" begin
    @testset "rsoc1_test" begin
        #=
        primal
            min 0a + 0b - 1x - 1y
        s.t.
            a    == 1/2 :w_2
            b    == 1   :w_1
            2a*b >= x^2+y^2  :w_3, w_4, w_5, w_6
        dual
            max w_2 + (1/2)w_1
        s.t.
            w_5 == -1
            w_6 == -1
            w_1 + w_3 == 0
            w_2 + w_4 == 0
            (w_3, w_4, w_5, w_6) \in RotatedSecondOrderCone
        =#
        primal_model = rsoc1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 2
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) ==
              Set([(VAF{Float64}, MOI.RotatedSecondOrderCone)],)
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VAF{Float64},MOI.RotatedSecondOrderCone}(),
        ) == 1
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        @test MOI.constant(obj) == 0.0
        @test Set(MOI.coefficient.(obj.terms)) == Set([1.0; 0.5])
    end

    @testset "rsoc2_test" begin
        #=
        primal
            min 0a + 0b - 1x - 1y
        s.t.
            a    == 1/2
            b    == 1
            2a*b >= x^2+y^2
        dual
            max w_2 + (1/2)w_1
        s.t.
            w_5 == -1
            w_6 == -1
            w_1 + w_3 == 0
            w_2 + w_4 == 0
            (w_3, w_4, w_5, w_6) \in RotatedSecondOrderCone
        =#
        primal_model = rsoc2_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 4
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraintTypesPresent())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VVF, MOI.RotatedSecondOrderCone)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.RotatedSecondOrderCone}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 2
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.terms) == [-1.0; -0.5]

        eq_con1, eq_con2 = MOI.get(
            dual_model,
            MOI.ListOfConstraintIndices{SAF{Float64},MOI.EqualTo{Float64}}(),
        )

        eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), eq_con1)
        eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), eq_con1)
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == -1.0
        eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), eq_con2)
        eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), eq_con2)
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == -1.0

        rsoc_con = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{VVF,MOI.RotatedSecondOrderCone}(1),
        )
        @test rsoc_con.variables == VI.(1:4)

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{VAF{Float64},MOI.RotatedSecondOrderCone}(
            1,
        )] == VI.(1:4)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] == eq_con1
        @test primal_var_dual_con[VI(2)] == eq_con2
    end
end


text:
# Transformations on dictionaries of xy-coordinates

using LinearAlgebra, Statistics
export transform, scale, rotate, translate, recenter

"""
`_transform(xy,A,b)`  appliex the affine transformation
`v |--> Av+b` to all entries in `xy`.
"""
function _transform(xy::Dict, A::Array{S,2}, b::Vector{T}) where {S<:Real,T<:Real}

    # apply the transformation x |--> Ax+b to all coordinates
    for k in keys(xy)
        xy[k] = A * xy[k] + b
    end
    nothing
end

"""
`_translate(xy,b)` adds the vector `b` to all vectors in `xy`.
"""
function _translate(xy::Dict{T,Vector{S}}, b::Vector{R}) where {T,S,R}
    A = Matrix{Float64}(I, 2, 2)
    _transform(xy, A, b)
end

"""
`_scale(xy,m)` multiplies all the vectors in `xy` by the scalar `m`.
"""
function _scale(xy::Dict{T,Vector{S}}, m) where {T,S}
    A = m * Matrix{Float64}(I, 2, 2)
    _transform(xy, A, [0, 0])
end

"""
`_rotate(xy,theta)` rotates all vectors in `xy`
by `theta`.
"""
function _rotate(xy::Dict{T,Vector{S}}, theta) where {T,S}
    A = zeros(2, 2)
    A[1, 1] = cos(theta)
    A[1, 2] = -sin(theta)
    A[2, 1] = sin(theta)
    A[2, 2] = cos(theta)
    b = [0, 0]
    _transform(xy, A, b)
end

"""
`_recenter(xy)` makes translates the points in `xy` so their
average is `[0,0]`.
"""
function _recenter(xy::Dict{T,Vector{S}}) where {T,S}
    n = length(xy)
    if n == 0
        return
    end
    b = [0.0, 0.0]
    for v in keys(xy)
        b += xy[v]
    end
    b /= n
    _translate(xy, -b)
end

"""
`recenter(G)` translates the graph's drawing so that the center of mass
of the vertices is at the origin.
"""
function recenter(G::SimpleGraph)
    ensure_embed(G)
    _recenter(G.cache[:xy])
end

"""
`transform(G,A,b)` applies an affine transformation to all coordinates
in the graph's drawing. Here `A` is 2-by-2 matrix and `b` is a 2-vector.
Each point `p` is mapped to `A*p+b`.
"""
function transform(G::SimpleGraph, A::Array{S,2}, b::Vector{T}) where {S<:Real,T<:Real}
    ensure_embed(G)
    _transform(G.cache[:xy], A, b)
end


text:
using FillArrays

"""
    ADI.SDIAlgorithm <: ADI.ADIAlgorithm

Spectral differential imaging (SDI) algorithms. These work on 4-D SDI tensors. To use these algorithms, simply treat them like functions (or call [`process`](@ref))

```julia
(::SDIAlgorithm)(data::AbstractArray{T,4}, angles, scales; [ref] kwargs...)
```

The data is expected to be laid out in `(nλ, nf, ny, nx)` format, so you may need to `permutedims` before processing the data. The `scales` correspond to the relative wavelength scales for each spectrum, and can be retrieved with `HCIToolbox.scale_list`.

# Algorithms

The current SDI implementations are
* [`SingleSDI`](@ref)
* [`DoubleSDI`](@ref)
* [`SliceSDI`](@ref)
"""
abstract type SDIAlgorithm <: ADIAlgorithm end

"""
    SingleSDI(alg)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in a single pass. This means that each channel will be scaled and then concatenated together, so an SDI tensor `(nλ, nf, y, x)` becomes a stack `(nλ * nf, y, x)` which is processed by the underlying `alg` like ADI data.

!!! tip
    `SingleSDI` is the default SDI mode. This means instead of writing
    ```julia
    SingleSDI(PCA(15))(data, angles, scales)
    ```
    you can just write
    ```julia
    PCA(15)(data, angles, scales)
    ```
"""
struct SingleSDI{ALG<:ADIAlgorithm} <: SDIAlgorithm
    alg::ALG
end

process(alg, spcube::AbstractArray{T,4}, angles, scales; kwargs...) where {T} =
    process(SingleSDI(alg), spcube, angles, scales; kwargs...)

function process(sdi::SingleSDI, spcube::AbstractArray{T,4}, angles, scales; method=:deweight, kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    big_cube = scale_and_stack(spcube, scales)
    angs = repeat(angles, inner=nλ)
    # do single-pass reconstruction
    if :ref in keys(kwargs)
        big_cube_ref = scale_and_stack(kwargs[:ref], scales)
        big_resid_cube = subtract(sdi.alg, big_cube; angles=angs, kwargs..., ref=big_cube_ref)
    else
        big_resid_cube = subtract(sdi.alg, big_cube; angles=angs, kwargs...)
    end

    # bin across spectral dim
    resid_cube = invscale_and_collapse(big_resid_cube, scales, frame_size)
    # derotate and combine
    return collapse!(resid_cube, angles; method=method, kwargs...)
end

"""
    DoubleSDI(alg)
    DoubleSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses `alg_spec` to reduce each spectral cube slice in the SDI tensor. Then, the spectral residual frames will be reduced using `alg_temp`, which will include the derotation and final combination.

The difference between [`DoubleSDI`](@ref) and [`SliceSDI`](@ref) is that `DoubleSDI` does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. `SliceSDI` does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

# Examples

```julia
julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = DoubleSDI(Classic(), GreeDS(15))(data, angles, scales)
```
"""
struct DoubleSDI{ALG<:ADIAlgorithm, ALG2<:ADIAlgorithm} <: SDIAlgorithm
    alg_spec::ALG
    alg_temp::ALG2
end

DoubleSDI(alg) = DoubleSDI(alg, alg)

function process(sdi::DoubleSDI, spcube::AbstractArray{T,4}, angles, scales; method=:deweight, kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    spec_resids = similar(spcube, n, ny, nx)
    # do first pass in spectral domain
    Threads.@threads for n in axes(spcube, 2)
        cube = @view spcube[:, n, :, :]
        scaled_cube = scale(cube, scales)
        angs = Zeros(nλ)
        if :ref in keys(kwargs)
            cube_ref = view(kwargs[:ref], :, n, :, :)
            scaled_cube_ref = scale(cube_ref, scales)
            R = subtract(sdi.alg_spec, scaled_cube; angles=angs, kwargs..., ref=scaled_cube_ref)
        else
            R = subtract(sdi.alg_spec, scaled_cube; angles=angs, kwargs...)
        end
        spec_resid = invscale(R, scales, frame_size)
        spec_resids[n, :, :] .= collapse(spec_resid)
    end
    # do second pass in temporal domain
    return sdi.alg_temp(spec_resids, angles; method=method, angles=angles, kwargs...)
end

"""
    SliceSDI(alg)
    SliceSDI(alg_spec, alg_temp)

A wrapper algorithm for spectral differential imaging (SDI) data reduced in two passes. The first pass uses `alg_temp` to reduce each temporal cube slice in the SDI tensor. These residuals will be rescaled and stacked into a new cube. Then, the temporal residual frames will be reduced using `alg_spec`, which will include the derotation and final combination.

The difference between [`SliceSDI`](@ref) and [`DoubleSDI`](@ref) is that `DoubleSDI` does its first pass in the spectral slice, effectively collapsing the slice before performing ADI on the residual cube. `SliceSDI` does its first pass in the temporal slice, collapsing it first before performing ADI on the residual cube.

# Examples

```julia
julia> data, angles, scales = # load data...

# Median subtraction for each spectral slice,
# GreeDS{PCA} subtraction on spectral residual cube
julia> res = SliceSDI(Classic(), GreeDS(15))(data, angles, scales)
```
"""
struct SliceSDI{ALG<:ADIAlgorithm, ALG2<:ADIAlgorithm} <: SDIAlgorithm
    alg_spec::ALG
    alg_temp::ALG2
end

SliceSDI(alg) = SliceSDI(alg, alg)

function process(sdi::SliceSDI, spcube::AbstractArray{T,4}, angles, scales; kwargs...) where T
    nλ, n, ny, nx = size(spcube)
    frame_size = (ny, nx)
    temp_resids = similar(spcube, nλ, ny, nx)
    # do first pass in temporal domain
    Threads.@threads for n in axes(spcube, 1)
        cube = @view spcube[n, :, :, :]
        if :ref in keys(kwargs)
            cube_ref = @view kwargs[:ref][n, :, :, :]
            temp_resids[n, :, :] .= sdi.alg_temp(cube, angles; kwargs..., ref=cube_ref)
        else
            temp_resids[n, :, :] .= sdi.alg_temp(cube, angles; kwargs...)
        end
    end
    # do second pass in temporal domain
    scaled_resid_cube = scale(temp_resids, scales)
    angs = Zeros(nλ)
    resid = sdi.alg_spec(scaled_resid_cube, angs; angles=angs, kwargs...)
    return invscale(resid, maximum(scales))
end


text:
# This file is a part of InverseFunctions.jl, licensed under the MIT License (MIT).

import Test
import InverseFunctions
import Documenter

Test.@testset "Package InverseFunctions" begin
    include("test_functions.jl")
    include("test_inverse.jl")

    # doctests
    Documenter.DocMeta.setdocmeta!(
        InverseFunctions,
        :DocTestSetup,
        :(using InverseFunctions);
        recursive=true,
    )
    Documenter.doctest(InverseFunctions)
end # testset


text:
function reveal_loaded_packages(mod)
    for Mod in Base.loaded_modules_array()
        if !Core.isdefined(mod, nameof(Mod))
            Core.eval(mod, Expr(:const, Expr(:(=), nameof(Mod), Mod)))
        end
    end
end

function import_packages!(packages,mod,deffered = Set{Symbol}())
    for p in packages
        Fezzik.reveal_loaded_packages(mod)
        if isdefined(mod,p)
            if typeof(Core.eval(mod,p)) === Module
                Core.eval(mod,p)
                println("[$p] already loaded")
                delete!(packages,p)
            else
                println("[$p] is not a Module")
            end
        else
            !(p in deffered) && println("trying to import $p")
            try
                Core.eval(mod, :(import $p))
                newmod = Core.eval(mod,p);
                isdefined(newmod,:__init__) && Core.eval(newmod,:__init__)()
                Fezzik.reveal_loaded_packages(mod)
                delete!(packages,p)
                println("[$p] loaded")
            catch e
                !(p in deffered) && println("failed to import $p deffering")
            end
        end
    end
    deepcopy(packages)
end

function brute_import_packages!(packages,mod,n = 5)
    deffered = Set{Symbol}()
    for i=1:n
        deffered = Fezzik.import_packages!(packages,mod,deffered)
    end
    for p in packages
        try
            occursin("_jll","$p") && continue
            Pkg.add("$p")
            Core.eval(mod, :(import $p))
            delete!(packages,p)
            for i=1:3
                deffered = Fezzik.import_packages!(packages,mod,deffered)
            end
        catch e
            @warn e
            @warn "could not import $p"
        end
    end
    !isempty(packages) && (@warn "failed to import packages" packages)
    nothing
end


text:
using DemoPackageXL
using Documenter

DocMeta.setdocmeta!(DemoPackageXL, :DocTestSetup, :(using DemoPackageXL); recursive=true)

makedocs(;
    modules=[DemoPackageXL],
    authors="Xubo Leng <hsupo.leng@gmail.com> and contributors",
    repo="https://github.com/HsupoLeng/DemoPackageXL.jl/blob/{commit}{path}#{line}",
    sitename="DemoPackageXL.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://HsupoLeng.github.io/DemoPackageXL.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/HsupoLeng/DemoPackageXL.jl",
    devbranch="main",
)


text:
import Mustache
import Distributions

type Report
   num::Integer
   name::String
   times::Vector{Float64}
   analysis::SampleAnalysis
   outliers::Outliers
end

function render_iter(s::String, iter)
    [{s => x} for x in iter]
end

function report_to_dict(report::Report)
   kde_times = Distributions.kde(report.times)
   {"name" => report.name,
    "number" => report.num,
    "times" => render_iter("x", report.times),
    "kdetimes" => render_iter("x", kde_times.x),
    "kdepdf" => render_iter("x", kde_times.density),
    "mean_point" => report.analysis.mean.point,
    "mean_confidence_level" => report.analysis.mean.confidence_level,
    "mean_lbound" => report.analysis.mean.lbound,
    "mean_ubound" => report.analysis.mean.ubound,
    "outlier_effect" => report.analysis.outlier_variance.effect,
    "outlier_fraction" => report.analysis.outlier_variance.frac,
    "std_lbound" => report.analysis.std.lbound,
    "std_point" => report.analysis.std.point,
    "std_ubound" => report.analysis.std.ubound}
end

function read_css(template_dir::String)
    readall(abspath(joinpath(template_dir, "criterion.css")))
end

function report(rfile::String, reports::Vector{Report})
    report_dicts = {report_to_dict(r) for r in reports}
    tmpl = Mustache.template_from_file("../templates/report.tpl")
    template_dir = abspath("../templates/")
    fh = open(rfile,"w")
    Mustache.render(fh, tmpl, {"name" => "aaa",
                    "include" => [{"dir" => string("file://", template_dir),
				   "css" => read_css(template_dir)}],
		    "report" => report_dicts})
    close(fh)
end

function report(rfile::String, r::Report)
    report(rfile, [r])
end

function report(rs::Vector{Report})
    report("criterion.html", rs)
end

function report(r::Report)
    report("criterion.html", [r])
end

# JUnit Output
function junit(reports::Vector{Report})

    attr_esc = (esc) -> begin
        if     esc == '\'' "'"
        elseif esc == '"'  """
        elseif esc == '<'  "<"
        elseif esc == '>'  ">"
        elseif esc == '&'  "&"
        else string(esc)
        end
    end

    esc_html = (str) -> join([attr_esc(c) for c in str], "")

    report_out = String[]
    ntests = length(reports)
    header = string("<?xml versions\"1.0\" encoding=\"UTF-8\"?>\n",
                    "<testsuite name=\"Criterion benchmarks\" tests=\"$ntests\">\n")
    push!(report_out, header)
    for report in reports
       name = esc_html(report.name)
       time = report.analysis.mean.point
       push!(report_out, "<testcase name=\"$name\" time=\"$time\"/>\n")
   end
   push!(report_out, "</testsuite>\n")
   join(report_out, "")
end

function junit(report::Report)
    junit([report])
end


text:
using Dispersal, Test, DimensionalData, Setfield, OffsetArrays, DynamicGrids, Random
using DynamicGrids: WritableGridData, setneighbor!, SimData, Extent, grids, neighborhood,
      mapsetneighbor!
using Dispersal: dispersalprob

struct TestFormulation <: KernelFormulation end
(f::TestFormulation)(d) = 1.0

# Dispersal in radius 2 neighborhood
hood = DispersalKernel{2}(; formulation=TestFormulation())

@testset "ConstructionBase compat" begin
    hood2 = @set hood.distancemethod = AreaToArea(10.0)
    @test hood2.distancemethod == AreaToArea(10.0)
    @test DynamicGrids.radius(hood2) == 2
    @test Dispersal.formulation(hood2) == TestFormulation()
    @test hood2.cellsize == 1.0
    @test hood2.kernel == hood.kernel
    @test hood2.kernel !== hood.kernel
    @test all(hood2.kernel .== 0.04)
end

@testset "DistanceMethod" begin
    cellsize = 1.0
    f = identity
    x, y = 10, 10

    @test dispersalprob(f, CentroidToCentroid(), x, y, cellsize) == 10sqrt(2)

    subsampling = (1.0, 2.0, 10.0, 20.0, 50.0, 70.0, 80.0)

    proba2c = map(subsampling) do ss
        dispersalprob(f, AreaToCentroid(ss), x, y, cellsize)
    end
    @test proba2c[1] == 10sqrt(2)
    @test issorted(proba2c)
    @test proba2c[end] ≈ proba2c[end-1] atol=1e-5

    proba2a = map(subsampling) do ss
        dispersalprob(f, AreaToArea(ss), x, y, cellsize)
    end
    @test proba2a[1] == 10sqrt(2)
    @test issorted(proba2a)
    @test proba2a[end] ≈ proba2a[end-1] atol=1e-5
end

@testset "buildkernel" begin
    #buildkernel(formulation, distancemethod, cellsize, r)
end

@testset "setneighbour!" begin
    hood_index = (2, 2)
    dest_index = (2, 2)
    rule = OutwardsBinaryDispersal(;
        neighborhood=hood,
        prob_threshold=0.0,
    )
    ruleset = Ruleset(rule)
    Dispersal.kernel(hood)

    state = true
    init = zeros(Bool, 5, 5)
    simdata = SimData(Extent(init=(_default_=init,), tspan=1:1), ruleset)
    griddata = WritableGridData(first(grids(simdata)))
    @test setneighbor!(griddata, neighborhood(rule), rule, state, hood_index, dest_index) == 1
    @test DynamicGrids.dest(griddata) == OffsetArray(
        [0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 1 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0], -1:7, -1:7)
    @test mapsetneighbor!(griddata, neighborhood(rule), rule, state, (5, 5)) == 24
    @test DynamicGrids.dest(griddata) == OffsetArray(
        [0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0
         0 0 0 1 0 0 0 0 0
         0 0 0 0 1 1 1 1 1
         0 0 0 0 1 1 1 1 1
         0 0 0 0 1 1 0 1 1
         0 0 0 0 1 1 1 1 1
         0 0 0 0 1 1 1 1 1], -1:7, -1:7)
end

init =  [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 81. 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

test2 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 9.0 9.0 9.0 0.0 0.0;
         0.0 0.0 0.0 9.0 9.0 9.0 0.0 0.0;
         0.0 0.0 0.0 9.0 9.0 9.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

test3 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 1.0 2.0 3.0 2.0 1.0 0.0;
         0.0 0.0 2.0 4.0 6.0 4.0 2.0 0.0;
         0.0 0.0 3.0 6.0 9.0 6.0 3.0 0.0;
         0.0 0.0 2.0 4.0 6.0 4.0 2.0 0.0;
         0.0 0.0 1.0 2.0 3.0 2.0 1.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
         0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]

@testset "inwards" begin
    hood = DispersalKernel{1}(; formulation=TestFormulation())
    ruleset = Ruleset(InwardsPopulationDispersal(;neighborhood=hood))
    output = ArrayOutput(init; tspan=1:3)
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] == test2
    @test output[3] == test3
end

@testset "outwards" begin
    hood = DispersalKernel{1}(; formulation=TestFormulation())
    ruleset = Ruleset(OutwardsPopulationDispersal(;neighborhood=hood))
    output = ArrayOutput(init; tspan=1:3)
    sim!(output, ruleset)
    @test output[1] == init
    @test output[2] == test2
    @test output[3] == test3
end

@testset "switched inwards" begin
    layerdata = cat([1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0  1.0  1.0  1.0],
                    [0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0
                     0.0  0.0  0.0  0.0  1.0  1.0  1.0],
                    [1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0
                     1.0  1.0  1.0  1.0  0.0  0.0  0.0]; dims=3)

    layer = DimensionalArray(layerdata, (X(1:7), Y(1:7), Ti(1:3)))

    init =  [0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 9.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 9.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0]

    test2 = [0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  9.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  0.0  0.0]

    test3 = [0.0  0.0  0.0  0.0  0.0  0.0  0.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             0.0  0.0  0.0  0.0  0.0  1.0  1.0;
             1.0  1.0  1.0  0.0  0.0  0.0  0.0;
             1.0  1.0  1.0  0.0  0.0  0.0  0.0;
             1.0  1.0  1.0  0.0  0.0  0.0  0.0]

    threshold = 0.5
    hood = DispersalKernel{1}(; formulation=TestFormulation())
    switched = SwitchedInwardsPopulationDispersal(; neighborhood=hood, threshold=threshold, layerkey=:layerdata)
    ruleset = Ruleset(switched)

    output = ArrayOutput(init; tspan=1:4)
    sim!(output, ruleset; aux=(layerdata=layerdata,))

    @test output[1] == init
    @test output[2] == test2
    @test output[3] == test3
end


text:
@testset "Structured Matrices" begin
    @testset "Diagonal" begin
        N = 3
        rrule_test(Diagonal, randn(N, N), (randn(N), randn(N)))
        D = Diagonal(randn( N))
        rrule_test(Diagonal, D, (randn(N), randn(N)))
        # Concrete type instead of UnionAll
        rrule_test(typeof(D), D, (randn(N), randn(N)))

        # TODO: replace this with a `rrule_test` once we have that working
        # see https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/24
        res, pb = rrule(Diagonal, [1, 4])
        @test pb(10*res) == (NO_FIELDS, [10, 40])
        comp = Composite{typeof(res)}(; diag=10*res.diag)  # this is the structure of Diagonal
        @test pb(comp) == (NO_FIELDS, [10, 40])
    end

    @testset "::Diagonal * ::AbstractVector" begin
        N = 3
        rrule_test(
            *,
            randn(N),
            (Diagonal(randn(N)), Diagonal(randn(N))),
            (randn(N), randn(N)),
        )
    end
    @testset "diag" begin
        N = 7
        rrule_test(diag, randn(N), (randn(N, N), randn(N, N)))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), randn(N, N)))
        rrule_test(diag, randn(N), (randn(N, N), Diagonal(randn(N))))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), Diagonal(randn(N))))
    end
    @testset "Symmetric" begin
        N = 3
        rrule_test(Symmetric, randn(N, N), (randn(N, N), randn(N, N)))
    end
    @testset "$f" for f in (Adjoint, adjoint, Transpose, transpose)
        n = 5
        m = 3
        rrule_test(f, randn(m, n), (randn(n, m), randn(n, m)))
        rrule_test(f, randn(1, n), (randn(n), randn(n)))
    end
    @testset "$T" for T in (UpperTriangular, LowerTriangular)
        n = 5
        rrule_test(T, T(randn(n, n)), (randn(n, n), randn(n, n)))
    end
end


text:
using Elly
using Base.Test

function test_hdfs(host="localhost", port=9000)
    hdfsclnt = HDFSClient(host, port)

    exists(hdfsclnt, "/tmp") || mkdir(hdfsclnt, "/tmp")

    println("listing files in root folder...")
    dirtree = readdir(hdfsclnt, "/")
    println(dirtree)
    @test "tmp" in dirtree

    println("getting server defaults...")
    defs = hdfs_server_defaults(hdfsclnt)
    for k in keys(defs)
        println("\t$k => $(defs[k])")
    end
    @test hdfs_default_block_size(hdfsclnt) > 0
    @test hdfs_default_replication(hdfsclnt) > 0

    println("getting file system status...")
    fs_status = hdfs_status(hdfsclnt)
    for k in keys(fs_status)
        println("\t$k => $(fs_status[k])")
    end
    @test hdfs_capacity(hdfsclnt) > 0

    println("/tmp should be a directory")
    @test isdir(hdfsclnt, "/tmp")

    println("du should be >= 0")
    @test du(hdfsclnt, "/") >= 0

    println("stat /tmp")
    st = stat(hdfsclnt, "/tmp")
    println(st)
    @test st.name == "/tmp"

    println("create a temporary dir...")
    cd(hdfsclnt, "/tmp")
    foo_dir = HDFSFile(hdfsclnt, "foo")
    mkdir(foo_dir)

    println("create a temporary file...")
    cd(hdfsclnt, "foo")
    bar_file = HDFSFile(hdfsclnt, "bar")
    teststr = "hello world\n"
    nloops = 1000
    tic()
    open(bar_file, "w") do f
        for idx in 1:nloops
            write(f, teststr)
        end
    end
    println("...done in $(toc()) secs")
    println("verify file size to be $(length(teststr)*nloops)...")
    @test filesize(bar_file) == length(teststr) * nloops

    println("delete file...")
    rm(bar_file)

    println("create a large file...")
    size_bytes = 128 * 1000 * 1000
    nloops = 5
    A = rand(UInt8, size_bytes)
    open(bar_file, "w") do f
        for idx in 1:nloops
            tic()
            write(f, A)
            println("...block written in $(toc()) secs")
        end
    end
    println("verify file size to be $(sizeof(A) * nloops)...")
    @test filesize(bar_file) == sizeof(A) * nloops

    println("read and verify...")
    B = Array(UInt8, size_bytes)
    open(bar_file, "r") do f
        for idx in 1:nloops
            tic()
            read!(f, B)
            println("...block read in $(toc()) secs")
            @test A == B
        end
    end

    println("blocks for $bar_file")
    blocks = hdfs_blocks(bar_file)
    println(blocks)

    println("setting replication factor for $bar_file")
    @test hdfs_set_replication(bar_file, 2)

    cd(hdfsclnt, "/tmp/foo")
    NFILES = 1000
    println("create many ($NFILES) files...")
    for idx in 1:NFILES
        bar_file = HDFSFile(hdfsclnt, "bar$idx")
        open(bar_file, "w") do f
            for idx in 1:nloops
                write(f, teststr)
            end
        end
    end
    allfiles = readdir(hdfsclnt, "/tmp/foo")
    println("delete many ($NFILES) files...")
    for idx in 1:NFILES
        bar_file = HDFSFile(hdfsclnt, "bar$idx")
        rm(bar_file)
    end
    @test length(allfiles) >= NFILES
end


text:
import Base: length, step, mean, rand, var

# model.jl implements AbstractModel and SDEModel.
include("model.jl")
# samplepath.jl implements AbstractSamplePath and SamplePath and
# extends Base.step, Base.length, to the appropriate parametric
# types.
include("samplepath.jl")
# basisfunctions.jl implements the Fourier and Faber-Schauder basis.
include("basisfunctions.jl")
# SDE.jl implements AbstractSDE, SDE and extends
#Distribution.rand to sample an SDE sample path.
include("sde.jl")
# GaussianVectors.jl implements Gaussian vectors in a flexible way.
include("gaussianvectors.jl")
# gaussianprocess.jl implements the abstract type AbstractGaussianProcess with
# subtypes GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
# Extends Base.length and Distributions.rand to AbstractGaussianProcess types.
include("gaussianprocess.jl")
# calculateposterior.jl implements calculateposterior methods to calculate the
# posterior.
include("calculateposterior.jl")

#samplepath.jl are the same in both versions.
#basisfunctions.jl are the same in both versions.
#sde.jl are the same in both versions.
#gaussianvectors.jl are the same in both versions.
#gaussianprocess.jl are the same in both versions.
#calculateposterior.jl are the same in both versions.


text:
function writecbfdata(filename,dat::CBFData, comments="")
    if endswith(filename,"cbf.gz")
        fd = gzopen(filename,"w")
    else
        @assert endswith(filename, "cbf")
        fd = open(filename,"w")
    end

    if comments == ""
        println(fd,"# Generated by ConicBenchmarkUtilities.jl")
    else
        println(fd,comments)
    end
    println(fd,"VER\n2\n")

    println(fd,"OBJSENSE")
    if dat.sense == :Min
        println(fd,"MIN")
    else
        @assert dat.sense == :Max
        println(fd,"MAX")
    end
    println(fd)

    if length(dat.psdvar) > 0
        println(fd, "PSDVAR")
        println(fd, length(dat.psdvar))
        for v in dat.psdvar
            println(fd, v)
        end
        println(fd)
    end

    println(fd,"VAR")
    println(fd,dat.nvar, " ", length(dat.var))
    for (cone,nvar) in dat.var
        println(fd, cone, " ", nvar)
    end
    println(fd)

    if length(dat.intlist) > 0
        println(fd, "INT")
        println(fd, length(dat.intlist))
        for k in dat.intlist
            println(fd, k-1)
        end
        println(fd)
    end

    if length(dat.psdcon) > 0
        println(fd, "PSDCON")
        println(fd, length(dat.psdcon))
        for v in dat.psdcon
            println(fd, v)
        end
        println(fd)
    end

    println(fd,"CON")
    println(fd,dat.nconstr, " ", length(dat.con))
    for (cone,ncon) in dat.con
        println(fd, cone, " ", ncon)
    end
    println(fd)

    if length(dat.objfcoord) > 0
        println(fd,"OBJFCOORD")
        println(fd,length(dat.objfcoord))
        for (a,b,c,d) in dat.objfcoord
            println(fd, a-1, " ", b-1, " ", c-1, " ", d)
        end
        println(fd)
    end

    if length(dat.objacoord) > 0
        println(fd,"OBJACOORD")
        println(fd,length(dat.objacoord))
        for (i,v) in dat.objacoord
            println(fd,i-1," ",v)
        end
        println(fd)
    end

    if length(dat.fcoord) > 0
        println(fd,"FCOORD")
        println(fd,length(dat.fcoord))
        for (a,b,c,d,e) in dat.fcoord
            println(fd,a-1," ",b-1," ",c-1," ",d-1," ",e)
        end
        println(fd)
    end

    if length(dat.acoord) > 0
        println(fd,"ACOORD")
        println(fd,length(dat.acoord))
        for (i,j,v) in dat.acoord
            println(fd,i-1," ",j-1," ",v)
        end
        println(fd)
    end

    if length(dat.hcoord) > 0
        println(fd,"HCOORD")
        println(fd,length(dat.hcoord))
        for (a,b,c,d,e) in dat.hcoord
            println(fd,a-1," ",b-1," ",c-1," ",d-1," ",e)
        end
        println(fd)
    end

    if length(dat.dcoord) > 0
        println(fd,"DCOORD")
        println(fd,length(dat.dcoord))
        for (a,b,c,d) in dat.dcoord
            println(fd, a-1, " ", b-1, " ", c-1, " ", d)
        end
        println(fd)
    end

    if length(dat.bcoord) > 0
        println(fd,"BCOORD")
        println(fd,length(dat.bcoord))
        for (i,v) in dat.bcoord
            println(fd,i-1," ",v)
        end
        println(fd)
    end

    close(fd)
end

max_stars_repo_name: mlubin/ConicBenchmarkUtilities

max_issues_repo_name: mlubin/ConicBenchmarkUtilities

max_forks_repo_name: mlubin/ConicBenchmarkUtilities


text:
using ApproxFun, BandedMatrices, BlockBandedMatrices, LinearAlgebra, BlockArrays, MultivariateOrthogonalPolynomials, Test
using SparseVolterraExamples

#############################
## Tests for Volterra integral equations of first and second kind
## with Clenshaw Jacobi expansion of the kernel.
####
## Includes tests for 1 to x and 1 to 1-x limits of integration.
#############################

@testset "Volterra Integral Equation 1st kind, could be monomial, flip=true" begin
    g(x) = sin(x)-x*cos(x);
    K(x,y) = y;
    @test triVolterraEQ1FullKernelSolver(g,K,30,true) ≈ Fun(x->sin(x),Jacobi(0,1, 0..1))
    #########################################################################
    g(x) = 2+x-2*exp(x)+x*exp(x)
    K(x,y) = x-y;
    @test triVolterraEQ1FullKernelSolver(g,K,20,true) ≈ Fun(x->x*exp(x),Jacobi(0,1, 0..1))
end

@testset "Volterra Integral Equation 1st kind, could be monomial, flip=false" begin
    N = 20
    g(x) = 3-x-2*exp(1-x)+(1-x)*exp(1-x);
    K(x,y) = 1-y-x;
    @test triVolterraEQ1FullKernelSolver(g,K,N,false) ≈ Fun(x -> x*exp(x), Jacobi(0,1,0..1))
end

@testset "Volterra Integral Equation 1st kind, not monomial, flip=false" begin
    N = 20;
    g(x) = 1/12*(4-4*(-1+x)^3+3*(-1+x)^4-4*cos((1-x)^3))
    K(x,y) = 1+y+sin(y^3);
    @test triVolterraEQ1FullKernelSolver(g,K,N,false) ≈ Fun(x->x^2,Jacobi(0,1,0..1))
end

@testset "Volterra Integral Equation 2nd kind, could be monomial, flip=false" begin
    N = 15;
    S = Jacobi(0,1,0..1);
    g = Fun(x -> (x-1)*cos(2*pi*(1-x))/(2*pi)+((4*pi^2+1)/(4*pi^2))*sin(2*pi*(1-x)), S);
    K(x,y)=y;
    @test triVolterraEQ2FullKernelSolver(g,K,N,false) ≈ Fun(x -> sin(2*pi*(1-x)), S)
    #########################################################################
    N = 15;
    S = Jacobi(0,1,0..1);
    g = Fun(x -> (-17/30)*x^6+2*x^5-(3/2)*x^4-(11/2)*x^3+(15/2)*x^2-(21/10)*x+1/6, S);
    K(x,y) = 7*x*y+5*y^2;
    @test triVolterraEQ2FullKernelSolver(g,K,N,false) ≈ Fun(x -> 1-x^3, S)
end

@testset "Volterra Integral Equation 2nd kind, not monomial, flip=true" begin
    g = Fun(x->x^2+2*x-exp(x)+2, Jacobi(0,1,0..1));
    K(x,y) = (y-x)^2;
    @test triVolterraEQ2FullKernelSolver(g,K,20,true) ≈ Fun(x->exp(x),Jacobi(0,1,0..1));
end


text:
using DataFrames, QuasiCopula, LinearAlgebra, GLM, RDatasets
using Test

Mmmec = dataset("mlmRev", "Mmmec");
df = Mmmec
y = :Deaths
grouping = :Region
covariates = [:UVB]

d = Poisson()
link = LogLink()

# forming CS model with Poisson base
Poisson_CS_model = CS_model(df, y, grouping, covariates, d, link)

@test typeof(Poisson_CS_model) == GLMCopulaCSModel{Float64, Poisson{Float64}, LogLink}
@test Poisson_CS_model.penalized == false

# forming CS model with Poisson base with L2 ridge penalty
Poisson_CS_model_penalized = CS_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(Poisson_CS_model_penalized) == GLMCopulaCSModel{Float64, Poisson{Float64}, LogLink}
@test Poisson_CS_model_penalized.penalized == true

# intercept only model
Poisson_CS_model_intercept_only = CS_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(Poisson_CS_model_intercept_only.data[1].X) == Poisson_CS_model_intercept_only.data[1].n

#### nb ####
d = NegativeBinomial()
link = LogLink()
nb_CS_model = CS_model(df, y, grouping, covariates, d, link)

@test typeof(nb_CS_model) == NBCopulaCSModel{Float64, NegativeBinomial{Float64}, LogLink}
@test nb_CS_model.penalized == false

# forming CS model with nb base with L2 ridge penalty
nb_CS_model_penalized = CS_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(nb_CS_model_penalized) == NBCopulaCSModel{Float64, NegativeBinomial{Float64}, LogLink}
@test nb_CS_model_penalized.penalized == true

# intercept only model
nb_CS_model_intercept_only = CS_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(nb_CS_model_intercept_only.data[1].X) == nb_CS_model_intercept_only.data[1].n

#### gaussian ####
d = Normal()
link = IdentityLink()
gaussian_CS_model = CS_model(df, y, grouping, covariates, d, link)

@test typeof(gaussian_CS_model) == GaussianCopulaCSModel{Float64}
@test gaussian_CS_model.penalized == false

# forming CS model with gaussian base with L2 ridge penalty
gaussian_CS_model_penalized = CS_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(gaussian_CS_model_penalized) == GaussianCopulaCSModel{Float64}
@test gaussian_CS_model_penalized.penalized == true

# intercept only model
gaussian_CS_model_intercept_only = CS_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(gaussian_CS_model_intercept_only.data[1].X) == gaussian_CS_model_intercept_only.data[1].n


text:
using LatticeBasisReduction
using Base.Test

# Custom rounding
@test LatticeBasisReduction.round_ties_down(1.3) == 1
@test LatticeBasisReduction.round_ties_down(1.5) == 1
@test LatticeBasisReduction.round_ties_down(1.6) == 2
@test LatticeBasisReduction.round_ties_down(2) == 2
@test LatticeBasisReduction.round_ties_down(-0.5) == -1

# Verify gram schmidt process with example 3.3 from Bremner
eps = 1e-3
@test_approx_eq_eps gram_schmidt([3 -1 5; -5 2 -1; -3 9 2]) [3 -1 5; -109/35 48/35 15/7; 1521/566 1859/283 -169/566] eps

# Verify LLL algorithm
x = [[-2 7 7 -5]
     [3 -2 6 -1]
     [2 -8 -9 -7]
     [8 -9 6 -4]]

@test_approx_eq_eps LLL(x, 1) [2 3 1 1; 2 0 -2 -4; -2 2 3 -3; 3 -2 6 -1] eps
@test_approx_eq_eps LLL(x, 1.) [2 3 1 1; 2 0 -2 -4; -2 2 3 -3; 3 -2 6 -1] eps
@test_approx_eq_eps LLL(Array{Float64}(x), 1) [2 3 1 1; 2 0 -2 -4; -2 2 3 -3; 3 -2 6 -1] eps

# Verify diophantine approximation for a single number
@test diophantine(float(π)) == 355//113

# Verify gauss for minimal basis of a lattice in R^2
x = [-56, 43]; y = [95, -73]
@test gauss(x,y) == ([1,1],[-1,2])
@test gauss(y,x) == ([1,1],[-1,2])

# Verify Simulataneous Diophantine Approximation
nums, dem = sda(sqrt(2), sqrt(3), sqrt(5), sqrt(7), ε = eps)
@test_approx_eq_eps nums[1]/dem sqrt(2) eps
@test_approx_eq_eps nums[2]/dem sqrt(3) eps
@test_approx_eq_eps nums[3]/dem sqrt(5) eps
@test_approx_eq_eps nums[4]/dem sqrt(7) eps


text:
## context

## A context stores objects where named values are looked up.
mutable struct Context
    view ## of what? a Dict, Module, CompositeKind, DataFrame.parent.
    parent ## a context or nothing
    _cache::Dict
end

function Context(view, parent=nothing)
    d = Dict()
    d["."] = isa(view, AnIndex) ? view.value : view
    Context(view, parent, d)
end

## used by renderTokens
function ctx_push(ctx::Context, view)
    Context(view, ctx) ## add context as a parent
end

function ctx_pop(ctx::Context)
    ctx.parent
end

# we have some rules here
# * Each part of a dotted name should resolve only against its parent.
# * Any falsey value prior to the last part of the name should yield ''.
# * The first part of a dotted name should resolve as any other name.
function lookup_dotted(ctx::Context, dotted)
    for key in split(dotted, ".")
        nctx = lookup(Context(ctx.view), key)
        falsy(nctx) && return nothing
        ctx = Context(nctx, ctx)
    end
    ctx.view
end

# look up key in context
function lookup(ctx::Context, key)
    if haskey(ctx._cache, key)
        return ctx._cache[key]
    end

    # use global lookup down
    global_lookup = false
    if startswith(key, "~")
        global_lookup = true
        key = replace(key, r"^~" => "")
    end

    # begin
    context = ctx
    value = nothing
    while context !== nothing
        value′ = nothing
        ## does name have a .?
        if occursin(r"\.", key)

            value′ = lookup_dotted(context, key)

            if value′ == nothing
                ## do something with "."
                ## we use .[ind] to refer to value in parent of given index;
                m = match(r"^\.\[(.*)\]$", key)
                m === nothing && break

                idx = m.captures[1]
                vals = context.parent.view

                # this has limited support for indices: "end", or a number, but no
                # arithmetic, such as `end-1`.
                ## This is for an iterable; rather than
                ## limit the type, we let non-iterables error.
                if true # isa(vals, AbstractVector)  || isa(vals, Tuple) # supports getindex(v, i)?
                if idx == "end"
                    value′ = AnIndex(-1, vals[end])
                else
                    ind = Base.parse(Int, idx)
                    value′ = AnIndex(ind, string(vals[ind]))
                end
                    #break
                end
            end
        else
            value′ = lookup_in_view(context.view, stripWhitespace(key))
        end

        if value′ !== nothing
            value = value′
            !global_lookup && break
        end

        context = context.parent
    end

    ## cache
    ctx._cache[key] = value
    return(value)
end

## Lookup value in an object by key
## This of course varies based on the view.
## After checking several specific types, if view is Tables compatible
## this will return "column" corresponding to the key
function lookup_in_view(view, key)

    val = _lookup_in_view(view, key)
    !falsy(val) && return val

    if Tables.istable(view)
        isempty(Tables.rows(view)) && return nothing
        sch = Tables.schema(Tables.rows(view))
        falsy(sch) && return nothing
        k = normalize(key)
        if k ∈ sch.names
            return [row[k] for row ∈ Tables.rows(view)]
        end
    # elseif  is_dataframe(view)

    #     if occursin(r"^:", key)  key = key[2:end] end
    #     key = Symbol(key)
    #     out = nothing
    #     if haskey(view, key)
    #             out = view[1, key] ## first element only
    #     end
    #     out
    else
        __lookup_in_view(view, key)
    end
end

# look up key in view, return `nothing` if not found
function _lookup_in_view(view::AbstractDict, key)
    get(view, normalize(key), nothing)
end

# support legacy use of `first` and `second` as variable names
# referring to piece, otherwise look up value
function _lookup_in_view(view::Pair, key)
    ## is it a symbol?
    key == "first" && return view.first
    key == "second" && return view.second
    view.first == normalize(key) ? view.second : nothing
end

function _lookup_in_view(view::NamedTuple, key)
    get(view, normalize(key), nothing)
end

function _lookup_in_view(view::Module, key)

    hasmatch = false
    re = Regex("^$key\$")
    for i in names(view, all=true)
        if occursin(re, string(i))
            hasmatch = true
            break
        end
    end

    out = nothing
    if hasmatch
        out = getfield(view, Symbol(key))  ## view.key
    end
    out

end

_lookup_in_view(view, key) = nothing

## Default lookup is likely not great,
function __lookup_in_view(view, key)

    k = normalize(key)

    # check propertyname, then fieldnames
    if k in propertynames(view)
        getproperty(view, k)
    else

        nms = fieldnames(typeof(view))

        # we match on symbol or string for fieldname
        if isa(k, Symbol)
            has_match = k in nms
        else
            has_match = Symbol(k) in nms
        end

        out = nothing
        if has_match
            out = getfield(view, Symbol(k))  ## view.key
        end

        out
    end
end


text:
# Understanding what the compiler does...

f(k) = 5k - 1

function g(k)
    for i=1:10
        k = f(k)
    end
    return k
end

r = g(1.)

code_native(g,(Float64,))

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
"""
    convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T

converts compact domain x ∈ [a,b] to compact domain out ∈ [c,d].
"""
function convertcompactdomain(x::T, a::T, b::T, c::T, d::T)::T where T <: Real

    return (x-a)*(d-c)/(b-a)+c
end

function convertcompactdomain(x::Vector{T}, a::T, b::T, c::T, d::T)::Vector{T} where T <: Real

    return collect( convertcompactdomain(x[i], a, b, c, d) for i = 1:length(x) )
end

#unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))

function runNLopt!(  opt,
    p0::Vector{T},
    obj_func,
    grad_func,
    p_lbs,
    p_ubs;
    max_iters = 10000,
    xtol_rel = 1e-12,
    ftol_rel = 1e-12,
    maxtime = Inf) where T

    @assert length(p0) == length(p_lbs) == length(p_ubs)

    opt.maxeval = max_iters
    opt.lower_bounds = p_lbs
    opt.upper_bounds = p_ubs
    opt.xtol_rel = xtol_rel
    opt.ftol_rel = ftol_rel
    opt.maxtime = maxtime

    opt.min_objective = (xx, gg)->genericNLoptcostfunc!(xx, gg, obj_func, grad_func)

    # optimize.
    (minf, minx, ret) = NLopt.optimize(opt, p0)

    N_evals = opt.numevals

    return minf, minx, ret, N_evals
end

function genericNLoptcostfunc!(x::Vector{T}, df_x, f, df)::T where T <: Real

    if length(df_x) > 0
        df_x[:] = df(x)
    end

    return f(x)
end


text:
import Random

#
# In this test we compute the average displacement of the x coordinates of the particles
#
function test1(;N=100_000,parallel=true,x=nothing)

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  end

  # Initialize auxiliary linked lists
  cl = CellList(x,box,parallel=parallel)

  # Function to be evalulated for each pair: sum of displacements on x
  f(x,y,avg_dx) = avg_dx + abs(x[1] - y[1])

  avg_dx = (N/(N*(N-1)/2)) * map_pairwise!(
    (x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),
    0.,box,cl,
    parallel=parallel
  )

  correct_result = 100.5799271448272
  return avg_dx ≈ correct_result, avg_dx
end

#
# In this test we compute the histogram of distances, expected to follow the
# function f(f) = ρ(4/3)π(r[i+1]^3 - r[i]^3) with ρ being the density of the system.
#
function test2(;N=100_000,parallel=true,x=nothing)

  # Number of particles, sides and cutoff
  sides = SVector{3,Float64}(250,250,250)
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  end

  # Initialize auxiliary linked lists
  cl = CellList(x,box,parallel=parallel)

  # Function to be evalulated for each pair: build distance histogram
  function build_histogram!(d2,hist)
    d = sqrt(d2)
    ibin = floor(Int,d) + 1
    hist[ibin] += 1
    return hist
  end

  # Preallocate and initialize histogram
  hist = zeros(Int,10)

  # Run pairwise computation
  hist = (N/(N*(N-1)/2)) * map_pairwise!(
    (x,y,i,j,d2,hist) -> build_histogram!(d2,hist),
    hist,box,cl,
    parallel=parallel
  )
  return hist

end

#
# In this test we compute the "gravitational potential", pretending that each particle
# has a different mass. In this case, the closure is used to pass the masses to the
# function that computes the potential.
#
function test3(;N=100_000,parallel=true,x=nothing)

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  end

  # Initialize auxiliary linked lists
  cl = CellList(x,box,parallel=parallel)

  # masses
  mass = [ 5 * x[i][1] for i in 1:N ]

  # Function to be evalulated for each pair: build distance histogram
  function potential(i,j,d2,u,mass)
    d = sqrt(d2)
    u = u - 9.8*mass[i]*mass[j]/d
    return u
  end

  # Run pairwise computation
  u = map_pairwise!(
    (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),
    0.0,box,cl,
    parallel=parallel
  )
  return u

end

#
# In this test we compute the "gravitational force", pretending that each particle
# has a different mass. In this case, the closure is used to pass the masses and
# the force vector to the function that computes the potential.
#
function test4(;N=100_000,parallel=true,x=nothing)

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  end

  # masses
  mass = [ 5 * x[i][1] for i in 1:N ]

  # Initialize auxiliary linked lists
  cl = CellList(x,box,parallel=parallel)

  # Function to be evalulated for each pair
  function calc_forces!(x,y,i,j,d2,mass,forces)
    G = 9.8*mass[i]*mass[j]/d2
    d = sqrt(d2)
    df = (G/d) * (x - y)
    forces[i] = forces[i] - df
    forces[j] = forces[j] + df
    return forces
  end

  # Preallocate and initialize forces
  forces = [ zeros(SVector{3,Float64}) for i in 1:N ]

  # Run pairwise computation
  forces = map_pairwise!(
    (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces),
    forces,box,cl,
    parallel=parallel
  )
  return forces

end

#
# In this test we compute the minimum distance between two independent sets of particles
#
function test5(;N1=1_500,N2=1_500_000,parallel=true,x=nothing,y=nothing)

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N1 ]
  end
  if y === nothing
    y = [ sides .* rand(SVector{3,Float64}) for i in 1:N2 ]
  end

  # Initialize auxiliary linked lists (largest set!)
  cl = CellList(x,y,box,parallel=parallel)

  # Function that keeps the minimum distance
  f(i,j,d2,mind) = d2 < mind[3] ? (i,j,d2) : mind

  # We have to define our own reduce function here (for the parallel version)
  function reduce_mind(output,output_threaded)
    mind = output_threaded[1]
    for i in 2:nthreads()
      if output_threaded[i][3] < mind[3]
        mind = output_threaded[i]
      end
    end
    return mind
  end

  # Initialize
  mind = ( 0, 0, +Inf )

  # Run pairwise computation
  mind = map_pairwise!(
    (x,y,i,j,d2,mind) -> f(i,j,d2,mind),
    mind,box,cl;reduce=reduce_mind, parallel=parallel
  )

  # Take the square root of the minimum distance to return
  return (mind[1],mind[2],sqrt(mind[3]))

end

#
# In this test we compute the minimum distance between two independent sets of particles,
# without periodic conditions
#
function test6(;N1=1_500,N2=1_500_000,parallel=true,x=nothing,y=nothing)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ rand(SVector{3,Float64}) for i in 1:N1 ]
  else
    N1 = length(x)
  end
  if y === nothing
    y = [ rand(SVector{3,Float64}) for i in 1:N2 ]
  else
    N2 = length(y)
  end

  # Boundaries
  xmin = [ +Inf, +Inf, +Inf ]
  xmax = [ -Inf, -Inf, -Inf ]
  for v in x
    @. xmin = min(xmin,v)
    @. xmax = max(xmax,v)
  end
  for v in y
    @. xmin = min(xmin,v)
    @. xmax = max(xmax,v)
  end

  # Obtain one upper bound for dmin by computing one distance for each element
  # of the smallest vector
  cutoff = +Inf
  for v in x
    iy = rand(1:N2)
    cutoff = min(CellListMap.norm(v-y[iy]),cutoff)
  end

  # Define box sides
  sides = (xmax - xmin) .+ cutoff
  box = Box(sides,cutoff)

  # Initialize auxiliary linked lists (largest set!)
  cl = CellList(x,y,box,parallel=parallel)

  # Function that keeps the minimum distance
  f(i,j,d2,mind) = d2 < mind[3] ? (i,j,d2) : mind

  # We have to define our own reduce function here (for the parallel version)
  function reduce_mind(output,output_threaded)
    mind = output_threaded[1]
    for i in 2:Threads.nthreads()
      if output_threaded[i][3] < mind[3]
        mind = output_threaded[i]
      end
    end
    return mind
  end

  # Initialize
  mind = ( 0, 0, +Inf )

  # Run pairwise computation
  mind = map_pairwise!(
    (x,y,i,j,d2,mind) -> f(i,j,d2,mind),
    mind,box,cl;reduce=reduce_mind,parallel=parallel
  )

  # Take the square root of the minimum distance to return
  return (mind[1],mind[2],sqrt(mind[3]))

end

#
# In this test we compute the complete neighbour list of particles, meaning all the pairs
# that are within the cutoff distance
#
function test7(;N=100_000,parallel=true,x=nothing)

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  Random.seed!(321)
  if x === nothing
    x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  end

  # Initialize auxiliary linked lists
  cl = CellList(x,box,parallel=parallel)

  # Function to be evalulated for each pair: push pair if d<cutoff
  function push_pair!(i,j,d2,pairs,cutoff)
    d = sqrt(d2)
    if d < cutoff
      push!(pairs,(i,j,d))
    end
    return pairs
  end

  # Reduction function
  function reduce_pairs(pairs,pairs_threaded)
    pairs = pairs_threaded[1]
    for i in 2:nthreads()
      append!(pairs,pairs_threaded[i])
    end
    return pairs
  end

  # Initialize output
  pairs = Tuple{Int,Int,Float64}[]

  # Run pairwise computation
  pairs = map_pairwise!(
    (x,y,i,j,d2,pairs) -> push_pair!(i,j,d2,pairs,cutoff),
    pairs,box,cl,
    reduce=reduce_pairs,
    parallel=parallel
  )
  return pairs

end

#
# florpi
#
function florpi(;N=100_000,cd=true,parallel=true)

  @inline dot(x::SVector{3,Float64},y::SVector{3,Float64}) = x[1]*y[1] + x[2]*y[2] + x[3]*y[3]

  function compute_pairwise_mean_cell_lists!(x,y,i,j,d2,hist,velocities,rbins)
    d = x - y
    r = sqrt(d2)
    ibin = searchsortedfirst(rbins, r) - 1
    hist[1][ibin] += 1
    hist[2][ibin] += dot(velocities[i]-velocities[j],d)/r
    return hist
  end

  function reduce_hist(hist,hist_threaded)
    for i in 1:Threads.nthreads()
      hist[1] .+= hist_threaded[i][1]
      hist[2] .+= hist_threaded[i][2]
    end
    return hist
  end

  n_halos = N

  if cd
    density = 10^5/250^3  # density of the original problem
    boxsize = (n_halos / density)^(1/3)
  else
    boxsize = 250.
  end

  Random.seed!(321)
  Lbox = [boxsize,boxsize,boxsize]
  positions = boxsize .* rand(Float64, 3, n_halos)
  velocities = rand(Float64, 3, n_halos)
  rbins = [0.,2.,4.,6.,8.,10.]
  r_max = maximum(rbins)

  n = size(positions)[2]
  positions = reshape(reinterpret(SVector{3,Float64},positions),n)
  velocities = reshape(reinterpret(SVector{3,Float64},velocities),n)

  box = Box(Lbox, r_max, UnitCellType=OrthorhombicCell, lcell=1)
  cl = CellList(positions,box,parallel=parallel)
  hist = (zeros(Int,length(rbins)-1), zeros(Float64,length(rbins)-1))

  # Needs this to stabilize the type of velocities and hist, probably
  function barrier(f!,velocities,rbins,hist,box,cl,reduce_hist,parallel)
    hist = map_pairwise!(
      (x,y,i,j,d2,hist) -> f!(x,y,i,j,d2,hist,velocities,rbins),
      hist, box, cl,
      reduce=reduce_hist,
      parallel=parallel
    )
    return hist
  end

  hist = barrier(compute_pairwise_mean_cell_lists!,
    velocities,rbins,hist,box,cl,reduce_hist,parallel)

  n_pairs = hist[1]
  mean_v_r = hist[2]
  mean_v_r[n_pairs .> 0] = mean_v_r[n_pairs .> 0]./n_pairs[n_pairs .> 0]

  correct_result = [
   0.0007883474482652579
  -0.0035662371878635722
  -0.00040742008823982926
   0.0003623877989466509
  -0.0010441334538614498
  ]
  return mean_v_r ≈ correct_result, mean_v_r

end


text:
# We explore the batch and fidelity here

@everywhere include("../INCLUDEME.jl")

@everywhere using WFLearning

ptask("nqubits-2", 8, 10, 40, 50, 0.01)

# @everywhere function main(nlayers; nbatch=10, epochs=200, learning_rate=0.1)
#     pmap(n->task("nqubits", n, nlayers, nbatch, epochs, learning_rate), 5:7)
# end
#
# main(5, epochs=200, learning_rate=0.1, nbatch=10)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using InteractBase, Colors
using CSSUtil
using InteractBulma, Blink
using Widgets
settheme!(Bulma())
using WebIO
w = Window()

#---
f = filepicker(label = "Upload");
body!(w, f)
observe(f)
#---

f = filepicker(multiple = true, accept = ".csv");
display(f)
observe(f)
#---
s = InteractBase.input("Write here")
body!(w, s)
observe(s)
#---
s = textbox("Write here", value = "A", className="is-danger")
body!(w, s)
observe(s)
#---
s = autocomplete(["Opt 1", "Option 2", "Opt 3"], "Write here")
body!(w, s)
observe(s)
opentools(w)
#---
v = Observable(RGB(1,0,0))
s = colorpicker(v)
body!(w, s)
observe(s)
#---
s1 = slider(1:20, style = Dict("width"=>"400px"))
sobs = observe(s1)
body!(w, vbox(s1, sobs));
#---
s1 = slider(vcat(0, logspace(-2, 2)), label = "Log slider")
sobs = observe(s1)
body!(w, vbox(s1, sobs));
#---
button1 = button("button one {{clicks}}")
num_clicks = observe(button1)
button2 = button("button two {{clicks}}", value = num_clicks)
body!(w, hbox(button1, button2, num_clicks));
#---
body!(w, button(label = "Press!"))
#---
v = checkbox(label = "Agree")
body!(w, v)
observe(v)
#---
v = toggle(true, "Agree", className="is-danger")
body!(w, v)
observe(v)
#---
v = checkboxes(["A", "B", "C"]);
body!(w, v)
observe(v)
#---
v = toggles(["A", "B", "C"], className="is-danger");
body!(w, v)
observe(v)
#---
using WebIO, Blink, Observables

width, height = 700, 300
colors = ["black", "gray", "silver", "maroon", "red", "olive", "yellow", "green", "lime", "teal", "aqua", "navy", "blue", "purple", "fuchsia"]
_color(i) = colors[i%length(colors)+1]
ui = @manipulate for nsamples in 1:200,
        sample_step in slider(0.01:0.01:1.0, value=0.1, label="sample step"),
        phase in slider(0:0.1:2pi, value=0.0, label="phase"),
        radii in 0.1:0.1:60,
        show in true,
        s in Observable(true)
    cxs_unscaled = [i*sample_step + phase for i in 1:nsamples]
    cys = sin.(cxs_unscaled) .* height/3 .+ height/2
    cxs = cxs_unscaled .* width/4pi
    show ? dom"svg:svg[width=$width, height=$height]"(
        (dom"svg:circle[cx=$(cxs[i]), cy=$(cys[i]), r=$radii, fill=$(_color(i))]"()
            for i in 1:nsamples)...
    ) : dom"div"("Nothing to see here")
end
body!(w, ui)
opentools(w)
#---

using InteractBase, Plots, CSSUtil, DataStructures

x = y = 0:0.1:30

freqs = OrderedDict(zip(["pi/4", "π/2", "3π/4", "π"], [π/4, π/2, 3π/4, π]))

mp = @manipulate for freq1 in freqs, freq2 in slider(0.01:0.1:4π; label="freq2")
    y = @. sin(freq1*x) * sin(freq2*x)
    plot(x, y)
end
body!(w, mp)

#---
s = dropdown(["a1", "a2nt", "a3"], label = "test")
body!(w, s)
observe(s)
#---

s = togglebuttons(["a1", "a2nt", "a3"], label = "x");
body!(w, s)
observe(s)
#---

s = radiobuttons(["a1", "a2nt", "a3"]);
display(s)
observe(s)
#---
# IJulia
ui = s
display(ui);
# Mux
using Mux
webio_serve(page("/", req -> ui))
#---


text:
"""
    add_objective!(
        community::CoreModel,
        objective_mets::Vector{String};
        objective_weights = Float64[],
        objective_column_index = 0,
    )

Add an objective to the `community` model. Supply the string names of the
objective metabolites in `objective_mets`. Optionally specify the weight to
assign each metabolite in the objective function, if unassigned then equal
weight is assumed. Also, optionally specify whether the objective already exists
in the model by assigning `objective_column_index`. If unassigned then an
objective column will be added, otherwise the column at `objective_column_index`
will be updated.

Note, the weights are negated inside the function so that the objective metabolites
are seen as reagents/substrates, not products in the reaction equation.

# Example
```
add_objective!(model, ["met1", "met2"]) # adds a new column with weights = [1,1]
add_objective!(model, ["met1", "met2"]; objective_weights=[0.1, 0.9]) # adds a new column
add_objective!(model, ["met1", "met2"]; objective_weights=[0.1, 0.9], objective_column_index=10) # updates column 10
```
"""
function add_objective!(
    community::CoreModel,
    objective_mets::Vector{String};
    objective_weights = Float64[],
    objective_column_index = 0,
)
    obj_inds = indexin(objective_mets, metabolites(community))
    if isempty(objective_weights)
        objective_weights = repeat([1.0], inner = length(objective_mets))
    end

    if objective_column_index == 0 # needs to be created
        nr, _ = size(community.S)
        objcol = spzeros(nr)
        objcol[obj_inds] .= -objective_weights

        # extend model by one reaction
        community.S = hcat(community.S, objcol)
        community.xl = [community.xl; 0.0]
        community.xu = [community.xu; 1000.0]
        community.rxns = [community.rxns; "community_biomass"]
        community.c = spzeros(size(community.S, 2))
        community.c[end] = 1.0
    else # only update
        community.S[:, objective_column_index] .= 0.0 # reset
        community.S[obj_inds, objective_column_index] .= -objective_weights
        dropzeros!(community.S)
        community.c = spzeros(size(community.S, 2))
        community.c[objective_column_index] = 1.0
        community.xl[objective_column_index] = 0.0
        community.xu[objective_column_index] = 1000.0
    end
end

"""
    add_model_with_exchanges(
        community::CoreModel,
        model::MetabolicModel,
        exchange_rxn_ids::Vector{String},
        exchange_met_ids::Vector{String};
        model_name = "unknown_species",
        biomass_id = nothing,
    )::CoreModel

Add `model` to `community`, which is a pre-existing community model with
`exchange_rxn_ids` and `exchange_met_ids`. The `model_name` is appended to
each reaction and metabolite, see [`join`](@ref). If `biomass_id` is specified
then a biomass metabolite for `model` is also added to the resulting model. The
column corresponding to the `biomass_id` reaction then produces this new biomass
metabolite with unit coefficient. Note, `exchange_rxn_ids` and
`exchange_met_ids` must already exist in the `community` model.

# Example
```
community = add_model_with_exchanges(community, model, exchange_rxn_ids, exchange_met_ids; model_name="species_2", biomass_id="BIOMASS_Ecoli_core_w_GAM")
```
"""
function add_model_with_exchanges(
    community::CoreModel,
    model::MetabolicModel,
    exchange_rxn_ids::Vector{String},
    exchange_met_ids::Vector{String};
    model_name = "unknown_species",
    biomass_id = nothing,
)::CoreModel

    exchange_met_community_inds = indexin(exchange_met_ids, metabolites(community))
    exchange_rxn_community_inds = indexin(exchange_rxn_ids, reactions(community))
    if any(isnothing.(exchange_met_community_inds)) ||
       any(isnothing.(exchange_rxn_community_inds))
        throw(
            DomainError(
                "exchange metabolite/reaction not found.",
                "Exchange metabolites/reactions must already be contained in the community model.",
            ),
        )
    end

    n_cmodel_rows, n_cmodel_cols = size(stoichiometry(community))
    n_model_rows, n_model_cols = size(stoichiometry(model))
    # A note on the variable names here.Suppose M is some sparse matrix, then I
    # = row indices, J = column indices and V = values at the associated
    # indices. So I[a] = i, J[a]=j and then M[i,j] = V[a]
    Iadd, Jadd, Vadd = findnz(stoichiometry(model))

    # shift to fit into community
    Iadd .+= n_cmodel_rows
    Jadd .+= n_cmodel_cols

    # when adding a single model not that many reactions, push! okay?
    exchange_rxn_model_idxs = indexin(exchange_rxn_ids, reactions(model))
    for i = 1:length(exchange_rxn_ids)
        isnothing(exchange_rxn_model_idxs[i]) && continue
        push!(Iadd, exchange_met_community_inds[i]) # already present ex met in community model
        push!(Jadd, n_cmodel_cols + exchange_rxn_model_idxs[i]) # new column hence the offset
        push!(Vadd, 1.0)
    end

    biomass_met = 0.0
    if biomass_id != "" # add biomass metabolite
        biomass_rxn = first(indexin([biomass_id], reactions(model)))
        push!(Iadd, n_model_rows + n_cmodel_rows + 1)
        push!(Jadd, biomass_rxn + n_cmodel_cols)
        push!(Vadd, 1.0)
        biomass_met = 1
    end

    n_metabolites_total = n_model_rows + n_cmodel_rows + biomass_met
    n_reactions_total = n_cmodel_cols + n_model_cols

    I, J, V = findnz(stoichiometry(community))
    I = [I; Iadd]
    J = [J; Jadd]
    V = [V; Vadd]
    S = sparse(I, J, V, n_metabolites_total, n_reactions_total)

    # A note on the variables here. The bounds are vectors of upper and lower
    # bounds for each reaction. So lbs = [lb_1, lb_2, lb_i, ...], ubs = [ub_1,
    # ub_2, ub_i, ...] for reaction i. See the bounds function for more
    # information
    lbsadd, ubsadd = bounds(model)
    lbs, ubs = bounds(community)
    lbs = [lbs; lbsadd]
    ubs = [ubs; ubsadd]

    rxnsadd = "$(model_name)_" .* reactions(model)
    if !isnothing(biomass_id)
        metsadd =
            ["$(model_name)_" .* metabolites(model); "$(model_name)_" * biomass_id]
    else
        metsadd = "$(model_name)_" .* metabolites(model)
    end
    rxns = [reactions(community); rxnsadd]
    mets = [metabolites(community); metsadd]

    # adds to the original community data, could possibly reset?
    I, V = findnz(balance(community))
    b = sparsevec(I, V, n_metabolites_total)
    I, V = findnz(objective(community))
    c = sparsevec(I, V, n_reactions_total)

    return CoreModel(S, b, c, lbs, ubs, rxns, mets)
end

"""
    join_with_exchanges(models::Vector{M},
        exchange_rxn_ids::Vector{String},
        exchange_met_ids::Vector{String};
        add_biomass_objective=false,
        biomass_ids::Vector{String},
        model_names=String[]
    )

Return a `CoreModel` representing the community model of `models` joined through
their `exchange_rxn_ids` and `exchange_met_ids`. These exchange reactions and
metabolites link to environmental metabolites and reactions.
Optionally specify `model_names` to append a specific name to each reaction
and metabolite of an organism for easier reference (default is `species_i` for
each model index i in `models`). Note, the bounds of the environmental variables
are all set to zero. Thus, to run a simulation you need to constrain them
appropriately. All the other bounds are inherited from the models used to
construct the community model.

If `add_biomass_objective` is true then `biomass_ids` needs to be supplied as
well. This creates a model with an extra reaction added to the end of the
stoichiometric matrix (last column) that can be assigned as the objective
reaction. It also creates biomass "metabolites" that can be used in this
objective reaction. Note, this reaction is unspecified, further action needs to
be taken to specify it, e.g. assign weights to the last column of the
stoichiometric matrix in the rows corresponding to the biomass metabolites.

To further clarify how this `join` works. Suppose you have 2 organisms with
stoichiometric matrices S₁ and S₂ and you want to link them with
`exchange_met_ids = [em₁, em₂, em₃, ...]` and `exchange_rxn_ids = [er₁, er₂,
er₃, ...]`. Then a new community stoichiometric matrix is constructed that looks
like this:
```
            _      er₁  er₂  er₃  ...  b_
           |S₁                           |
           |   S₂                        |
        em₁|                             |
S   =   em₂|                             |
        em₃|                             |
        ...|                             |
        bm₁|                             |
        bm₂|_                           _|

```
The exchange reactions in each model get linked to environmental metabolites,
`emᵢ`, and these get linked to environmental exchanges, `erᵢ`. These `erᵢ`
behave like normal single organism exchange reactions. When
`add_biomass_objective` is true each model's biomass becomes a pseudo-metabolite
(`bmᵢ`). These can be weighted in column `b`, called the `community_biomass`
reaction in the community model, if desired. Refer to the tutorial if this is
unclear.

# Example
```
m1 = load_model(core_model_path)
m2 = load_model(CoreModel, core_model_path)

# need to list ALL the exchanges that will form part of the entire model
exchange_rxn_ids = filter(looks_like_exchange_reaction, reactions(m1))
exchange_met_ids = [first(keys(reaction_stoichiometry(m1, ex_rxn))) for ex_rxn in exchange_rxn_ids]

biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"]

community = join_with_exchanges(
    [m1, m2],
    exchange_rxn_ids,
    exchange_met_ids;
    add_biomass_objective = true,
    biomass_ids = biomass_ids,
)
```
"""
function join_with_exchanges(
    models::Vector{M},
    exchange_rxn_ids::Vector{String},
    exchange_met_ids::Vector{String};
    add_biomass_objective = true,
    biomass_ids = String[],
    model_names = String[],
)::CoreModel where M <: MetabolicModel

    if add_biomass_objective && isempty(biomass_ids)
        throw(
            DomainError(
                "biomass_ids",
                "Please add supply the string ids of the biomass functions when `add_biomass_objective` is true.",
            ),
        )
    end

    if length(exchange_met_ids) != length(exchange_rxn_ids)
        throw(
            DomainError(
                "Exchange identifiers are misspecified",
                "The lenghts of the exchange metabolite and reactions are different.",
            ),
        )
    end

    # get offsets to construct community S
    reaction_lengths = [n_reactions(model) for model in models]
    metabolite_lengths = [n_metabolites(model) for model in models]
    reaction_offset = [0; cumsum(reaction_lengths[1:end-1])]
    metabolite_offset = [0; cumsum(metabolite_lengths[1:end-1])]

    # get each model's S matrix (needed for the size calculations)
    stoichs = [stoichiometry(model) for model in models]
    nnzs = [findnz(stoich) for stoich in stoichs] # nonzero indices per model

    # size calculations
    column_add = add_biomass_objective ? 1 : 0 # objective rxn
    row_add = add_biomass_objective ? length(models) : 0 # biomass as metabolites
    nnz_add = add_biomass_objective ? (1 + length(models)) : 0
    nnz_total =
        sum(length(first(nnz)) for nnz in nnzs) +
        length(models) * length(exchange_rxn_ids) +
        length(exchange_met_ids) + nnz_add
    n_reactions_metabolic = sum(reaction_lengths)
    n_reactions_total = n_reactions_metabolic + length(exchange_rxn_ids) + column_add
    n_metabolites_metabolic = sum(metabolite_lengths)
    n_metabolites_total = n_metabolites_metabolic + length(exchange_met_ids) + row_add

    # Create empty storage vectors
    lbs = spzeros(n_reactions_total)
    ubs = spzeros(n_reactions_total)
    rxns = Array{String,1}(undef, n_reactions_total)
    mets = Array{String,1}(undef, n_metabolites_total)
    I = Array{Int,1}(undef, nnz_total)
    J = Array{Int,1}(undef, nnz_total)
    V = Array{Float64,1}(undef, nnz_total)

    # build metabolic components, block diagonals
    kstart = 1
    for i = 1:length(models)
        kend = kstart + length(nnzs[i][3]) - 1
        rng = kstart:kend
        I[rng] .= nnzs[i][1] .+ metabolite_offset[i]
        J[rng] .= nnzs[i][2] .+ reaction_offset[i]
        V[rng] .= nnzs[i][3]
        kstart += length(nnzs[i][3])
    end

    # build environmental - exchange links
    for i = 1:length(models)
        exchange_rxn_inds = indexin(exchange_rxn_ids, reactions(models[i]))
        exchange_met_inds = indexin(exchange_met_ids, metabolites(models[i]))
        for (n, (ex_rxn, ex_met)) in enumerate(zip(exchange_rxn_inds, exchange_met_inds)) # each exchange rxn has one exchange met
            isnothing(ex_rxn) && continue
            isnothing(ex_met) && continue
            # connect environmental metabolite with exchange metabolite
            I[kstart] = n_metabolites_metabolic + n
            J[kstart] = ex_rxn + reaction_offset[i]
            V[kstart] = -stoichs[i][ex_met, ex_rxn] # ex is normally negative, make positive
            kstart += 1
        end
    end

    # # build diagonal environmental exchanges
    for i = 1:length(exchange_rxn_ids)
        I[kstart] = n_metabolites_metabolic + i
        J[kstart] = n_reactions_metabolic + i
        V[kstart] = -1.0
        kstart += 1
    end

    if add_biomass_objective
        n_before_biomass_row = n_metabolites_metabolic + length(exchange_met_ids)
        for i = 1:length(models)
            biomass_ind = first(indexin([biomass_ids[i]], reactions(models[i])))
            I[kstart] = i + n_before_biomass_row
            J[kstart] = biomass_ind + reaction_offset[i]
            V[kstart] = 1.0
            kstart += 1
        end
    end

    S = sparse(
        I[1:kstart-1],
        J[1:kstart-1],
        V[1:kstart-1],
        n_metabolites_total,
        n_reactions_total,
    ) # could be that some microbes don't have all the exchanges, hence kstart-1

    _reaction_offsets = cumsum(reaction_lengths)
    _metabolite_offsets = cumsum(metabolite_lengths)
    for i = 1:length(models)
        species = isempty(model_names) ? "species_$(i)" : model_names[i]
        tlbs, tubs = bounds(models[i])
        lbs[reaction_offset[i]+1:_reaction_offsets[i]] .= tlbs
        ubs[reaction_offset[i]+1:_reaction_offsets[i]] .= tubs
        rxns[reaction_offset[i]+1:_reaction_offsets[i]] =
            "$(species)_" .* reactions(models[i])
        mets[metabolite_offset[i]+1:_metabolite_offsets[i]] =
            "$(species)_" .* metabolites(models[i])
    end
    mets[_metabolite_offsets[end]+1:_metabolite_offsets[end]+length(exchange_met_ids)] .=
        exchange_met_ids
    rxns[_reaction_offsets[end]+1:_reaction_offsets[end]+length(exchange_rxn_ids)] .=
        exchange_rxn_ids

    if add_biomass_objective
        rxns[end] = "community_biomass"
        for i = 1:length(models)
            species = isempty(model_names) ? "species_$(i)" : model_names[i]
            mets[end-length(biomass_ids)+i] = "$(species)_" .* biomass_ids[i]
        end
    end

    return CoreModel(S, spzeros(size(S, 1)), spzeros(size(S, 2)), lbs, ubs, rxns, mets)
end