text:
import Base: @deprecate

mutable struct Bea
    url::AbstractString
    key::AbstractString
end

function Bea()
    Base.depwarn("The `Bea` struct is deprected; " *
        "`key` is now autoloaded on package startup " *
        "or use `USER_ID = key`.", :Bea)
    key = ""
    if "BEA_KEY" in keys(ENV)
        key = ENV["BEA_KEY"]
    elseif isfile(joinpath(homedir(), ".beadatarc"))
        open(joinpath(homedir(),".beadatarc"), "r") do f
            key = read(f, String)
        end
        key = rstrip(key)
    else
        error("No API key found, connection not initialized")
    end
    println("API key loaded.")

    return Bea(API_URL, key)
end

@deprecate Bea(key) USER_ID = key

function bea_query_dep(url, querydict)
    r = try
        HTTP.request("GET", url, query = querydict)
    catch e
        if isa(e, Base.IOError)
            println("Caught $e")
            HTTP.request("GET", url, query = querydict)
        else
            rethrow(e)
        end
    end
    r_body = String(r.body)
    r_json = JSON.parse(r_body)
    if haskey(r_json["BEAAPI"]["Results"], "Error")
        error(r_json["BEAAPI"]["Results"]["Error"]["APIErrorDescription"])
    end
    return r_json
end

function get_bea_datasets(b::Bea)
    Base.depwarn("`get_bea_datasets` is deprecated, use `bea_datasets` instead.",
        :get_bea_datasets)
    url = b.url
    key = b.key
    bea_method = "GetDataSetList"

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "ResultFormat" => "JSON")

    response_json = bea_query_dep(url, querydict)
    response_dict = response_json["BEAAPI"]["Results"]["Dataset"]
    datasets = Dict(d["DatasetName"] => d["DatasetDescription"] for d in response_dict)

    return datasets
end

function get_bea_parameterlist(b::Bea, dataset::String)
    Base.depwarn("`get_bea_parameterlist` is deprecated, use `bea_parameterlist` instead.",
        :get_bea_parameterlist)
    url = b.url
    key = b.key
    bea_method = "GetParameterList"

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "DatasetName" => dataset,
                     "ResultFormat" => "JSON")

    response_json = bea_query_dep(url, querydict)
    response_dict = response_json["BEAAPI"]["Results"]["Parameter"]

    return response_dict
end

function get_nipa_table(b::Bea, TableName::AbstractString, frequency::AbstractString, startyear::Int, endyear::Int)
    Base.depwarn("`get_nipa_table` is deprecated, use `bea_table` instead.",
        :bea_table)
    url = b.url
    key = b.key
    dataset = "NIPA"
    bea_method = "GetData"
    years = join(collect(startyear:1:endyear), ",")

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "DatasetName" => dataset,
                     "TableName" => TableName,
                     "Frequency" => frequency,
                     "Year" => years,
                     "ResultFormat" => "JSON")

    response_json = bea_query(url, querydict)

    # Extract the vector of Dicts containing the data
    # data_dicts = response_json["BEAAPI"]["Results"]["Data"]
    data_dicts = response_json["Data"] # new format
    # Retrieve data from each Dict
    all_data = [parse_data_dict_dep(d) for d in data_dicts]
    values = [tup[1] for tup in all_data]
    dates = [tup[2] for tup in all_data]
    linenums = [tup[3] for tup in all_data]
    # Create data frame of values, reshape, and create names
    dflong = DataFrame(date = dates, line = linenums, value = values)
    df = unstack(dflong, :date, :line, :value) # Convert to "wide" format
    newnames = [Symbol(string("line", name)) for name in names(df)[2:end]]
    rename!(df, [:date; newnames])

    # Create OrderedDict of line descriptions
    linedesc = [tup[4] for tup in all_data]
    linekeys = OrderedDict{AbstractString, AbstractString}(zip(linenums, linedesc))

    # Extract metadata for the table
    # tablenotes = response_json["BEAAPI"]["Results"]["Notes"]
    tablenotes = response_json["Notes"]
    tableinfo = split(tablenotes[1]["NoteText"])
    tablenum = rstrip(tableinfo[2], '.')
    tabledesc = join(tableinfo[3:end], " ")
    if length(tablenotes) > 1
        noteref = [note["NoteRef"] for note in tablenotes[2:end]]
        notetext = [note["NoteText"] for note in tablenotes[2:end]]
        notes = OrderedDict{AbstractString, AbstractString}(zip(noteref, notetext))
    else
        notes = ""
    end

    frequency == "Q" ? freq = "Quarterly" : freq = "Annual"
    out = BeaTable("NIPA", tablenum, tabledesc, "", "", DataFrame(), DataFrame(),
        freq, startyear, endyear, TableName, "", df)
    return out
end

function parse_data_dict_dep(dict::Dict)
    linenum = dict["LineNumber"]
    linedesc = dict["LineDescription"]
    # Add footnote indicator to line description
    if haskey(dict, "NoteRef")
        fn = dict["NoteRef"]
        linedesc = string(linedesc, " (", fn, ")")
    end
    # Remove commas from data values
    dataval = parse(Float64, replace(dict["DataValue"], "," => ""))

    # Change date from string to Date() type
    timeperiod = dict["TimePeriod"]
    year = parse(Int, timeperiod[1:4])
    # Quarterly data
    if occursin(r"Q", timeperiod)
        quarter = parse(Int, timeperiod[end])
        date = Date(year, (quarter*3), 1)
    # Annual data
    else
        date = Date(year, 12, 1)
    end

    return (dataval, date, linenum, linedesc)
end


text:
using Parameters: @with_kw, @unpack
using LinearAlgebra: eigvals
using ForwardDiff
using QuadGK: quadgk

# Omnivory preference functions
function adapt_pref(u, p, t)
    return p.ω * u[1] / (p.ω * u[1] + (1 - p.ω) * u[2])
end

function fixed_pref(u, p, t)
    return p.Ω
end

# Parameters need to be defined after above functions so they can have defaults
@with_kw mutable struct ModelPar
    # Logistic Parameters
    r = 2.0
    # `K_base` measures the underyling K outside of any forcing applied
    K_base = 3
    K = 3
    # Consumer Parameters
    a_RC = 1
    h_RC = 0.4
    e_RC = 0.8
    m_C = 0.4
    # Predator Parameters
    a_CP = 0.5
    h_CP = 0.3
    e_CP = 0.6
    m_P = 0.2
    # Omnivory Parameters
    a_RP = 0.2
    h_RP = 0.6
    e_RP = 0.4
    # if Ω = 0 then we have a food chain
    Ω = 0.1
    # Forcing Function
    # 'fixed_pref' for passive omnivory
    # 'adapt_pref' for active omnivory
    pref::Function = fixed_pref
    # Used in the adaptive forcing to bias towards C or R
    ω = 0.5
end

#NOTE: we could use these in the model, but I am scared of all the function call
#      overhead, likely could be fixed with inlining, but will just leave it for now
function f_RP(u, p, t)
    R, C, P = u
    Ω = p.pref(u, p, t)
    return Ω * p.a_RP * R * P / (1 + Ω * p.a_RP * p.h_RP * R + (1 - Ω) * p.a_CP * p.h_CP * C)
end

function f_CP(u, p, t)
    R, C, P = u
    Ω = p.pref(u, p, t)
    return (1 - Ω) * p.a_CP * C * P / (1 + Ω * p.a_RP * p.h_RP * R + (1 - Ω) * p.a_CP * p.h_CP * C)
end

degree_omnivory(u, p) = f_RP(u, p, 0.0) / (f_RP(u, p, 0.0) + f_CP(u, p, 0.0))

function model!(du, u, p, t)
    @unpack r, K = p
    @unpack a_RC, h_RC, e_RC, m_C = p
    @unpack a_CP, h_CP, e_CP, m_P = p
    @unpack a_RP, h_RP, e_RP = p
    R, C, P = u

    # setup the density dependent preference
    Ω = p.pref(u, p, t)

    # compute those once for all
    int_RC = a_RC * R * C / (1 + a_RC * h_RC * R)
    denom_RCP = 1 + Ω * a_RP * h_RP * R + (1 - Ω) * a_CP * h_CP * C
    num_CP = (1 - Ω) * a_CP * C * P
    num_RP = Ω * a_RP * R * P

    # ODE
    du[1] = r * R * (1 - R / K) - int_RC - num_RP / denom_RCP
    du[2] = e_RC * int_RC - num_CP / denom_RCP - m_C * C
    du[3] = (e_RP * num_RP + e_CP * num_CP) / denom_RCP - m_P * P

    return du
end

# Utilities for doing eigenvalue analysis
function rhs(u, p)
    du = similar(u)
    model!(du, u, p, zero(u))
    return du
end

find_eq(u, p) = nlsolve((du, u) -> model!(du, u, p, zero(u)), u).zero
cmat(u, p) = ForwardDiff.jacobian(x -> rhs(x, p), u)

"""M is the community matrix, we can be calculated with `cmat(u, p)`"""
λ1_stability(M) = maximum(real.(eigvals(M)))

"""M is the community matrix, we can be calculated with `cmat(u, p)`"""
function λ1_stability_imag(M)
    ev = eigvals(M)
    imag.(ev)[findall(real.(ev) .== maximum(real.(ev)))[1]]
end

"""M is the community matrix, we can be calculated with `cmat(u, p)`
Note: `\nu` is the what to input `ν` which looks a bit too much like `v` for my taste
"""
ν_stability(M) = λ1_stability((M + M') / 2)

# overshoot and oscillation range

abs_sol(sol, t, eq) = abs.(sol(t) .- eq)

function overshoot(sol, eq, spc, t_beg, t_end)
    return quadgk(t -> abs_sol(sol, t, eq)[spc], t_beg, t_end)[1]
end

# Calculate max-min metric
function min_max(sol, spc, t_beg, t_end, len = 100000)
    return maximum(sol(range(t_beg, t_end, length = len))[spc, :]) -
    minimum(sol(range(t_beg, t_end, length = len))[spc, :])
end

# find first time equilibrium is hit
function find_times_hit_equil_press(res)
    eq = res[1, end], res[2, end], res[3, end]
    times = zeros(3)
    for spc in 1:3
        for i in 20:length(res)
            # cannot be too strict here otherwise the value of the
            # first ht time varies a lort which will have serious
            # impact on min and max (overshoot too) leading to major oscillations lenth of the ts must be high enough too.
            if isapprox(res[spc, i], eq[spc], atol = 0.01)
                times[spc] = res.t[i]
                break
            end
        end
    end
    return times
end

# Basic plot for Sensitivity analysis (well akin to sa)
function plot_sa_unit(rg, res, id, leg = false, xlb = "", ylb = "", ann = 1)
    # az = 'a':'z'
    cols = ["#000000" "#555555" "#cccccc"]
    labs = ["FC" "PO" "RO"]
    lty = ["solid" "dashed" "solid"]
    for j in 1:3
        plot(rg, [res[i][j][id] for i in eachindex(rg)],
            label = labs[j], color = cols[j], linestyle = lty[j])
    end
    if leg
        legend()
    end
    xlabel(xlb)
    ylabel(ylb)
    # did not work
    # annotate(string(az[ann]), [.1;.1], xycoords="figure fraction",
    #     annotation_clip = false)
end

# Basic plot for Sensitivity analysis (well akin to sa)
function plot_sa_unit_rev(rg, res, id, leg = false, xlb = "", ylb = "")
    # az = 'a':'z'
    cols = ["#000000" "#555555" "#cccccc"]
    labs = ["FC" "PO" "RO"]
    lty = ["solid" "dashed" "solid"]
    for j in 1:3
        plot(rg, [res[i][j][id] for i in reverse(eachindex(rg))], label = labs[j], color = cols[j], linestyle = lty[j])
    end
    if leg
        legend()
    end
    xlabel(xlb)
    ylabel(ylb)
end

# Helper function to check out degree of omnivory
function check_doo(res, rgs)
    type = ["Fixed" "Responsive"]
    phase = ["Equilibrium", "Pulse", "Transient", "New Equilibrium"]
    # Check values of degree of omnivory
    for i in 2:3
        printstyled(type[i - 1], "---------\n", color = :green)
        for j in eachindex(rgs)
            printstyled(phase[j], " --> ", color = :blue)
            # see press_unit output
            sol = res[i][14](rgs[j])
            println(
                maximum(
                    [degree_omnivory(
                    sol[:, k],
                    res[i][13]
                    ) for k in eachindex(sol)]
                )
            )
            # maximum degree of omnivory when the C is min (food web is top heavy)
            printstyled("where topheaviest --> ", color = :blue)
            idmx = argmin([sol[2, k] for k in eachindex(sol)])
            # println(idmx)
            print(
                degree_omnivory(
                    sol[idmx],
                    res[i][13]
                )
            )
            ti = rgs[j][idmx]
            println("(time = $ti)")
        end
    end
end

ext: jl

max_stars_repo_name: McCannLab/Labmnivory

max_issues_repo_name: McCannLab/Labmnivory

max_forks_repo_name: McCannLab/Labmnivory


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Union arrays
# need a custom representation of Union types since arrow unions
# are ordered, and possibly indirected via separate typeIds array
# here, T is Meta.UnionMode.Dense or Meta.UnionMode.Sparse,
# typeIds is a NTuple{N, Int32}, and U is a Tuple{...} of the
# unioned types
struct UnionT{T, typeIds, U}
end

unionmode(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = T
typeids(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = typeIds
Base.eltype(::Type{UnionT{T, typeIds, U}}) where {T, typeIds, U} = U

ArrowTypes.ArrowType(::Type{<:UnionT}) = ArrowTypes.UnionType()

# iterate a Julia Union{...} type, producing an array of unioned types
function eachunion(U::Union, elems=nothing)
    if elems === nothing
        return eachunion(U.b, Type[U.a])
    else
        push!(elems, U.a)
        return eachunion(U.b, elems)
    end
end

function eachunion(T, elems)
    push!(elems, T)
    return elems
end

# produce typeIds, offsets, data tuple for DenseUnion
isatypeid(x::T, ::Type{types}) where {T, types} = isatypeid(x, fieldtype(types, 1), types, 1)
isatypeid(x::T, ::Type{S}, ::Type{types}, i) where {T, S, types} = x isa S ? i : isatypeid(x, fieldtype(types, i + 1), types, i + 1)

"""
    Arrow.DenseUnion

An `ArrowVector` where the type of each element is one of a fixed set of types, meaning its eltype is like a julia `Union{type1, type2, ...}`.
An `Arrow.DenseUnion`, in comparison to `Arrow.SparseUnion`, stores elements in a set of arrays, one array per possible type, and an "offsets"
array, where each offset element is the index into one of the typed arrays. This allows a sort of "compression", where no extra space is
used/allocated to store all the elements.
"""
struct DenseUnion{T, S} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    arrow2::Vector{UInt8} # if arrow blob is compressed, need a 2nd reference for uncompressed offsets bytes
    typeIds::Vector{UInt8}
    offsets::Vector{Int32}
    data::S # Tuple of ArrowVector
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::DenseUnion) = size(s.typeIds)
nullcount(x::DenseUnion) = 0 # DenseUnion has no validity bitmap; only children do

@propagate_inbounds function Base.getindex(s::DenseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    @inbounds off = s.offsets[i]
    @inbounds x = s.data[typeId + 1][off + 1]
    return x
end

@propagate_inbounds function Base.setindex!(s::DenseUnion{UnionT{T, typeIds, U}}, v, i::Integer) where {T, typeIds, U}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    vtypeId = Int8(typeids[isatypeid(v, U)])
    if typeId == vtypeId
        @inbounds off = s.offsets[i]
        @inbounds s.data[typeId +1][off + 1] = v
    else
        throw(ArgumentError("type of item to set $(typeof(v)) must match existing item $(fieldtype(U, typeid))"))
    end
    return v
end

# convenience wrappers for signaling that an array shoudld be written
# as with dense/sparse union arrow buffers
struct DenseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Dense, nothing, U}}
    itr::T
end

DenseUnionVector(x::T) where {T} = DenseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:DenseUnionVector}) = Base.IndexLinear()
Base.size(x::DenseUnionVector) = (length(x.itr),)
Base.iterate(x::DenseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::DenseUnionVector, i::Int) = getindex(x.itr, i)

function todense(::Type{UnionT{T, typeIds, U}}, x) where {T, typeIds, U}
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    len = length(x)
    types = Vector{UInt8}(undef, len)
    offsets = Vector{Int32}(undef, len)
    data = Tuple(Vector{i == 1 ? Union{Missing, fieldtype(U, i)} : fieldtype(U, i)}(undef, 0) for i = 1:fieldcount(U))
    for (i, y) in enumerate(x)
        typeid = y === missing ? 0x00 : UInt8(typeids[isatypeid(y, U)])
        @inbounds types[i] = typeid
        @inbounds offsets[i] = length(data[typeid + 1])
        push!(data[typeid + 1], y)
    end
    return types, offsets, data
end

struct SparseUnionVector{T, U} <: AbstractVector{UnionT{Meta.UnionMode.Sparse, nothing, U}}
    itr::T
end

SparseUnionVector(x::T) where {T} = SparseUnionVector{T, Tuple{eachunion(eltype(x))...}}(x)
Base.IndexStyle(::Type{<:SparseUnionVector}) = Base.IndexLinear()
Base.size(x::SparseUnionVector) = (length(x.itr),)
Base.iterate(x::SparseUnionVector, st...) = iterate(x.itr, st...)
Base.getindex(x::SparseUnionVector, i::Int) = getindex(x.itr, i)

# sparse union child array producer
# for sparse unions, we split the parent array into
# N children arrays, each having the same length as the parent
# but with one child array per unioned type; each child
# should include the elements from parent of its type
# and other elements can be missing/default
function sparsetypeids(::Type{UnionT{T, typeIds, U}}, x) where {T, typeIds, U}
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    len = length(x)
    types = Vector{UInt8}(undef, len)
    for (i, y) in enumerate(x)
        typeid = y === missing ? 0x00 : UInt8(typeids[isatypeid(y, U)])
        @inbounds types[i] = typeid
    end
    return types
end

struct ToSparseUnion{T, A} <: AbstractVector{T}
    data::A
end

ToSparseUnion(::Type{T}, data::A) where {T, A} = ToSparseUnion{T, A}(data)

Base.IndexStyle(::Type{<:ToSparseUnion}) = Base.IndexLinear()
Base.size(x::ToSparseUnion) = (length(x.data),)

Base.@propagate_inbounds function Base.getindex(A::ToSparseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(A, i)
    @inbounds x = A.data[i]
    return @inbounds x isa T ? x : ArrowTypes.default(T)
end

function compress(Z::Meta.CompressionType, comp, x::A) where {A <: DenseUnion}
    len = length(x)
    nc = nullcount(x)
    typeIds = compress(Z, comp, x.typeIds)
    offsets = compress(Z, comp, x.offsets)
    buffers = [typeIds, offsets]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, comp, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

"""
    Arrow.SparseUnion

An `ArrowVector` where the type of each element is one of a fixed set of types, meaning its eltype is like a julia `Union{type1, type2, ...}`.
An `Arrow.SparseUnion`, in comparison to `Arrow.DenseUnion`, stores elements in a set of arrays, one array per possible type, and each typed
array has the same length as the full array. This ends up with "wasted" space, since only one slot among the typed arrays is valid per full
array element, but can allow for certain optimizations when each typed array has the same length.
"""
struct SparseUnion{T, S} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    typeIds::Vector{UInt8}
    data::S # Tuple of ArrowVector
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(s::SparseUnion) = size(s.typeIds)
nullcount(x::SparseUnion) = 0

@propagate_inbounds function Base.getindex(s::SparseUnion{T}, i::Integer) where {T}
    @boundscheck checkbounds(s, i)
    @inbounds typeId = s.typeIds[i]
    @inbounds x = s.data[typeId + 1][i]
    return x
end

@propagate_inbounds function Base.setindex!(s::SparseUnion{UnionT{T, typeIds, U}}, v, i::Integer) where {T, typeIds, U}
    @boundscheck checkbounds(s, i)
    typeids = typeIds === nothing ? (0:(fieldcount(U) - 1)) : typeIds
    vtypeId = Int8(typeids[isatypeid(v, U)])
    @inbounds s.typeIds[i] = vtypeId
    @inbounds s.data[vtypeId + 1][i] = v
    return v
end

arrowvector(U::Union, x, i, nl, fi, de, ded, meta; denseunions::Bool=true, kw...) =
    arrowvector(denseunions ? DenseUnionVector(x) : SparseUnionVector(x), i, nl, fi, de, ded, meta; denseunions=denseunions, kw...)

arrowvector(::UnionType, x::Union{DenseUnion, SparseUnion}, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::UnionType, x, i, nl, fi, de, ded, meta; kw...)
    UT = eltype(x)
    if unionmode(UT) == Meta.UnionMode.Dense
        x = x isa DenseUnionVector ? x.itr : x
        typeids, offsets, data = todense(UT, x)
        data2 = map(y -> arrowvector(y[2], i, nl + 1, y[1], de, ded, nothing; kw...), enumerate(data))
        return DenseUnion{UT, typeof(data2)}(UInt8[], UInt8[], typeids, offsets, data2, meta)
    else
        x = x isa SparseUnionVector ? x.itr : x
        typeids = sparsetypeids(UT, x)
        data3 = Tuple(arrowvector(ToSparseUnion(fieldtype(eltype(UT), j), x), i, nl + 1, j, de, ded, nothing; kw...) for j = 1:fieldcount(eltype(UT)))
        return SparseUnion{UT, typeof(data3)}(UInt8[], typeids, data3, meta)
    end
end

function compress(Z::Meta.CompressionType, comp, x::A) where {A <: SparseUnion}
    len = length(x)
    nc = nullcount(x)
    typeIds = compress(Z, comp, x.typeIds)
    buffers = [typeIds]
    children = Compressed[]
    for y in x.data
        push!(children, compress(Z, comp, y))
    end
    return Compressed{Z, A}(x, buffers, len, nc, children)
end

function makenodesbuffers!(col::Union{DenseUnion, SparseUnion}, fieldnodes, fieldbuffers, bufferoffset, alignment)
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # typeIds buffer
    push!(fieldbuffers, Buffer(bufferoffset, len))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    bufferoffset += padding(len, alignment)
    if col isa DenseUnion
        # offsets buffer
        blen = sizeof(Int32) * len
        push!(fieldbuffers, Buffer(bufferoffset, blen))
        @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
        bufferoffset += padding(blen, alignment)
    end
    for child in col.data
        bufferoffset = makenodesbuffers!(child, fieldnodes, fieldbuffers, bufferoffset, alignment)
    end
    return bufferoffset
end

function writebuffer(io, col::Union{DenseUnion, SparseUnion}, alignment)
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    # typeIds buffer
    n = writearray(io, UInt8, col.typeIds)
    @debug 1 "writing array: col = $(typeof(col.typeIds)), n = $n, padded = $(padding(n, alignment))"
    writezeros(io, paddinglength(n, alignment))
    if col isa DenseUnion
        n = writearray(io, Int32, col.offsets)
        @debug 1 "writing array: col = $(typeof(col.offsets)), n = $n, padded = $(padding(n, alignment))"
        writezeros(io, paddinglength(n, alignment))
    end
    for child in col.data
        writebuffer(io, child, alignment)
    end
    return
end

max_stars_repo_name: timkpaine/arrow

max_issues_repo_name: timkpaine/arrow

max_forks_repo_name: XpressAI/arrow


text:
using StatisticalRethinking

println( "All models are wrong, but some are useful." );

x = 1:3

x = x*10

x = log.(x)

x = sum(x)

x = exp(x)
x = x*10
x = log(x)
x = sum(x)
x = exp(x)

[log(0.01^200) 200 * log(0.01)]

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# Script for creating Figure 4 in scale manuscript

using IcanProj
using CSV
using DataFrames
using NearestNeighbors
using NetCDF
using Dates
using Statistics
using PyPlot
using JFSM2

# Function for computing error statistics

function compute_error(path_res, cfg, df_obs, idxs, time_sim)

    file_swe = joinpath(path_res, "results_$(cfg)", "swe_1km.nc")
    file_hs = joinpath(path_res, "results_$(cfg)", "snowdepth_1km.nc")

    swe = ncread(file_swe, "swe")
    hs = ncread(file_hs, "snowdepth")

    swe_obs, swe_sim = Float64[], Float64[]
    hs_obs, hs_sim = Float64[], Float64[]

    for ix in 1:length(idxs)

        idx = idxs[ix]

        time_obs = DateTime(df_obs.year[ix], df_obs.month[ix], df_obs.day[ix])

        it = findfirst(t -> t == time_obs, time_sim)

        if (it != nothing)

            push!(swe_obs, 1000*df_obs.swe[ix])
            push!(swe_sim, swe[idx, it])

            push!(hs_obs, df_obs.hs[ix])
            push!(hs_sim, hs[idx, it])

        end

    end

    @info "Number of measurements: $(length(swe_obs))"

    bias_swe = mean(swe_sim .- swe_obs)
    rmse_swe = sqrt(mean((swe_sim .- swe_obs).^2))
    mean_swe_obs = mean(swe_obs);

    bias_hs = mean(hs_sim .- hs_obs)
    rmse_hs = sqrt(mean((hs_sim .- hs_obs).^2))
    mean_hs_obs = mean(hs_obs);

    return bias_swe, rmse_swe, mean_swe_obs, bias_hs, rmse_hs, mean_hs_obs

end

# Settings

path_res = "/data04/jmg/fsm_simulations/netcdf/fsmres_open"

path_figure = joinpath(dirname(pathof(IcanProj)), "..", "plots", "snow_evaluation")

# Read tables

df_links = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "df_links.csv")) |> DataFrame

df_obs = CSV.File(joinpath(dirname(pathof(IcanProj)), "..", "data", "snowdata.csv")) |> DataFrame

# Load metadata from simulations

file = joinpath(path_res, "results_32/swe_1km.nc")

time_sim = ncread(file, "time_array")

time_sim = [DateTime(time_sim[1, i], time_sim[2, i], time_sim[3, i], time_sim[4, i]) for i in 1:size(time_sim, 2)]

df_tmp = DataFrame(ind_senorge = convert.(Int, ncread(file, "id")))

df_links = join(df_tmp, df_links, on = :ind_senorge, kind = :left)

# Search for indicies for the model

coord_model = permutedims([df_links.xcoord df_links.ycoord])

coord_obs = permutedims([df_obs.xcoord df_obs.ycoord])

coord_model[ismissing.(coord_model)] .= -1000000

tree = BruteTree(coord_model)

idxs, dists = knn(tree, coord_obs, 1, true)

idxs = [x[1] for x in idxs]

dists = [x[1] for x in dists]

# Compute error statistics for all configurations

bias_swe_all, rmse_swe_all, mean_swe_obs_all = [], [], []
bias_hs_all, rmse_hs_all, mean_hs_obs_all = [], [], [], []

for cfg in 1:32

    println(cfg)

    bias_swe, rmse_swe, mean_swe_obs, bias_hs, rmse_hs, mean_hs_obs = compute_error(path_res, cfg, df_obs, idxs, time_sim)

    push!(bias_swe_all, bias_swe)
    push!(rmse_swe_all, rmse_swe)
    push!(mean_swe_obs_all, mean_swe_obs)

    push!(bias_hs_all, bias_hs)
    push!(rmse_hs_all, rmse_hs)
    push!(mean_hs_obs_all, mean_hs_obs)

end

# Save to tables

df_res = DataFrame(bias_swe_all = bias_swe_all,
                   rmse_swe_all = rmse_swe_all,
                   mean_swe_obs_all = mean_swe_obs_all,
                   bias_hs_all = bias_hs_all,
                   rmse_hs_all = rmse_hs_all,
                   mean_hs_obs_all = mean_hs_obs_all)

df_res |> CSV.write(joinpath(path_figure, "table_snow_evaluation.csv"))

# Plot results

ax = figure(figsize = (10, 4))
plot(bias_swe_all ./ mean_swe_obs_all, marker="o", linestyle="None")
xlabel("Model configuration")
ylabel("Normalized bias (-)")
title("SWE")
xticks(collect(0:31), collect(1:32))
savefig(joinpath(path_figure, "bias_swe.png"), dpi = 200)
close()

ax = figure(figsize = (10, 4))
plot(rmse_swe_all ./ mean_swe_obs_all, marker="o", linestyle="None")
xlabel("Model configuration")
ylabel("Normalized RMSE (-)")
title("SWE")
xticks(collect(0:31), collect(1:32))
savefig(joinpath(path_figure, "rmse_swe.png"), dpi = 200)
close()

ax = figure(figsize = (10, 4))
plot(bias_hs_all ./ mean_hs_obs_all, marker="o", linestyle="None")
xlabel("Model configuration")
ylabel("Normalized bias (-)")
title("Snow depth")
xticks(collect(0:31), collect(1:32))
savefig(joinpath(path_figure, "bias_hs.png"), dpi = 200)
close()

ax = figure(figsize = (10, 4))
plot(rmse_hs_all ./ mean_hs_obs_all, marker="o", linestyle="None")
xlabel("Model configuration")
ylabel("Normalized RMSE (-)")
title("Snow depth")
xticks(collect(0:31), collect(1:32))
savefig(joinpath(path_figure, "rmse_hs.png"), dpi = 200)
close()

# Box plots

cfgs = cfg_table()

fig, axes = plt[:subplots](nrows = 3, ncols = 2)

fig[:set_size_inches](8, 6)

fig[:text](0.04, 0.5, "Bias (%)", va="center", rotation="vertical")

fig[:text](1-0.04, 0.5, "NRMSE (-)", va="center", rotation="vertical")

fig[:text](0.5, 0.04, "Parametrization option", ha = "center")

# Bias for snow water equivalent

variable = 100 .* df_res.bias_swe_all ./ df_res.mean_swe_obs_all

explanation = cfgs.exchng

axes[1][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[1][:xaxis][:set_ticklabels]([])

axes[1][:annotate]("Exchng", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

explanation = cfgs.hydrol

axes[2][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[2][:xaxis][:set_ticklabels]([])

axes[2][:annotate]("Hydraul", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

explanation = cfgs.albedo

axes[3][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[3][:set_xticklabels](["0", "1"])

axes[3][:annotate]("Albedo", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

# RMSE for snow water equivalent

variable = df_res.rmse_swe_all ./ df_res.mean_swe_obs_all

explanation = cfgs.exchng

axes[4][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[4][:xaxis][:set_ticklabels]([])

axes[4][:yaxis][:tick_right]()

axes[4][:yaxis][:set_label_position]("right")

axes[4][:annotate]("Exchng", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

explanation = cfgs.hydrol

axes[5][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[5][:xaxis][:set_ticklabels]([])

axes[5][:yaxis][:tick_right]()

axes[5][:yaxis][:set_label_position]("right")

axes[5][:annotate]("Hydraul", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

explanation = cfgs.albedo

axes[6][:boxplot]([variable[explanation .== 0] variable[explanation .== 1]])

axes[6][:set_xticklabels](["0", "1"])

axes[6][:yaxis][:tick_right]()

axes[6][:yaxis][:set_label_position]("right")

axes[6][:annotate]("Albedo", xy=[0.03; 0.8], xycoords="axes fraction", fontsize=10.0)

# Save figure

savefig(joinpath(path_figure, "error_vs_parameterization.pdf"), dpi = 600)


text:
@testset "arraymath" for i = 1:5

x = CuArray(T,10,5)
log(x)
exp(x)

end


text:
struct FieldNameAndValuePair
    name::String
    value
end

function AbstractTrees.printnode(io::IO, x::FieldNameAndValuePair)
    if x.value isa Ptr && x.value == C_NULL
        Base.print(io, string(x.name, " -> ( 0x0 )"))
    else
        Base.print(io, x.name)
    end
end

function AbstractTrees.children(x::FieldNameAndValuePair)
    x.value isa Ptr && x.value == C_NULL && return ()
    return (simple_type_value(x.value), )
end

node_tag_enum_names = map(x -> x[1], CEnum.name_value_pairs(NodeTag))

const IGNORE_LIST = [
    :T_Invalid,
    :T_TIDBitmap, #Cvoid
    :T_List,
    :T_SeqScan, #duplicate of Scan
    :T_DistinctExpr, #duplicate of OpExpr
    :T_NullIfExpr, #duplicate of OpExpr
    :T_Value, # special union struct
    :T_String,
    :T_Null,
    :T_Integer,
    :T_Float,
]

target_node_tag_enum_names = filter(x -> !(x in IGNORE_LIST), node_tag_enum_names)

for enum_name in target_node_tag_enum_names
    enum_name_str = string(enum_name)
    @assert startswith(enum_name_str, "T_")
    tp_name_str = SubString(enum_name_str, 3)
    tp_name = Symbol(tp_name_str)
    try
        tp = eval(Meta.parse("PGQuery.$tp_name"))
        tp == Nothing && continue
        local printnode_body = quote
            hasfield($tp, :type) && $(enum_name != :T_Expr) && @assert x.type == $enum_name "x.type is $(x.type), but expected to be $($enum_name)"
            Base.print(io, $tp_name_str)
            Base.print(io, "(")
        end
        local children_body = quote
            hasfield($tp, :type) && $(enum_name != :T_Expr) && @assert x.type == $enum_name "x.type is $(x.type), but expected to be $($enum_name)"
            children = []
        end
        is_first_simple_field = true
        for (i, (fld, fld_tp)) in enumerate(zip(fieldnames(tp), tp.types))
            fld == :type && continue
            fld_name = string(fld)
            if is_simple_type(fld_tp)
                printnode_body = quote
                    $printnode_body
                    $(!is_first_simple_field) && Base.print(io, ", ")
                    Base.print(io, $fld_name)
                    Base.print(io, "=")
                    Base.print(io, simple_type_value(x.$fld))
                end
                is_first_simple_field = false
            else
                children_body = quote
                    $children_body
                    push!(children, FieldNameAndValuePair($fld_name, x.$fld))
                end
            end
        end

        q = (quote
            function convert_to_node_type(::Val{$enum_name})
                return $tp_name
            end
            function AbstractTrees.printnode(io::IO, x::$tp_name)
                $printnode_body
                Base.print(io, ")")
            end

            function AbstractTrees.children(x::$tp_name)
                $children_body
                return children
            end
        end)
        eval(q)
    catch e
        #swallow the error
    end
end

function convert_to_node_type(::Val{T}) where T
    str_name = string(T)
    @assert startswith(str_name, "T_")
    type_name = SubString(str_name, 3)
    return eval(Meta.parse("PGQuery.$type_name"))
end

function convert_to_node_type(::Val{T_List})
    return List
end

function AbstractTrees.printnode(io::IO, x::List)
    @assert x.type == T_List
    Base.print(io, string("List(length=",x.length,")"))
end

function AbstractTrees.children(x::Ptr)
    return ()
end

function AbstractTrees.printnode(io::IO, x::Value)
    actual_value = convert_to_actual_value(x)
    if actual_value === nothing
        Base.print(io, "null")
    else
        AbstractTrees.printnode(io, actual_value === nothing ? C_NULL : actual_value)
    end
end


text:
import Entropies: SymbolicPermutation, CountOccurrences, symbolize, symbolize!
export SymbolicPermutation

function transferentropy(s, t, est::SymbolicPermutation;
        base = 2, q = 1, τT = -1, τS = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1)

    # Generalized embeddings of each marginal separately, using the same dimension for each
    # marginal. Then symbolize the embedded time series, which gives integer symbol sequences.
    symb_s = symbolize(s, est)
    symb_t = symbolize(t, est)

    # After symbolization, use the general interface with occurrence frequency entropy
    # estimator (probabilities are obtained from histograms of the symbols).
    transferentropy(symb_s, symb_t, CountOccurrences(),
        τT = τT, τS = τS, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯)
end

function transferentropy!(symb_s, symb_t, s, t, est::SymbolicPermutation; base = 2, q = 1,
        τT = -1, τS = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1)

    symbolize!(symb_s, s, est)
    symbolize!(symb_t, t, est)

    transferentropy(symb_s, symb_t, CountOccurrences(),
        τT = τT, τS = τS, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯)
end

function transferentropy(s, t, c, est::SymbolicPermutation; base = 2, q = 1,
        τT = -1, τS = -1, τC = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1, dC = 1)

    symb_s = symbolize(s, est)
    symb_t = symbolize(t, est)
    symb_c = symbolize(c, est)

    transferentropy(symb_s, symb_t, symb_c, CountOccurrences(),
        τT = τT, τS = τS, τC = τC, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯, dC = dC)
end

function transferentropy!(symb_s, symb_t, symb_c, s, t, c, est::SymbolicPermutation;
        base = 2, q = 1,
        τT = -1, τS = -1, τC = -1, η𝒯 = 1, dT = 1, dS = 1, d𝒯 = 1, dC = 1)

    symbolize!(symb_s, s, est)
    symbolize!(symb_t, t, est)
    symbolize!(symb_c, c, est)

    transferentropy(symb_s, symb_t, symb_c, CountOccurrences(),
        τT = τT, τS = τS, τC = τC, η𝒯 = η𝒯, dT = dT, dS = dS, d𝒯 = d𝒯, dC = dC)
end


text:
module Binary_Cross_Entropy_Loss
    using LoopVectorization

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss_matrix = zeros(Float32, size(output_matrix))
        @avxt for i in axes(loss_matrix, 1), j in axes(loss_matrix, 2)
            loss_matrix[i,j] = -(sample*log(max(output_matrix[i,j], 1e-8))+(1-sample_matrix[i,j])*log(max(1-output_matrix[i,j], 1e-8)))
        end
        return loss_matrix
    end

    function prop!(δ::Array{Float32}, output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        @avxt for i in eachindex(δ)
            δ[i] = (output_matrix[i]-sample_matrix[i])/(output_matrix[i]-output_matrix[i]^2+1e-8)
        end
    end
end


text:
import MittagLefflerFunctions.Quadrature.H2
import MittagLefflerFunctions.Quadrature.H3
using PyPlot

α = 1.7
β = 1.2
sep = 0.02
x = 2.0
M = 200
t = range(-0.2, 0.2, length=M)
θ = π/4

w = x^(1/α) .+ t * exp(Complex(0.0, θ))
y = H2.(α, β, w, x, sep)

figure(1)
H0 = ( 1 + α - 2β ) / ( 2α * x^(β/α) )
plot(t, real.(y), t, imag.(y), [0], [H0], "o")
grid(true)

γ₊ = x^(1/α) * exp(Complex(0.0, π/α))
γ₋ = x^(1/α) * exp(Complex(0.0, -π/α))
w = γ₊ .+ t * exp(Complex(0.0, θ))

y1 = zeros(Complex{Float64}, M)
y2 = zeros(Complex{Float64}, M)
y3 = zeros(Complex{Float64}, M)
for m = 1:M
    H₊, H₋ = H3(α, β, w[m], x, sep)
    y1[m] = H₊
    y2[m] = H₋
    y3[m] = ( ( w[m] - γ₊ ) * ( w[m] - γ₋ )
             / ( ( w[m]^α + x ) * ( 2w[m] - γ₊ - γ₋ ) ) )
end

figure(2)
H0₊ = ( ( ( 1 + α - 2β ) * ( γ₊ - γ₋ ) - 2γ₊ )
      / ( 2α * γ₊^β * ( γ₊ - γ₋ ) ) )
plot(t, real.(y1), t, imag.(y1),
     [0], [real(H0₊)], "o", [0], [imag(H0₊)], "o")
grid(true)

figure(3)
H0₋ = ( sinpi((1-β)/α) / ( α * sinpi(1/α) ) ) / x^(β/α)
plot(t, real.(y2), t, imag.(y2), [0], [H0₋], "o")
grid(true)

figure(4)
ϕ0 = γ₊^(1-α) / α
plot(t, real.(y3), t, imag.(y3),
     [0], [ real(ϕ0) ], "o", [0], [ imag(ϕ0) ], "o")
grid(true)


text:
import Base: ==

########## Display ##########

@enum Segment a=Int('a') b=Int('b') c=Int('c') d=Int('d') e=Int('e') f=Int('f') g=Int('g')

mutable struct Display
    signal::Vector{Segment}
    digit::Union{Int, Missing}
end

signal(display::Display) = display.signal
digit(display::Display) = display.digit
decoded(display::Display) = !ismissing(digit(display))

Display() = Display(Vector{Segment}(undef, 0), missing)

function Display(signal_str::String)
    signal = Segment.(Int.(collect(signal_str)))
    n_segments = length(signal)
    n_segments == 2 ? digit = 1 :
    n_segments == 3 ? digit = 7 :
    n_segments == 4 ? digit = 4 :
    n_segments == 7 ? digit = 8 :
    digit = missing
    Display(signal, digit)
end

Base.Int(display::Display) =
    decoded(display) ? digit(display) : error("display not decoded")

function contains(display, sub_display)
    for segment in signal(sub_display)
        !(segment in signal(display)) && return false
    end
    return true
end

display1::Display == display2::Display =
    contains(display1, display2) && contains(display2, display1)

########## Entry ##########

mutable struct Entry
    patterns::Vector{Display}
    output::Vector{Display}
end

patterns(entry::Entry) = entry.patterns
output(entry::Entry) = entry.output

function Entry(entry_str::String)
    patterns_str, output_str = String.(split(entry_str, " | "))
    patterns = Display.(String.(split(patterns_str, " ")))
    output = Display.(String.(split(output_str, " ")))
    Entry(patterns, output)
end

Base.Int(entry::Entry) = parse(Int, join(string.(Int.(output(entry)))))

function get_pattern(entry::Entry, digit::Int)
    for pattern in patterns(entry)
        if decoded(pattern)
            Int(pattern) == digit && return pattern
        end
    end
    error("digit not found")
end

function decode!(entry::Entry)

    # d_9 => d_4 & 6 segments
    pattern_9 = Display()
    pattern_4 = get_pattern(entry, 4)
    for pattern in patterns(entry)
        if !decoded(pattern) && length(signal(pattern)) == 6 && contains(pattern, pattern_4)
            pattern.digit = 9
            pattern_9 = pattern
            break
        end
    end

    # e => d_8 - d_9
    segment_e = nothing
    pattern_8 = get_pattern(entry, 8)
    for segment in signal(pattern_8)
        if !(segment in signal(pattern_9))
            segment_e = segment
            break
        end
    end

    # d_2 => e & 5 segments
    pattern_2 = Display()
    for pattern in patterns(entry)
        if !decoded(pattern) && length(signal(pattern)) == 5 && segment_e in signal(pattern)
            pattern.digit = 2
            pattern_2 = pattern
            break
        end
    end

    # c => in d_1 & in d_2
    pattern_1 = get_pattern(entry, 1)
    segment_c = only(signal(pattern_1) ∩ signal(pattern_2))

    # d_0 => c & 6 segments
    for pattern in patterns(entry)
        if !decoded(pattern) && length(signal(pattern)) == 6 && segment_c in signal(pattern)
            pattern.digit = 0
            break
        end
    end

    # d_6 => 6 segments
    pattern_6 = Display()
    for pattern in patterns(entry)
        if !decoded(pattern) && length(signal(pattern)) == 6
            pattern.digit = 6
            pattern_6 = pattern
            break
        end
    end

    # d_5 => in d_6
    for pattern in patterns(entry)
        if !decoded(pattern) && contains(pattern_6, pattern)
            pattern.digit = 5
            break
        end
    end

    # d_3 => only one left
    for pattern in patterns(entry)
        if !decoded(pattern)
            pattern.digit = 3
            break
        end
    end

    # match outputs with patterns
    for display in output(entry)
        for pattern in patterns(entry)
            if display == pattern
                display.digit = pattern.digit
                break
            end
        end
    end

    return Int(entry)

end

########## Script ##########

function main()

    # open file and store content into an input string
    input = open("inputs//day8.txt", "r") do io
        read(io, String)
    end

    # split input into a data vector
    data = String.(split(input, "\r\n"))

    # parse data as a vector of entries
    entries = Entry.(data)

    # part 1
    result = 0
    for entry in entries
        for display in entry.output
            result += decoded(display)
        end
    end
    println("Part 1: ", result)

    # part 2
    outputs_number = decode!.(entries)
    println("Part 2: ", sum(outputs_number))

end

main()


text:
# Functions to calculate basic quantities for stellar spectra.
# For these models this is only the Q_0, <hν>_i, L_i, L_n, β quantities.
#
using Unitful, UnitfulRecipes, UnitfulAstro
using Dierckx

"Interpolate a spectrum. The routine strips units so the user needs to take care with that"
function interpolate_spectrum(λ, flux)
    itp = Spline1D(ustrip(λ), ustrip(flux));
    return itp
end
"""
  res = specQuantities(λ, itp)

Calculate various quantities for an interpolated spectrum.

    λ: Wavelength array, assumed to be in Angstrom
  itp: Interpolated spectrum returned from the interpolate_spectrum function above.

Returns:
  res: Named tuple with the following keys:
     Ln: Luminosity in photons with E<13.6 eV  [erg/s]
     Li: Luminosity in photons with E>13.6 eV  [erg/s]
     Q0: Number of photons with E>13.6 eV      [1/s]
   hnui: Mean energy in photons with E>13.6eV  [eV]
      β: Ln/Li
"""
function specQuantities(λ, itp)

    max_l = maximum(λ)
    min_l = minimum(λ)
    lambda_i = 911.267

    norm, err = quadgk(itp, min_l, max_l)
    Ln, err = quadgk(itp, lambda_i, max_l)
    Li, err = quadgk(itp, min_l, lambda_i)

    β = Ln/Li
    # The units are expected to be erg/s/angstrom for the flux so
    # this gives the units for the luminosities
    Ln = Ln*u"erg/s"
    Li = Li*u"erg/s"

    m1_int(x) = itp(x)/x
    mm1_int(x) = x*itp(x)
    norm, err = quadgk(itp, min_l, lambda_i)
    pre = Unitful.h*Unitful.c
    hnui, err = quadgk(m1_int, min_l, lambda_i)
    hnui = pre*hnui/norm
    Q0, err = quadgk(mm1_int, min_l, lambda_i)
    Q0 = Q0*u"erg/s"/pre

    # Finally adjust units because the integrals are done unitless so
    # we need to multiply my Angstrom (except for Ln/Li which are done above)
    hnui = uconvert(u"eV", hnui*1u"1/angstrom")
    Q0 = upreferred(Q0*1u"angstrom")

    return (Ln=Ln, Li=Li, hnui=hnui, Q0=Q0, β=β)
end

#
# Black-body functions below
#

function Bnu(ν, T)
    z = Unitful.h.*ν./(Unitful.k.*T)
    2 .*Unitful.h.*ν^3/(Unitful.c^2 .*(exp(z).-1))
end

function Blambda(λ, T)
    z = Unitful.h.*Unitful.c ./(Unitful.k.*T*λ)
    2 .*Unitful.h.*Unitful.c^2/(λ^5 .*(exp(z).-1))
end

"""
   PlanckAstro(λ, T; convert_to_L=true)

The Planck function in astronomical units and following planck.pro in IDL

Arguments:
    λ: Wavelength (must have units)
    T: Temperature (must have units)

Keywords:
    convert_to_L: Whether the result should be converted to a luminosity.
                If set to true, the radius-T_eff relationship from
                Sternberg, Hoffmann & Pauldrach (2003) is used

"""
function PlanckAstro(λ, T; convert_to_L=true)
    Bl = Blambda(λ, T)

    # Now convert to erg/s/cm^3/A
    Bl = π*uconvert(u"erg/s/cm^2/angstrom", Bl)

    # And then to luminosities if requested
    if (convert_to_L)
        Bl=uconvert(u"erg/s/angstrom", Bl*4*π*radius_for_T(T)^2)
    end
    return ustrip(Bl)
end

#
# Finally we need a T_eff-radius relation to calculate luminosities.
#
"The radius relation for dwarfs from Sternberg, Hoffmann & Pauldrach (2003) "
radius_for_T(T) = (0.2*ustrip(T)/1e3+2)*u"Rsun"


text:
"""
parse_options(kwargs)

Internal function. Takes the keyword arguments from the main function and parses it into a
    usable Dict object

# Examples
```julia-repl
julia> parse_options(ex::Expr)
Dict{String,Any} with 2 entries:
  "screen_name" => "jack"
  ...
```
"""
function parse_options(kwargs)
    options = Dict{String, Any}()
    for arg in kwargs
        options[string(arg[1])] = arg[2]
    end
    options
end

"""
parse_results(cursorable, newdata::Dict, api_options, data_holder, cur_count)

Internal function. parses Twitter API results by determining the type of data and organizing
    for cursorized processing.

    There are two methods, this one takes a DICT object, indicating a set of user IDs or search results.
    returns cursorable, newdata, api_options, cur_count.

# Examples
```julia-repl
julia> parse_results(cursorable, newdata::Dict, api_options, data_holder, cur_count)
  ...
```
"""
function parse_results(cursorable, newdata::Dict, api_options, data_holder, cur_count)
    # Handle data type to cursor
    if haskey(newdata, "ids")
        newdata["ids"] = vcat(data_holder, newdata["ids"])
        cur_count += length(newdata["ids"])
        cursorable = (newdata["next_cursor"] != 0) & (cur_count < api_options["count"] )
        api_options["cursor"] = newdata["next_cursor"]
    elseif haskey(newdata, "statuses")
        out = [Tweets(x) for x in newdata["statuses"]]
        newdata["statuses"] = vcat(data_holder, out)
        cur_count += length(newdata["statuses"])
        cursorable = cur_count < api_options["count"]
        api_options["max_id"] = minimum(x.id for x in newdata["statuses"])
    else
        cur_count += length(newdata)
        cursorable = cur_count < api_options["count"]
    end
    cursorable, newdata, api_options, cur_count
end

"""
parse_results(cursorable, newdata::Array, api_options, data_holder, cur_count)

Internal function. parses Twitter API results by determining the type of data and organizing appropriately.
    There are two methods, this one takes an ARRAY object, indicating a set of Tweets.
    returns cursorable, newdata, api_options, cur_count

# Examples
```julia-repl
julia> parse_results(cursorable, newdata::Array, api_options, data_holder, cur_count)
  ...
```
"""
function parse_results(cursorable, newdata::Array, api_options, data_holder, cur_count)
    newdata = Tweets[Tweets(x) for x in newdata]
    length(newdata) == 0 && return false, data_holder, api_options, cur_count
    # tree of options for max_id or since id
    cur_count += length(newdata)
    cursorable = cur_count < api_options["count"]
    api_options["max_id"] = minimum([x.id for x in newdata])-1  # get min id
    newdata = vcat(data_holder, newdata)
    cursorable, newdata, api_options, cur_count
end
"""
cursor(cursorable::Bool, newdata::Dict, options::Dict, endp::String, cur_count::Integer)

Internal function method for gathering IDS. Takes a tuple and returns a tuple of equal size, calls the Twitter API
    until the desired count of records is recovered or until the API exhausts its limits.

    Note: when a DICT object is provided as the data, this function assumes you are gathering
    follower or friends ids.

# Examples
```julia-repl
julia>  while cursorable & (length(newdata["ids"]) < min_records)
            cursorable, newdata, options, kwargs, endp = cursor(cursorable, newdata, options, kwargs, endp)
        end
```
"""
################# cursor when new data is a Dict object - like followers or friends IDS
function cursor(cursorable::Bool, newdata::Dict, options::Dict, endp::String, cur_count::Integer)
    cursorable == false && return cursorable, newdata, options, endp, cur_count
    data_holder = haskey(newdata, "ids") ? newdata["ids"] : haskey(newdata, "statuses") ? newdata["statuses"] : [] # save existing ids
    api_options = copy(options) # the get_oauth overwrites options, so store the correct data here
    cur_alloc = reconnect("$endp") # start reconnect loop
    remaining_calls = cur_alloc["remaining"]
    @debug "$remaining_calls calls left on this endpoint."
    r = get_oauth("https://api.twitter.com/1.1/$endp", options)
    if r.status == 200
        newdata = JSON.parse(String(r.body))
        cursorable, newdata, api_options, cur_count = parse_results(cursorable, newdata, api_options, data_holder, cur_count)
        cursorable, newdata, api_options, endp, cur_count
    else
        error("Twitter API returned $(r.status) status")
    end
end

"""
cursor(cursorable::Bool, newdata::Dict, options::Dict, endp::String, cur_count::Integer)

Internal function for gathering . Takes a tuple and returns a tuple of equal size, calls the Twitter API
    until the desired count of records is recovered or until the API exhausts its limits.

    Note: when an ARRAY object is provided as the data, this function assumes you are gathering
    a tweet timeline.

    Note: when  retrieving tweets, the API always starts with the most recent. Therefore,
    if you want a chunk of older tweets, you must specify both since_id, and max_id when cursoring.

# Examples
```julia-repl
julia>  while cursorable & (length(newdata["ids"]) < count)
            cursorable, newdata, options, endp = cursor(cursorable, newdata, options, endp)
        end
```
"""
function cursor(cursorable::Bool, newdata::Array, options::Dict, endp::String, cur_count::Integer)
    cursorable == false && return cursorable, newdata, options, endp, cur_count
    data_holder = copy(newdata) # save existing ids
    api_options = copy(options) # the get_oauth overwrites options, so store the correct data here
    cur_alloc = reconnect("$endp") # start reconnect loop
    remaining_calls = cur_alloc["remaining"]
    @debug "$remaining_calls calls left on this endpoint."
    r = get_oauth("https://api.twitter.com/1.1/$endp", options)
    if r.status == 200
        # parse and put into proper type form
        newdata = JSON.parse(String(r.body))
        cursorable, newdata, api_options, cur_count = parse_results(cursorable, newdata, api_options, data_holder, cur_count)
        cursorable, newdata, api_options, endp, cur_count
    else
        error("Twitter API returned $(r.status) status")
    end
end

########### EXPORTED FUNCTIONS.......

"""
get_followers_ids(; kwargs...)
Get a Dict object of follower ids from a particular Twitter user. This function will call the API as
many times as allowed or until the desired `max_records` is reached, whichever comes first.
# Examples
```julia-repl
julia> get_followers_ids(screen_name = "jack", count = 10_000)
Dict{String,Any} with 6 entries:
  "previous_cursor_str" => "0"
  ...
```
"""
function get_followers_ids(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "followers/ids.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1 # default to one record
    end

    cur_count = 0
    cursorable = true
    newdata = Dict{String,Any}()
    newdata["ids"] = [] #Array{String,1}[]

    while cursorable & (length(newdata["ids"]) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

"""
get_friends_ids(; kwargs...)
Get a Dict object of follower ids from a particular Twitter user. This function will call the API as
until the desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_friends_ids(screen_name = "barackobama", count = 1000)
```
"""
function get_friends_ids(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "friends/ids.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1 # default to one record
    end

    cur_count = 0
    cursorable = true
    newdata = Dict{String,Any}()
    newdata["ids"] = [] #Array{String,1}[]

    while cursorable & (length(newdata["ids"]) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

########################## OTHER TYPE:

"""
get_mentions_timeline(; kwargs...)
Get an array object of mentions for a particular Twitter user. This function will call the API until the
desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_mentions_timeline(screen_name = "twitter", count = 1000)
```
"""
function get_mentions_timeline(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "statuses/mentions_timeline.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1
    end
    cur_count = 0
    # make the first call to the API
    cursorable = true
    newdata = Tweets[]

    while cursorable & (length(newdata) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

"""
get_user_timeline(; kwargs...)
Get an array object of timeline tweets from a particular Twitter user. This function will call the API until the
desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_user_timeline(screen_name = "twitter", count = 1000)
```
"""
function get_user_timeline(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "statuses/user_timeline.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1
    end
    cur_count = 0
    # make the first call to the API
    cursorable = true
    newdata = Tweets[]

    while cursorable & (length(newdata) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

"""
get_home_timeline(; kwargs...)
Get an array object of timeline tweets from the owning user. This function will call the API until the
desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_home_timeline(count = 1000)
```
"""
function get_home_timeline(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "statuses/home_timeline.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1
    end
    cur_count = 0
    # make the first call to the API
    cursorable = true
    newdata = Tweets[]

    while cursorable & (length(newdata) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

"""
get_retweets_of_me(; kwargs...)
Get an array object of retweets from the owning user. This function will call the API until the
desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_retweets_of_me(count = 1000)
```
"""
function get_retweets_of_me(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "statuses/retweets_of_me.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1
    end
    cur_count = 0
    # make the first call to the API
    cursorable = true
    newdata = Tweets[]

    while cursorable & (length(newdata) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end

"""
get_retweets_of_me(; kwargs...)
Get an array object of retweets from the owning user. This function will call the API until the
desired `count` is reached or the API runs out, whichever comes first.
# Examples
```julia-repl
julia> get_retweets_of_me(count = 1000)
```
"""
function get_search_tweets(; kwargs...)
    # Could be doing some pre-allocation here to optimize performance,
    # but since this is an API function that only deals with 25K records at most...
    endp = "search/tweets.json"
    options = parse_options(kwargs)

    if "count" ∈ keys(options)
        count = options["count"]
    else
        options["count"] = 1
        count = 1
    end
    cur_count = 0
    # make the first call to the API
    cursorable = true
    newdata = Tweets[]

    while cursorable & (length(newdata) < count)
        cursorable, newdata, options, endp, cur_count = cursor(cursorable, newdata, options, endp, cur_count)
    end
    newdata
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  S'(t) = -b * S(t) * In(t) / N(t),
  E'(t) = b * S(t) * In(t) / N(t) - nu * E(t),
  In'(t) = nu * E(t) - a * In(t),
  N'(t) = 0,
  Cu'(t) = nu * E(t),
  y1(t) = Cu(t),
  y2(t) = N(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
### File for testing that the equilibrate function finds the proper steady states and stabilities.

#Tests on the standard networks.
for rn in reaction_networks_standard                                                #For every reaction network in the set.
    p_vals = map(amp->amp*rand(length(rn.params)),[1.,10.,100.])                    #Generates three random sets of parameter values.
    for p in p_vals                                                                 #For each parameter value set.
        fps = steady_states(rn,p,stability_type=true)                               #Finds the fixed points..
        stab = map(fp->stability(fp[1],rn,p), fps)                                  #Calculates stability using the stability calculation function
        @test [entry[2] for entry in fps]==stab                                     #Tests that the stability values corresponds.
        for fp in fps                                                               #For each fixed point found.
            fp[2][1] ? (tend = 10) : (tend = -10)                                   #Checks if fixed point is unstable, in that case we will run solver in negative time.
            end_point = solve(ODEProblem(rn,fp[1],(0.,tend),p),Rosenbrock23())[end] #Starts a simulation in the fixed point found.
            @test maximum(end_point-fp[1])<0.0001                                   #Checks that we are still at the same point.
        end
    end
end

#Tests on the networks containing hill functions.
for rn in reaction_networks_hill
    p_vals = map(amp->rand(1:amp,length(rn.params)),[3,5,10])
    for p in p_vals
        fps = steady_states(rn,p,stability_type=true)
        stab = map(fp->stability(fp[1],rn,p), fps)
        @test [entry[2] for entry in fps]==stab
        for fp in fps
            fp[2][1] ? (tend = 10) : (tend = -10)
            end_point = solve(ODEProblem(rn,fp[1],(0.,tend),p),Rosenbrock23())[end]
            @test maximum(end_point-fp[1])<0.0001
        end
    end
end

#Tests on the networks requiring fixed concentrations.
for i = 1:length(reaction_networks_fixed_conc)
    p_vals = map(amp->amp*rand(length(reaction_networks_fixed_conc[i].params)),[1.,10.,100.])
    for p in p_vals
        fps = steady_states(reaction_networks_fixed_conc[i],p,stability_type=true,fixed_conc[i]...)
        stab = map(fp->stability(fp[1],reaction_networks_fixed_conc[i],p), fps)
        @test [entry[2] for entry in fps]==stab
        for fp in fps
            fp[2][1] ? (tend = 10) : (tend = -10)
            end_point = solve(ODEProblem(reaction_networks_fixed_conc[i],fp[1],(0.,tend),p),Rosenbrock23())[end]
            @test maximum(end_point-fp[1])<0.0001
        end
    end
end

max_stars_repo_name: TorkelE/BRNEquilibrate

max_issues_repo_name: TorkelE/BRNEquilibrate

max_forks_repo_name: TorkelE/BRNEquilibrate


text:
using Gurobi, MathProgBase

env = Gurobi.Env()

# Test supplying environment manually

solver = GurobiSolver(env)
# Check that the env for each model is the same
m1 = MathProgBase.LinearQuadraticModel(solver)
@assert m1.inner.env === env
m2 = MathProgBase.LinearQuadraticModel(solver)
@assert m2.inner.env === env
# Check that finalizer doesn't touch env when manually provided
finalize(m1.inner)
@assert Gurobi.is_valid(env)

# Test creating environment automatically

solver = GurobiSolver()
# Check that the env for each model is different
m3 = MathProgBase.LinearQuadraticModel(solver)
@assert m3.inner.env !== env
m4 = MathProgBase.LinearQuadraticModel(solver)
@assert m4.inner.env !== env
@assert m3.inner.env !== m4.inner.env
# Check that env is finalized with model when not supplied manually
finalize(m3.inner)
@assert !Gurobi.is_valid(m3.inner.env)


text:
# ==============================================================================
# Parseur lisant une instance de probleme de partitionnement d'ensembles bi-objectif (2-SPA)

function loadInstance2SPA(fname::String)

    f = open("../SPA/instances/"*"bio"*fname)
    nbctr, nbvar = parse.(Int, split(readline(f))) # nombre de contraintes , nombre de variables
    A = zeros(Int, nbctr, nbvar)                   # matrice des contraintes
    c1 = zeros(Int, nbvar)                         # vecteur des couts
    c2 = zeros(Int, nbvar)                         # deuxième vecteur des couts
    nb = zeros(Int, nbvar)
    for i in 1:nbvar
        flag = 1
        for valeur in split(readline(f))
            if flag == 1
                c1[i] = parse(Int, valeur)
                flag +=1
            elseif flag == 2
                c2[i] = parse(Int, valeur)
                flag +=1
            elseif flag == 3
                nb[i] = parse(Int, valeur)
                flag +=1
            else
                j = parse(Int, valeur)
                A[j,i] = 1
            end
        end
    end
    close(f)
    return c1, c2, A
end

# ==============================================================================
# Parseur lisant un fichier de Y_N (dans dossier SPA/Y) pour le 2-SPA traite

function loadNDPoints2SPA(fname::String)
    fname = "../SPA/Y/Y_N_"*fname
    f = open(fname)
    temps = parse.(Float64, readline(f))
    len = parse.(Int, readline(f))
    xN = Vector{Int64}(undef, len)
    yN = Vector{Int64}(undef, len)
    for i in 1:len
        x, y =  parse.(Float64, split(readline(f)))
        xN[i],yN[i] = round(Int,x),round(Int,y)
    end
    close(f)
    return xN, yN
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
include("src/toolbox.jl")
include("src/augmentation.jl")
include("src/solver.jl")
include("src/cvxreg.jl")
include("src/testdata.jl")
using DataFrames, JLD,CSV,JSON
using DataStructures

n = 10000;
d = 4;
logrho = -3;
rho = 10.0^logrho;
func_type = "A";
random_state = 22;
ntrain = n;
ntest = 10000;
r = (d <= 3) ? 0.8 : 0.9;
SNR = 3;

LBFGS = false;
rule = 4;
LS = true;
WS = false;
inexact=false;

xtrain, ytrain, xtest, ytest, xtest_bd, ytest_bd, xtest_int, ytest_int = generate_training_test_samples(func_type,ntrain,ntest,d,r,SNR,1)

maxiter=1000
maxtime = 1800
verbose = 1
block = 0
Xmat = xtrain;
Y = ytrain;
Xflat = mat_to_flat(Xmat,n,d);
violTOL = 1.0e-4;
innerTOL = 1.0e-6;
outerTOL = 1.0e-4;
dropzero = true

params1 = AlgorithmParameters(inexact=true, violTOL=violTOL, outerTOL=outerTOL)
params2 = AlgorithmParameters(inexact=inexact, violTOL=1.0e-8, outerTOL=outerTOL)
aug1 = ActiveSetAugmentation(rule,n,block,params1.violTOL)
aug2 = ActiveSetAugmentation(rule,n,block,params2.violTOL)
runtime = @elapsed phi,xi,lamb, W, Wlen, all_objs,all_obj_times, stages,Wlen_list,searching_times,pow_times,f_evals_list, L_list, optimization_starting_steps =
two_stage_active_set_profiling(Xmat,Y,rho;
    verbose = verbose,random_state=random_state, maxiter=maxiter, maxtime=maxtime, block=block, dropzero=dropzero,  params_list = [params1, params2], augs = [aug1, aug2]);

if block == 0
    J,II,_= findnz(W)
else
    II,J,_= findnz(W)
end

pinfeas, max_viol = compute_infeasibility(Xflat,n,d,phi,xi)
gap,phi_f,xi_f = get_duality_gap(Xflat,Y,rho,II,J, Wlen, lamb.nzval,phi,xi, n, d,return_feasible_soln=true)

lamb_sparse = dropzeros(lamb)
if block ==0
    J_sparse,I_sparse,_ = findnz(lamb_sparse)
else
    I_sparse,J_sparse,_ = findnz(lamb_sparse)
end
Wlen_sparse = length(lamb_sparse.nzval)
W_sparse = sparse(J_sparse,I_sparse,fill(true,Wlen_sparse))

ypred_is = cvx_predict(xtrain, Xflat,phi_f,xi_f, n, d);
ypred = cvx_predict(xtest, Xflat,phi_f,xi_f, n, d);
ypred_bd = cvx_predict(xtest_bd, Xflat,phi_f,xi_f, n, d);
ypred_int = cvx_predict(xtest_int, Xflat,phi_f,xi_f, n, d);

rmse_is = rmse(ypred_is, ytrain)
rsq_is = rsquared(ypred_is, ytrain)
rmse_oos = rmse(ypred, ytest)
rsq_oos = rsquared(ypred, ytest)
rmse_bd = rmse(ypred_bd, ytest_bd)
rsq_bd = rsquared(ypred_bd, ytest_bd)
rmse_int = rmse(ypred_int, ytest_int)
rsq_int = rsquared(ypred_int, ytest_int)


text:
using SegyIO
using DelimitedFiles

##
# choose the vertical resampling interval
yres = 4.0;
max_depth = 500.;

##
infile = raw"C:\Users\u77932\Documents\LEB\data\AEM\Frome\galei\lines\Frome_LEB_galei.dat";
outdir = raw"C:\Users\u77932\Documents\LEB\data\AEM\Frome\galei\segy";

inversion_output = readdlm(infile, Float64);

# define column numbers
easting_col = 7;
northing_col = 8;
elevation_col = 9;
conductivity_cols = [24,53] #ranger for 30 layer model;
thickness_cols = [54,83];
#layer_top_elevation_cols = [14 43]
line_col = 5;

lines = round.(Int, inversion_output[:,line_col]);

##

function get_line_indices(line)
    inds = findall(x -> x in line, lines)
    return inds
end

function line2segy(line_number)
    # find our actual line indices
    line_inds = get_line_indices(line_number)

    # our first step is to define our datum as 20 m above the highest point in survey

    elevations = inversion_output[line_inds,elevation_col]
    ymax = ceil(maximum(elevations)/10.) *10.

    # create an elevation grid

    ymin = floor(minimum(elevations)/10.)*10. - max_depth

    grid_elevations = collect(range(ymax, ymin, length = Int(ceil((ymax - ymin)/yres))))

    # get conductivity
    sigma = transpose(inversion_output[line_inds,conductivity_cols[1]:conductivity_cols[2]])
    thickness = transpose(inversion_output[line_inds,thickness_cols[1]:thickness_cols[2]])
    #layer_top_elevation = inversion_output[line_inds,layer_top_elevation_cols[1]:layer_top_elevation_cols[2]]

    ## get coordinates

    easting = inversion_output[line_inds,easting_col]
    northing = inversion_output[line_inds,northing_col]
    ## create our array of interpolated AEM conductivity, negative indicates free air
    interpolated = -1*ones(Float32, length(grid_elevations), length(line_inds))

    # iterate
    for i in 1:length(line_inds)

        layer_top_elev = elevations[i] * ones(size(thickness)[1])
        layer_top_elev[2:end] = layer_top_elev[2:end] - cumsum(thickness[1:end-1,i])
        #layer_top_elev = layer_top_elevation[i,:]
        sigma_trace = sigma[:,i]
        # iterate through our layers
        for j in 1:(length(layer_top_elev) - 1)
            interpolated[(grid_elevations .< layer_top_elev[j]) .& (grid_elevations .>= layer_top_elev[j + 1]),i] .= sigma_trace[j]
        end
    end

    ## try writing these out in depth space

    block = SeisBlock(interpolated)

    set_header!(block, "SourceX", round.(Int, easting))
    set_header!(block, "SourceY", round.(Int, northing))
    set_header!(block, "GroupX", round.(Int, easting))
    set_header!(block, "GroupY", round.(Int, northing))
    set_header!(block, "ElevationScalar", round.(Int, elevations))
    set_header!(block, "TraceNumWithinLine", Array(1:length(line_inds)))
    set_header!(block, "TraceNumWithinFile", Array(1:length(line_inds)))
    #set_header!(block, "TraceIDCode", round.(Int, fiducial * 10))
    set_header!(block, "RecGroupElevation", round.(Int, elevations))
    set_header!(block, "DelayRecordingTime", Int(ymax) * -1)
    set_header!(block, "dt", Int(yres * 1e3))

    ##
    segy_write(joinpath(outdir,string(line_number) * ".segy"), block)
end
##
for l in unique(lines)

    line2segy(l)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""OwnerReference contains enough information to let you identify an owning object. An owning object must be in the same namespace as the dependent, or be cluster-scoped, so there is no namespace field.

    IoK8sApimachineryPkgApisMetaV1OwnerReference(;
        apiVersion=nothing,
        blockOwnerDeletion=nothing,
        controller=nothing,
        kind=nothing,
        name=nothing,
        uid=nothing,
    )

    - apiVersion::String : API version of the referent.
    - blockOwnerDeletion::Bool : If true, AND if the owner has the \"foregroundDeletion\" finalizer, then the owner cannot be deleted from the key-value store until this reference is removed. Defaults to false. To set this field, a user needs \"delete\" permission of the owner, otherwise 422 (Unprocessable Entity) will be returned.
    - controller::Bool : If true, this reference points to the managing controller.
    - kind::String : Kind of the referent. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - name::String : Name of the referent. More info: http://kubernetes.io/docs/user-guide/identifiers#names
    - uid::String : UID of the referent. More info: http://kubernetes.io/docs/user-guide/identifiers#uids
"""
mutable struct IoK8sApimachineryPkgApisMetaV1OwnerReference <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    blockOwnerDeletion::Any # spec type: Union{ Nothing, Bool } # spec name: blockOwnerDeletion
    controller::Any # spec type: Union{ Nothing, Bool } # spec name: controller
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid

    function IoK8sApimachineryPkgApisMetaV1OwnerReference(;apiVersion=nothing, blockOwnerDeletion=nothing, controller=nothing, kind=nothing, name=nothing, uid=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("blockOwnerDeletion"), blockOwnerDeletion)
        setfield!(o, Symbol("blockOwnerDeletion"), blockOwnerDeletion)
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("controller"), controller)
        setfield!(o, Symbol("controller"), controller)
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApimachineryPkgApisMetaV1OwnerReference, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1OwnerReference

const _property_map_IoK8sApimachineryPkgApisMetaV1OwnerReference = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("blockOwnerDeletion")=>Symbol("blockOwnerDeletion"), Symbol("controller")=>Symbol("controller"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("uid")=>Symbol("uid"))
const _property_types_IoK8sApimachineryPkgApisMetaV1OwnerReference = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("blockOwnerDeletion")=>"Bool", Symbol("controller")=>"Bool", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("uid")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1OwnerReference }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1OwnerReference))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1OwnerReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1OwnerReference[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1OwnerReference }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1OwnerReference[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1OwnerReference)
    (getproperty(o, Symbol("apiVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("uid")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1OwnerReference }, name::Symbol, val)
end


text:
"""
    objective_function_template(nobs::Function,
                                obj_contribution::Function)

A constructor of objects of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) for defining an `objective_function_template`.

Arguments
===
+ `nobs`: a function of `data` that computes the number of observations of the particular data type,
+ `observation`: a function of `data` and an indexing value `i` which extracts the `i`th observation from the data
+ `obj_contribution`: a function of the parameters `theta` observation (calculated from `data`, and index `i` using the `observation` above) that returns a `Float64`.

Result
===
An `objective_function_template` object with fields `nobs` and `obj_contributions`.

"""
struct objective_function_template
    nobs::Function
    observation::Function
    obj_contribution::Function
end

"""
    objective_function(theta::Vector,
                       data::Any,
                       template::objective_function_template,
                       br::Bool = false)

Evaluates the objective function at `theta` by adding up all contributions in
`data`, according to [`objective_function_template`](@ref).

Arguments
===
+ `theta`: a `Vector` of parameter values at which to evaluate the objective function
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `objective_function`.
+ `template`: an [`objective_function_template`](@ref) object
+ `br`: a `Bool`. If `false` (default), the objective function is constructed by adding up all contributions in
`data`, according to [`objective_function_template`](@ref), before it is evaluated at `theta`. If `true` then the bias-reducing penalty in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) is computed and added to the objective function.

Result
===
A `Float64`.

Details
===
`data` can be used to pass additional constants other than the actual data to the objective.
"""
function objective_function(theta::Vector,
                            data::Any,
                            template::objective_function_template,
                            br::Bool=false)
    p = length(theta)
    n_obs = template.nobs(data)
    objective = 0
    for i in 1:n_obs
        obs_i = template.observation(data, i)
        objective += template.obj_contribution(theta, obs_i)
    end
    if (br)
        objective + obj_quantities(theta, data, template, true)[1]
    else
        objective
    end
end

function obj_quantities(theta::Vector,
                        data::Any,
                        template::objective_function_template,
                        penalty::Bool=false)

    objective_i(pars, obs) = template.obj_contribution(pars, obs)

    function objective_gradient(theta::Vector{T}, obs) where {T <: Real}
        x = (theta, obs)
        f = (x0, x1) -> objective_i(x0, x1)
        if !(haskey(GCACHE, T))
            tape = ReverseDiff.compile(ReverseDiff.GradientTape(f, x))
            GCACHE[T] = (tape, (zeros(T, length(theta)), zeros(T, length(obs))))
        end
        tape, y = GCACHE[T]
        return ReverseDiff.gradient!(y, tape, x)[1]
    end

    function objective_hessian(pars, obs)
        ForwardDiff.jacobian(eta -> objective_gradient(eta, obs), pars)
    end
    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)
    for i in 1:n_obs
        obs_i = template.observation(data, i)
        cpsi = objective_gradient(theta, obs_i)
        psi += cpsi
        emat += cpsi * cpsi'
        jmat += objective_hessian(theta, obs_i)
    end
    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv)
    if (penalty)
        br_penalty = - tr(jmat_inv * emat) / 2
        [br_penalty, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end

function obj_quantities_old(theta::Vector,
                            data::Any,
                            template::objective_function_template,
                            penalty::Bool=false)
    function gr_i(eta::Vector, i::Int)
        out = similar(eta)
        ForwardDiff.gradient!(out, beta -> template.obj_contribution(beta, data, i), eta)
    end
    function he_i(eta::Vector, i::Int)
        out = similar(eta, p, p)
        ForwardDiff.hessian!(out, beta -> template.obj_contribution(beta, data, i), eta)
    end
    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)
    for i in 1:n_obs
        cpsi = gr_i(theta, i)
        psi += cpsi
        emat += cpsi * cpsi'
        jmat += - he_i(theta, i)
    end
    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv)
    if (penalty)
        br_penalty = - tr(jmat_inv * emat) / 2
        # br_penalty = n_obs * log(det(Matrix{Float64}(I * n_obs, p, p) -
        #                              jmat_inv * emat)) / 2
        # br_penalty = + log(det(sum(njmats))) / 2 - log(det(emat)) / 2
        [br_penalty, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end

"""
    estimating_function_template(object::objective_function_template)

Construct an [`estimating_function_template`](@ref) from an [`objective_function_template`](@ref).

Arguments
===
+ `object`: an [`objective_function_template`](@ref)

Details
===
The field `ef_contribution` of the result is computed by differentiating (using the [**ForwardDiff**](https://github.com/JuliaDiff/ForwardDiff.jl) package) `object.obj_contribution` with respect to `theta`.

Result
===
A [`estimating_function_template`](@ref) object.

"""
function estimating_function_template(obj_template::objective_function_template)
    function ef_contribution(theta::Vector{T}, obs) where {T <: Real}
        x = (theta, obs)
        f = obj_template.obj_contribution
        if !(haskey(EOCACHE, T))
            tape = ReverseDiff.compile(ReverseDiff.GradientTape(f, x))
            EOCACHE[T] = (tape, (zeros(T, length(theta)), zeros(T, length(obs))))
        end
        tape, y = EOCACHE[T]
        return ReverseDiff.gradient!(y, tape, x)[1]
    end
    return estimating_function_template(obj_template.nobs, obj_template.observation, ef_contribution)
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/CH4forcing.jl")

    add_comp!(m, ch4forcing, :ch4forcing)

    set_param!(m, :ch4forcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/$valdir/c_n2oconcentration.csv"))
    set_param!(m, :ch4forcing, :c_CH4concentration, readpagedata(m,"test/validationdata/$valdir/c_ch4concentration.csv"))

    ##running Model
    run(m)

    @test !isnan(m[:ch4forcing, :f_CH4forcing][10])

    forcing=m[:ch4forcing,:f_CH4forcing]
    forcing_compare=readpagedata(m,"test/validationdata/$valdir/f_ch4forcing.csv")

    @test forcing ≈ forcing_compare rtol=1e-3
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using Revise
using JDF

a=gf()
b = Array(a.ORIGIN_STATE_ABR)
io = iow()
@which compress_then_write(eltype(b), b, io)
@time metadata = compress_then_write(eltype(b), b, io)
close(io)

io = ior()
@which column_loader(eltype(b), b, io)

@time okhl = column_loader(eltype(b), io, metadata)
close(io)

all(skipmissing(b) .== skipmissing(okhl))

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
export GramMatrix, SparseGramMatrix

import MultivariateMoments: trimat, SymMatrix, getmat
export gram_operate, getmat

abstract type AbstractGramMatrix{T, B, U} <: MP.AbstractPolynomialLike{U} end

function MP.monomialtype(::Union{AbstractGramMatrix{T, B},
                                 Type{<:AbstractGramMatrix{T, B}}}) where {T, B}
    return MP.monomialtype(B)
end
MP.polynomialtype(::Union{AbstractGramMatrix{T, B, U}, Type{<:AbstractGramMatrix{T, B, U}}}) where {T, B, U} = MP.polynomialtype(B, U)

Base.:(+)(x::MP.APL, y::AbstractGramMatrix) = x + MP.polynomial(y)
Base.:(+)(x::AbstractGramMatrix, y::MP.APL) = MP.polynomial(x) + y
Base.:(+)(x::AbstractGramMatrix, y::AbstractGramMatrix) = MP.polynomial(x) + MP.polynomial(y)
Base.:(-)(x::MP.APL, y::AbstractGramMatrix) = x - MP.polynomial(y)
Base.:(-)(x::AbstractGramMatrix, y::MP.APL) = MP.polynomial(x) - y
Base.:(-)(x::AbstractGramMatrix, y::AbstractGramMatrix) = MP.polynomial(x) - MP.polynomial(y)
Base.:(==)(p::MP.APL, q::AbstractGramMatrix) = p == MP.polynomial(q)
Base.:(==)(p::AbstractGramMatrix, q::AbstractGramMatrix) = iszero(p - q)

"""
    struct GramMatrix{T, B, U, MT <: AbstractMatrix{T}} <: AbstractGramMatrix{T, B, U}
        Q::MT
        basis::B
    end

Gram matrix ``x^\\top Q x`` where `Q` is a symmetric matrix indexed by the
vector of polynomials of the basis `basis`.
"""
struct GramMatrix{T, B, U, MT <: AbstractMatrix{T}} <: AbstractGramMatrix{T, B, U}
    Q::MT
    basis::B
end
GramMatrix{T, B, U}(Q::AbstractMatrix{T}, basis::B) where {T, B<:AbstractPolynomialBasis, U} = GramMatrix{T, B, U, typeof(Q)}(Q, basis)
GramMatrix{T, B}(Q::AbstractMatrix{T}, basis::B) where {T, B<:AbstractPolynomialBasis} = GramMatrix{T, B, _promote_sum(eltype(Q))}(Q, basis)
function GramMatrix(
    # We don't use `AbstractMatrix` to avoid clash with method below with `σ`.
    Q::Union{MultivariateMoments.SymMatrix{T}, MultivariateMoments.VectorizedHermitianMatrix{T}},
    basis::AbstractPolynomialBasis) where T
    return GramMatrix{T, typeof(basis)}(Q, basis)
end
# When taking the promotion of a GramMatrix of JuMP.Variable with a Polynomial JuMP.Variable, it should be a Polynomial of AffExpr
MP.constantmonomial(p::GramMatrix) = MP.constantmonomial(MP.monomialtype(p))
MP.variables(p::GramMatrix) = MP.variables(p.basis)
MP.nvariables(p::GramMatrix) = MP.nvariables(p.basis)

Base.zero(::Type{GramMatrix{T, B, U, MT}}) where {T, B, U, MT} = GramMatrix{T, B, U, MT}(similar(MT, (0, 0)), MultivariateBases.empty_basis(B))
Base.iszero(p::GramMatrix) = iszero(MP.polynomial(p))

Base.getindex(p::GramMatrix, I...) = getindex(p.Q, I...)
Base.copy(p::GramMatrix) = GramMatrix(copy(p.Q), copy(p.basis))

MultivariateMoments.getmat(p::GramMatrix{T}) where {T} = p.Q

function GramMatrix{T}(f::Function, basis::AbstractPolynomialBasis, σ=1:length(basis)) where T
    GramMatrix{T, typeof(basis)}(trimat(T, f, length(basis), σ), basis)
end
function GramMatrix{T}(f::Function, monos::AbstractVector) where T
    σ, sorted_monos = MP.sortmonovec(monos)
    return GramMatrix{T}(f, MonomialBasis(sorted_monos), σ)
end

function GramMatrix(Q::AbstractMatrix{T}, basis::AbstractPolynomialBasis, σ=1:length(basis)) where T
    return GramMatrix{T}((i, j) -> Q[σ[i], σ[j]], basis)
end
function GramMatrix(Q::AbstractMatrix, monos::AbstractVector)
    σ, sorted_monos = MP.sortmonovec(monos)
    return GramMatrix(Q, MonomialBasis(sorted_monos), σ)
end

#function Base.convert{T, PT <: AbstractPolynomial{T}}(::Type{PT}, p::GramMatrix)
#    # coefficienttype(p) may be different than T and MP.polynomial(p) may be different than PT (different module)
#    convert(PT, MP.polynomial(p))
#end
function MP.polynomial(p::GramMatrix{T, B, U}) where {T, B, U}
    return MP.polynomial(p, U)
end
function MP.polynomial(p::GramMatrix, ::Type{S}) where S
    return MP.polynomial(getmat(p), p.basis, S)
end

change_basis(p::GramMatrix{T, B}, ::Type{B}) where {T, B<:AbstractPolynomialBasis} = p
function change_basis(p::GramMatrix, B::Type{<:AbstractPolynomialBasis})
    return GramMatrix(MultivariateBases.change_basis(p.Q, p.basis, B)...)
end

"""
    gram_operate(::typeof(+), p::GramMatrix, q::GramMatrix)

Computes the Gram matrix equal to the sum between `p` and `q`. On the opposite,
`p + q` gives a polynomial equal to `p + q`. The polynomial `p + q` can also be
obtained by `polynomial(gram_operate(+, p, q))`.
"""
function gram_operate(::typeof(+), p::GramMatrix{S, B, US, SymMatrix{S}}, q::GramMatrix{T, B, UT, SymMatrix{T}}) where {S, US, T, UT, B}
    basis, Ip, Iq = MultivariateBases.merge_bases(p.basis, q.basis)
    U = MA.promote_operation(+, S, T)
    n = length(basis)
    Qvec = Vector{U}(undef, div(n * (n + 1), 2))
    MA.mutable_operate!(zero, Qvec)
    Q = SymMatrix(Qvec, n)
    for j in 1:n
        for i in 1:j
            if !iszero(Ip[j]) && !iszero(Ip[i])
                MultivariateMoments.symmetric_setindex!(
                    Q, MA.add!(Q[i, j], p.Q[Ip[i], Ip[j]]), i, j)
            end
            if !iszero(Iq[j]) && !iszero(Iq[i])
                MultivariateMoments.symmetric_setindex!(
                    Q, MA.add!(Q[i, j], q.Q[Iq[i], Iq[j]]), i, j)
            end
        end
    end
    return GramMatrix(Q, basis)
end
function gram_operate(::typeof(+), p::GramMatrix{S, Bp}, q::GramMatrix{T, Bq}) where {S, T, Bp, Bq}
    B = promote_type(Bp, Bq)
    return gram_operate(+, change_basis(p, B), change_basis(q, B))
end

"""
    gram_operate(/, p::GramMatrix, α)

Computes the Gram matrix equal to `p / α`. On the opposite,
`p / α` gives a polynomial equal to `p / α`. The polynomial `p / α` can also be
obtained by `polynomial(gram_operate(/, p, α))`.
"""
function gram_operate(::typeof(/), q::GramMatrix, α)
    return GramMatrix(map(x -> x / α, q.Q), q.basis)
end

function Base.isapprox(p::GramMatrix, q::GramMatrix; kwargs...)
    p.basis == q.basis && isapprox(p.Q, q.Q; kwargs...)
end

struct SparseGramMatrix{T, B, U, MT} <: AbstractGramMatrix{T, B, U}
    sub_gram_matrices::Vector{GramMatrix{T, B, U, MT}}
end

Base.zero(::Type{SparseGramMatrix{T, B, U, MT}}) where {T, B, U, MT} = SparseGramMatrix(GramMatrix{T, B, U, MT}[])
function MP.polynomial(p::SparseGramMatrix)
    return mapreduce(identity, MA.add!, p.sub_gram_matrices, init = zero(MP.polynomialtype(p)))
end


text:
__precompile__()

module Configs
    include("utils.jl")
    using JSON, Memoize
    export  getconfig,
            setconfig!,
            hasconfig,
            initconfig

    configs = nothing
    const configs_defaultorder = [
        "default.json",
        "custom-environment-variables.json"
    ]

    function initconfig(;deployment_key = "DEPLOYMENT", configs_directory = "configs")::NamedTuple
        global configs = Dict{String, Any}()
        configs_order = copy(configs_defaultorder)
        configs_directory = parseenvkey("CONFIGS_DIRECTORY", configs_directory)
        deployment_key = parseenvkey("DEPLOYMENT_KEY", deployment_key)
        configs_files = getfiles(configs_directory)
        deployment = parseenvkey(deployment_key, false)
        deployment != false && insert!(configs_order, 2, lowercase(deployment) * ".json")

        filter!((file)-> file in configs_files, configs_order)
        for file in configs_order
            open(joinpath(configs_directory, file), "r") do filepath
                file_content = String(read(filepath))
                newtree = JSON.parse(file_content)
                if file === "custom-environment-variables.json"
                    parsecustomenv!(newtree)
                end
                override!(configs, newtree)
            end
        end
        (; configs_directory = configs_directory, deployment_key = deployment_key, configs_order = configs_order)
    end

    @memoize Dict{Tuple{String}, Any} function getconfig(path::String = "")
        global configs
        configs === nothing && initconfig()
        configs isa Dict && (configs = makeimmutable(configs))
        path === "" && return configs
        epath = Meta.parse("configs.$path")
        try
            eval(epath)
        catch err
            throw(Configserror("no such config: " * path))
        end
    end

    function setconfig!(path::String, value::String)
        try
            value = JSON.parse(value)
        catch err
        end
        value = pathtodict(path, value)
        setconfig!(value)
    end
    function setconfig!(path::String, value)
        value = pathtodict(path, value)
        setconfig!(value)
    end
    function setconfig!(path::String, value::Union{Tuple, AbstractArray, NamedTuple, Dict})
        value = json(value) |> JSON.parse
        value = pathtodict(path, value)
        setconfig!(value)
    end
    function setconfig!(value::Dict)
        configs === nothing && initconfig()
        configs isa NamedTuple && throw(Configserror("""config is immutable. Please set all values before calling "get" """))
        override!(configs, value)
    end

    function hasconfig(path::String)::Bool
        path === "" && throw(Configserror("a path is required to query a config"))
        configs === nothing && initconfig()
        subpaths = split(path, ".")
        ref = configs
        for i in eachindex(subpaths)
            subpath = Symbol(subpaths[i])
            if haskey(ref, subpath)
                length(subpaths) === i && return true
                ref = ref[subpath]
            else
                return false
            end
        end
    end
end

max_stars_repo_name: tisztamo/Configs

max_issues_repo_name: tisztamo/Configs

max_forks_repo_name: tisztamo/Configs