text:
#################### Model Simulation ####################
"""
    gettune(m::Model, block::Integer)
"""
function gettune(m::Model, block::Integer)
    block == 0 && return gettune(m)
    m.samplers[block].tune
end
"""
    gettune(m::Model)

Get block-sampler tuning parameters.

Returns a `Vector{Any}` of all block-specific tuning parameters without `block` input, and turning parameters for the specified `block` otherwise.
"""
function gettune(m::Model)
    Any[gettune(m, i) for i = 1:length(m.samplers)]
end
"""
    settune!(m::Model, tune, block::Integer)
"""
function settune!(m::Model, tune, block::Integer)
    block == 0 && return settune!(m, tune)
    m.samplers[block].tune = tune
end
"""
    settune!(m::Model, tune::Vector{Any})

Set tuning parameters for one or all blocks.

Assigns desired tune values to model.

* `m` : model containing the nodes of interest.

* `tune` : tune values to be assigned to models; if no `block` value is input, `tune` must be a Vector with length equal to `m.samplers`.

* `block` : Integer denoting which block's tune value is to be reassigned.
"""
function settune!(m::Model, tune::Vector{Any})
    nsamplers = length(m.samplers)
    ntune = length(tune)
    nsamplers == ntune || throw(
        DimensionMismatch("tried to assign $ntune tune elements to $nsamplers samplers"),
    )

    for i = 1:nsamplers
        settune!(m, tune[i], i)
    end
end

"""
    gradlogpdf(m::Model, block::Integer=0, transform::Bool=false;
                    dtype::Symbol=:forward)

"""
function gradlogpdf(
    m::Model,
    block::Integer = 0,
    transform::Bool = false;
    dtype::Symbol = :forward,
)
    x0 = unlist(m, block, transform)
    value = gradlogpdf!(m, x0, block, transform, dtype = dtype)
    relist!(m, x0, block, transform)
    value
end

"""
    gradlogpdf!(m::Model, x::AbstractVector{T}, block::Integer=0,
                      transform::Bool=false; dtype::Symbol=:forward) where {T<:Real}

"""
function gradlogpdf!(
    m::Model,
    x::AbstractVector{T},
    block::Integer = 0,
    transform::Bool = false;
    dtype::Symbol = :forward,
) where {T<:Real}
    f = y -> logpdf!(m, y, block, transform)
    if dtype == :Zygote
        mylogpdf!(m, x, block, transform)
        #Zygote.gradient(f, x)
    else
        FiniteDiff.finite_difference_gradient(f, x)
    end

end

"""
    logpdf(m::Model, block::Integer=0, transform::Bool=false)
"""
function logpdf(m::Model, block::Integer = 0, transform::Bool = false)
    params = keys(m, :block, block)
    targets = keys(m, :target, block)
    logpdf(m, params, transform) + logpdf(m, setdiff(targets, params))
end
"""
    logpdf(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false)

Compute the sum of log-densities for stochastic nodes.

Returns the resulting numeric value of summed log-densities.

* `m`: model containing the stochastic nodes for which to evaluate log-densities.

* `block` : sampling block of stochastic nodes over which to sum densities (default: all stochastic nodes).

* `nodekeys` : nodes over which to sum densities.

* `x` : value (possibly different than the current one) at which to evaluate densities.

* `transform` : whether to evaluate evaluate log-densities of block parameters on the link–transformed scale.
"""
function logpdf(m::Model, nodekeys::Vector{Symbol}, transform::Bool = false)
    lp = 0.0
    for key in nodekeys
        lp += logpdf(m[key], transform)
        isfinite(lp) || break
    end
    #m.likelihood = isnan(lp) ? -Inf : lp
    #m.likelihood
    lp
end

function pseudologpdf(m::Model, nodekeys::Vector{Symbol}, y, transform::Bool = false)
    lp = 0.0
    for key in nodekeys
        lp += pseudologpdf(m[key], y, transform)
        isfinite(lp) || break
    end
    lp
end

function conditional_likelihood(m::Model, nodekeys::Vector{Symbol}, args...)
    conditional_likelihood(m[nodekeys[1]], args...)
end

function rand(m::Model, nodekeys::Vector{Symbol}, x::Int64)
    rand(m[nodekeys[1]], x)
end

"""
    logpdf(m::Model, x::AbstractArray{T}, block::Integer=0,
            transform::Bool=false) where {T<:Real}
"""
function logpdf(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:Real}
    x0 = unlist(m, block)
    lp = logpdf!(m, x, block, transform)
    relist!(m, x0, block)
    lp
end

"""
    logpdf!(m::Model, x::N, block::Integer=0,
              transform::Bool=false) where N<:GeneralNode
"""
function logpdf!(
    m::Model,
    x::N,
    block::Integer = 0,
    transform::Bool = false,
) where {N<:GeneralNode}
    params = keys(m, :block, block)
    targets = keys(m, :target, block)
    m[params] = relist(m, x, params, transform)
    lp = logpdf(m, setdiff(params, targets), transform)
    for key in targets
        isfinite(lp) || break
        node = m[key]
        update!(node, m)
        lp += key in params ? logpdf(node, transform) : logpdf(node)
    end
    lp
end
"""
    gradlogpdf!(m::Model, x::AbstractArray{T}, block::Integer=0,transform::Bool=false)
     where T<:GeneralNode

"""
function gradlogpdf!(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:GeneralNode}
    gradlogpdf!(m, x, block, transform)
end

"""
    gradlogpdf(m::Model, targets::Array{Symbol, 1})::Tuple{Float64, Array{Float64}}

Compute the gradient of log-densities for stochastic nodes.

Returns the resulting gradient vector. Method `gradlogpdf!()` additionally updates model `m` with supplied values `x`.

* `m` : model containing the stochastic nodes for which to compute the gradient.

* `block` : sampling block of stochastic nodes for which to compute the gradient (default: all stochastic nodes).

* `x`: value (possibly different than the current one) at which to compute the gradient.

* `transform`: whether to compute the gradient of block parameters on the link–transformed scale.

* `dtype` : type of differentiation for gradient calculations. Options are
  * `:central` : central differencing.

  * `:forward` : forward differencing.
"""
function gradlogpdf(m::Model, targets::Array{Symbol,1})::Tuple{Float64,Array{Float64}}
    vp = 0.0
    gradp = Array[]
    for key in targets
        node = m[key]
        update!(node, m)
        v, grad = gradlogpdf(node)
        vp += v
        push!(gradp, grad)
    end
    m.likelihood = vp
    gr = .+(gradp...)

    vp, gr
end
"""
    gradlogpdf!(m::Model, x::N, block::Integer=0,transform::Bool=false)::Tuple{Float64, Vector{Float64}}
        where N<:GeneralNode

Returns the resulting gradient vector. Method `gradlogpdf!()` additionally updates model `m` with supplied values `x`.
"""
function gradlogpdf!(
    m::Model,
    x::N,
    block::Integer = 0,
    transform::Bool = false,
)::Tuple{Float64,Vector{Float64}} where {N<:GeneralNode}
    params = keys(m, :block, block)
    targets = keys(m, :target, block)
    # likelihood
    v, grad = gradlogpdf(m, targets)
    # prior
    vp, gradp = gradlogpdf(m[params[1]], x)
    gr = gradp .+ grad
    gr .= any(isnan.(gr)) ? -Inf : gr[:]
    vp + v, gr
end
# """
#     logpdf!(m::Model, x::AbstractArray{T}, block::Integer=0,
#                   transform::Bool=false) where {T<:Real}

# Compute the sum of log-densities for stochastic nodes.

# The resulting numeric value of summed log-densities. Method `logpdf!()` additionally updates model `m` with supplied values `x`.
# """
# function logpdf!(
#     m::Model,
#     x::AbstractArray{T},
#     block::Integer = 0,
#     transform::Bool = false,
# ) where {T<:Real}
#     params = keys(m, :block, block)
#     targets = keys(m, :target, block)
#     m[params] = relist(m, x, params, transform)
#     lp = logpdf(m, setdiff(params, targets), transform)
#     lp = isnan(lp) ? -Inf : lp
#     for key in targets
#         isfinite(lp) || break
#         node = m[key]
#         update!(node, m)
#         lp += key in params ? logpdf(node, transform) : logpdf(node)
#     end
#     lp
# end

function logpdf!(
    m::Model,
    x::AbstractArray{T},
    params::Array{Symbol},
    targets::Array{Symbol},
    transform::Bool = false,
) where {T<:Real}
    m[params] = relist(m, x, params, transform)

    df = [i for i in params if i ∉ targets]
    lp = logpdf(m, df, transform)

    isnan(lp) && return -Inf

    for key in targets
        isfinite(lp) || break
        isnan(lp) && return -Inf
        node = m[key]
        un = update!(node, m)
        #m[key] = node
        m.nodes[key] = un
        lp += key in params ? logpdf(node, transform) : logpdf(node)
    end
    lp
end

function mylogpdf!(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:Real}
    params = keys(m, :block, block)
    targets = keys(m, :target, block)
    diff = setdiff(params, targets)
    myfun(y) = begin
        #v = relist(m, x, params, transform)
        m[params] = relist(m, y, params, transform)
        #lp = sum(y)

        lp = logpdf(m, diff, transform)
        # lp = isnan(lp) ? -Inf : lp
        #lp = 0.0
        for key in targets
            isfinite(lp) || break
            node = m[key]
            update!(node, m)
            lp += key in params ? logpdf(node, transform) : logpdf(node)
        end
        #lp = 0.0
        lp
    end
    lp = myfun(x)

    p = pullback(myfun, x)
    d1 = p[2](1.0)
    lp
end

function rand!(m::Model, x::Int64, block::Integer = 0)
    params = keys(m, :block, block)
    res = rand(m, params, x)
    res
end

"""
    sample!(m::Model, block::Integer=0)

Generate one MCMC sample of values for a specified model.

Returns the model updated with the MCMC sample and, in the case of `block=0`, the `iter` field incremented by 1.

* `m` : model specification.

* `block` : block for which to sample values (default: all blocks).

"""
function sample!(m::Model, block::Integer = 0)
    m.iter += 1
    #isoneblock = block != 0
    #blocks = isoneblock ? block : 1:length(m.samplers)
    #for b in blocks
    for sampler in m.samplers
        #sampler = m.samplers[b]
        #value = sampler.eval(m::Model, b::Int)
        res = sample!(sampler, m)
        if res !== nothing
            m[sampler.params] = res
       #     update!(m, b)
        end
    end
    #m.iter# -= isoneblock
#    m.likelihood = final_likelihood(m)
    m
end

function sample!(s::Sampler, m::Model)
    s.value = unlist(m, s.params, transform=s.transform)
    lpdf(x) = s.tune.logf(m, x, s.params, s.targets, s.transform)
    sample!(s, lpdf, adapt=m.iter < m.burnin, gen=m.iter, model=m)
    relist(m, s.value, s.params, s.transform)
end

function my_relist(m, v, p, s)
    m[p] = relist(m, v, p, s)
end

function pseudologpdf!(
    m::Model,
    x::AbstractArray{T},
    y::AbstractArray,
    params::Vector{Symbol},
    targets::Vector{Symbol},
    transform::Bool = false,
) where {T<:Real}

    m[params] = relist(m, x, params, transform)
    lp = 0.0
    for key in targets
        isfinite(lp) || break
        node = m[key]
        update!(node, m)
        lp += key in params ? pseudologpdf(node, y, transform) : pseudologpdf(node, y)
    end
    lp
end

function conditional_likelihood!(
    m::Model,
    x::AbstractArray{T},
    params::Vector{Symbol},
    targets::Vector{Symbol},
    args...,
) where {T<:Real}
    m[targets] = relist(m, x, targets)
    conditional_likelihood(m, targets, args...)
end

function final_likelihood(model::Model)::Float64
    logpdf(model, keys_output(model))
end
"""
    unlist(m::Model, block::Integer=0, transform::Bool=false)
"""
function unlist(m::Model, block::Integer = 0; transform::Bool = false)
    unlist(m, keys(m, :block, block), transform=transform)
end

"""
    unlist(m::Model, monitoronly::Bool)
"""
function unlist(m::Model, monitoronly::Bool)
    f = let m = m, monitoronly = monitoronly
        key -> begin
            node = m[key]
            lvalue =
                isa(node, TreeVariate) ? unlist_tree(node) : unlist(node)
            monitoronly ? lvalue[node.monitor] : lvalue
        end
    end
    r = vcat(map(f, keys(m, :dependent))..., m.likelihood)
    r = [isa(i, ForwardDiff.Dual) ? ForwardDiff.value(i) : i for i in r]
    r
end
"""
    unlist(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false)

Convert (unlist) sets of logical and/or stochastic node values to vectors.

Returns vectors of concatenated node values.

* `m` : model containing nodes to be unlisted or relisted.

* `block` : sampling block of nodes to be listed (default: all blocks).

* `nodekeys` : node(s) to be listed.

* `transform` : whether to apply a link transformation in the conversion.
"""
function unlist(m::Model, nodekeys::Vector{Symbol}; transform::Bool = false)
    f = let m = m, transform = transform
        key -> unlist(m[key], transform)
    end
    vcat(map(f, nodekeys)...)
end

"""
    relist(m::Model, x::AbstractArray{T}, block::Integer=0,
            transform::Bool=false) where {T<:Real}
"""
function relist(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:Real}
    relist(m, x, keys(m, :block, block), transform)
end
"""
    relist(m::Model, x::AbstractArray{T}, block::Integer=0,
            transform::Bool=false) where {T<:GeneralNode}
"""
function relist(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:GeneralNode}

    relist(m, x, keys(m, :block, block), transform)
end

"""
    relist(m::Model, x::AbstractArray{T},
            nodekeys::Vector{Symbol}, transform::Bool=false) where {T<:Any}
"""
function relist(
    m::Model,
    x::AbstractArray{T},
    nodekeys::Vector{Symbol},
    transform::Bool = false,
) where {T<:Any}
    values = Dict{Symbol,Union{Any,Real}}()

    N = length(x)
    offset = 0
    for key in nodekeys
        value, n = relistlength(m[key], view(x, (offset+1):N), transform)
        values[key] = value
        offset += n
    end
    offset == length(x) ||
        throw(ArgumentError("incompatible number of values to put in nodes"))
    values
end

"""
    relist(m::Model, x::N, nodekeys::Vector{Symbol}, transform::Bool=false) where N<:GeneralNode

Reverse of unlist; ie. Converts vectors to sets of logical and/or stochastic node values. Same inputs and return values as unlist.
"""
function relist(
    m::Model,
    x::N,
    nodekeys::Vector{Symbol},
    transform::Bool = false,
) where {N<:GeneralNode}
    values = Dict{Symbol,Any}()
    offset = 0
    for key in nodekeys
        value, n = relistlength(m[key], x, transform)
        values[key] = value
        offset += n
    end

    values
end

"""
    relist!(m::Model, x::AbstractArray{T}, block::Integer=0,
              transform::Bool=false) where {T<:Any}
"""
function relist!(
    m::Model,
    x::AbstractArray{T},
    block::Integer = 0,
    transform::Bool = false,
) where {T<:Any}
    nodekeys = keys(m, :block, block)
    values = relist(m, x, nodekeys, transform)
    for key in nodekeys
        assign!(m, key, values[key])
    end
    update!(m, block)
end

function assign!(m::Model, key::Symbol, value::T) where {T<:Real}
    m[key].value = value
end

function assign!(m::Model, key::Symbol, value::T) where {T<:GeneralNode}
    m[key].value = value
end

function assign!(m::Model, key::Symbol, value::T) where {T<:Array}
    if isa(m[key], TreeVariate)
        @assert (2,) == size(value)
        d = Array{Float64,1}(undef, 2)
    else
        @assert size(m[key].value) == size(value)
        d = similar(m[key].value)
    end
    for (ind, elem) in enumerate(value)
        d[ind] = elem
    end
    m[key].value = d
end

"""
    relist!(m::Model, x::AbstractArray{T}, nodekey::Symbol,
                  transform::Bool=false) where {T<:Real}

Reverse of unlist; ie. Converts vectors to sets of logical and/or stochastic node values. Same inputs as unlist.

Returns `m`, with values copied to the nodes.

"""
function relist!(
    m::Model,
    x::AbstractArray{T},
    nodekey::Symbol,
    transform::Bool = false,
) where {T<:Real}
    node = m[nodekey]
    m[nodekey] = relist(node, x, transform)
    update!(m, node.targets)
end

"""
    update!(m::Model, block::Integer=0)
"""
function update!(m::Model, block::Integer = 0)
    nodekeys = block == 0 ? keys(m, :dependent) : m.samplers[block].targets
    update!(m, nodekeys)
end
"""
    update!(m::Model, nodekeys::Vector{Symbol})

Update values of logical and stochastic model node according to their relationship with others in a model.

Returns the model with updated nodes.

* `m` : mode with nodes to be updated.

* `block` : sampling block of nodes to be updated (default: all blocks).

* `nodekeys` : nodes to be updated in the given order.
"""
function update!(m::Model, nodekeys::Vector{Symbol})
    for key in nodekeys
        update!(m[key], m)
    end
    m
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
# Autogenerated wrapper script for nv_codec_headers_jll for any
export dynlink_cuda_h, dynlink_cuviddec_h, dynlink_loader_h, dynlink_nvcuvid_h, nvEncodeAPI_h

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `dynlink_cuda_h`
const dynlink_cuda_h_splitpath = ["include", "ffnvcodec", "dynlink_cuda.h"]

# This will be filled out by __init__() for all products, as it must be done at runtime
dynlink_cuda_h_path = ""

# dynlink_cuda_h-specific global declaration
# This will be filled out by __init__()
dynlink_cuda_h = ""

# Relative path to `dynlink_cuviddec_h`
const dynlink_cuviddec_h_splitpath = ["include", "ffnvcodec", "dynlink_cuviddec.h"]

# This will be filled out by __init__() for all products, as it must be done at runtime
dynlink_cuviddec_h_path = ""

# dynlink_cuviddec_h-specific global declaration
# This will be filled out by __init__()
dynlink_cuviddec_h = ""

# Relative path to `dynlink_loader_h`
const dynlink_loader_h_splitpath = ["include", "ffnvcodec", "dynlink_loader.h"]

# This will be filled out by __init__() for all products, as it must be done at runtime
dynlink_loader_h_path = ""

# dynlink_loader_h-specific global declaration
# This will be filled out by __init__()
dynlink_loader_h = ""

# Relative path to `dynlink_nvcuvid_h`
const dynlink_nvcuvid_h_splitpath = ["include", "ffnvcodec", "dynlink_nvcuvid.h"]

# This will be filled out by __init__() for all products, as it must be done at runtime
dynlink_nvcuvid_h_path = ""

# dynlink_nvcuvid_h-specific global declaration
# This will be filled out by __init__()
dynlink_nvcuvid_h = ""

# Relative path to `nvEncodeAPI_h`
const nvEncodeAPI_h_splitpath = ["include", "ffnvcodec", "nvEncodeAPI.h"]

# This will be filled out by __init__() for all products, as it must be done at runtime
nvEncodeAPI_h_path = ""

# nvEncodeAPI_h-specific global declaration
# This will be filled out by __init__()
nvEncodeAPI_h = ""

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"nv_codec_headers")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global dynlink_cuda_h_path = normpath(joinpath(artifact_dir, dynlink_cuda_h_splitpath...))

    global dynlink_cuda_h = dynlink_cuda_h_path
    global dynlink_cuviddec_h_path = normpath(joinpath(artifact_dir, dynlink_cuviddec_h_splitpath...))

    global dynlink_cuviddec_h = dynlink_cuviddec_h_path
    global dynlink_loader_h_path = normpath(joinpath(artifact_dir, dynlink_loader_h_splitpath...))

    global dynlink_loader_h = dynlink_loader_h_path
    global dynlink_nvcuvid_h_path = normpath(joinpath(artifact_dir, dynlink_nvcuvid_h_splitpath...))

    global dynlink_nvcuvid_h = dynlink_nvcuvid_h_path
    global nvEncodeAPI_h_path = normpath(joinpath(artifact_dir, nvEncodeAPI_h_splitpath...))

    global nvEncodeAPI_h = nvEncodeAPI_h_path
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()

ext: jl


text:
module DM

using DataFrames: getproperty
using Format, FileIO, Reexport
using HDF5, DocStringExtensions
import CSV, JSON

@reexport using CSVFiles, DataFrames, Query

include("config.jl")
include("utilities.jl")
include("data_entry/base.jl")
include("data_entry/iterate.jl")
include("data_entry/registry.jl")
include("format_utilities/hdf5_util.jl")
include("format_utilities/csv_util.jl")
include("file_operation.jl")

export DataEntry
export load_config_from_json, convert_recursive_array, parse_string_to_num
export activate_param_set, get_param_set
export save, load, check, load_file_array, delete, writeattr, readattr, read
export walk_entry_file
export load_registry, register_entry, write_registry, query_registry

end # end module


text:
# Copyright Peter Luschny. License is MIT.

module BasicSwingFactorial
export SwingFactorial

"""
Return the factorial of ``n``. Implementation of the swing algorithm using no
primes. An advanced version based on prime-factorization which is much faster
is available as the prime-swing factorial. However the claim is that this is
the fastest algorithm not using prime-factorization. It has the same recursive
structure as his big brother.
"""
function SwingFactorial(n::Int)::BigInt

    smallOddFactorial = BigInt[        0x0000000000000000000000000000001,
    0x0000000000000000000000000000001, 0x0000000000000000000000000000001,
    0x0000000000000000000000000000003, 0x0000000000000000000000000000003,
    0x000000000000000000000000000000f, 0x000000000000000000000000000002d,
    0x000000000000000000000000000013b, 0x000000000000000000000000000013b,
    0x0000000000000000000000000000b13, 0x000000000000000000000000000375f,
    0x0000000000000000000000000026115, 0x000000000000000000000000007233f,
    0x00000000000000000000000005cca33, 0x0000000000000000000000002898765,
    0x00000000000000000000000260eeeeb, 0x00000000000000000000000260eeeeb,
    0x0000000000000000000000286fddd9b, 0x00000000000000000000016beecca73,
    0x000000000000000000001b02b930689, 0x00000000000000000000870d9df20ad,
    0x0000000000000000000b141df4dae31, 0x00000000000000000079dd498567c1b,
    0x00000000000000000af2e19afc5266d, 0x000000000000000020d8a4d0f4f7347,
    0x000000000000000335281867ec241ef, 0x0000000000000029b3093d46fdd5923,
    0x0000000000000465e1f9767cc5866b1, 0x0000000000001ec92dd23d6966aced7,
    0x0000000000037cca30d0f4f0a196e5b, 0x0000000000344fd8dc3e5a1977d7755,
    0x000000000655ab42ab8ce915831734b, 0x000000000655ab42ab8ce915831734b,
    0x00000000d10b13981d2a0bc5e5fdcab, 0x0000000de1bc4d19efcac82445da75b,
    0x000001e5dcbe8a8bc8b95cf58cde171, 0x00001114c2b2deea0e8444a1f3cecf9,
    0x0002780023da37d4191deb683ce3ffd, 0x002ee802a93224bddd3878bc84ebfc7,
    0x07255867c6a398ecb39a64b83ff3751, 0x23baba06e131fc9f8203f7993fc1495]

    function oddProduct(m, len)
        if len < 24
            p = BigInt(m)
            for k in 2:2:2(len-1)
                p *= (m - k)
            end
            return p
        end
        hlen = len >> 1
        oddProduct(m - 2 * hlen, len - hlen) * oddProduct(m, hlen)
    end

    function oddFactorial(n)
        if n < 41
            oddFact = smallOddFactorial[1+n]
            sqrOddFact = smallOddFactorial[1+div(n, 2)]
        else
            sqrOddFact, oldOddFact = oddFactorial(div(n, 2))
            len = div(n - 1, 4)
            (n % 4) != 2 && (len += 1)
            high = n - ((n + 1) & 1)
            oddSwing = div(oddProduct(high, len), oldOddFact)
            oddFact = sqrOddFact^2 * oddSwing
        end
        (oddFact, sqrOddFact)
    end

    n < 0 && ArgumentError("n must be ≥ 0")
    if n == 0 return 1 end
    sh = n - count_ones(n)
    oddFactorial(n)[1] << sh
end

#START-TEST-################################################
using Test

function main()
    @testset "SwingFactorial" begin
        for n in 0:999
            S = SwingFactorial(n)
            B = Base.factorial(BigInt(n))
            @test S == B
        end
    end

    GC.gc()
    n = 1000000
    @time SwingFactorial(n)
end

main()

end # module

max_stars_repo_name: NaldoDj/Fast-Factorial-Functions

max_issues_repo_name: NaldoDj/Fast-Factorial-Functions

max_forks_repo_name: NaldoDj/Fast-Factorial-Functions


text:
# https://juliadocs.github.io/Documenter.jl/stable/man/guide/#Package-Guide
# push!(LOAD_PATH,"../src/")

# Run these locally to build docs/build folder:
# PS C:\Users\akua\repos\github\ModelBaseEcon.jl> julia --color=yes --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
# PS C:\Users\akua\repos\github\ModelBaseEcon.jl> julia --project=docs/ docs/make.jl

using Documenter, ModelBaseEcon

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(sitename = "ModelBaseEcon.jl",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         modules = [ModelBaseEcon],
         doctest = false,
         pages = [
        "Home" => "index.md",
        "Examples" => "examples.md"
    ]
)

deploydocs(
    repo = "github.com/bankofcanada/ModelBaseEcon.jl.git",
)


text:
function gen_hamiltonian(omegaQ,omegaR,g,cutoffN, strong_interaction, Braak)
  # This function generates a times invariant matrix for the Jaynes-Cummings
  # Hamiltonian. The parameters omegaR, omegaQ and g should be given in
  # rad/s. The cutoffN parameter is the highest energy level of the harmonic
  # oscillator that is considered. Higher levels are truncated.

  # Some constant parameters for the Jaynes-Cummings
  # Some constant parameters for the Jaynes-Cummings
  identity = [1 0 ; 0 1]; sigmaZ = 0.5*[1 0 ; 0 -1];
  sigmaMinus = [0 1 ; 0 0]; sigmaPlus = [0 0 ; 1 0];
  global h_bar
	#info(a_dagger(cutoffN) * a(cutoffN))
	half_matrix = Braak ? zeros(cutoffN, cutoffN) : eye(cutoffN)*0.5
  # The Jaynes-Cummings Hamiltonian is written with the following
  # conventions: omegaR is the resonator frequency, omegaQ is the qubit
  # frequency and g is the coupling.
  jaynesCummings = zeros(Complex128,2*cutoffN,2*cutoffN)
  jaynesCummings = jaynesCummings +
     kron(h_bar .* omegaR .* ((a_dagger(cutoffN) * a(cutoffN))+half_matrix), identity) -	# harmonic oscillator term
     h_bar .*omegaQ .* kron(eye(cutoffN), sigmaZ) +						# qubit term
     h_bar .* g .* ( kron(a(cutoffN), sigmaPlus) + kron(a_dagger(cutoffN), sigmaMinus) )	# interaction term
	if Braak
		jaynesCummings-=h_bar .*omegaQ .* kron(eye(cutoffN), sigmaZ)
	end
  if strong_interaction
  	jaynesCummings+= h_bar .* g .* (kron(a(cutoffN), sigmaMinus) + kron(a_dagger(cutoffN), sigmaPlus))# strong interaction term
 	end
  return jaynesCummings
end

#This function returns the desired B block
function get_B_block(hamiltonian, index)
	low = (index+1)*2
	return sub(hamiltonian, low:low+1, low:low+1)
end


text:
module TensorNetworkCodes

# imports
using Combinatorics: combinations
using GraphPlot: gplot
using ITensors: Index, ITensor, dim, hastags, inds, tags
using LightGraphs: Graph, add_edge!
using Random: AbstractRNG, MersenneTwister, RandomDevice, GLOBAL_RNG, rand
using StatsBase: ProbabilityWeights, sample

# exports
export QuantumCode, SimpleCode, CodeGraph, TensorNetworkCode
include("types.jl")
export pauli_are_commuting, pauli_are_independent, pauli_commutation
export pauli_pow, pauli_product, pauli_product_pow
export pauli_random_operator, pauli_rep_change, pauli_weight
include("pauli_functions.jl")
export num_qubits, verify
export find_distance_logicals, find_pure_error, find_pure_errors, find_syndrome
export gauge, permute, purify
include("code_functions.jl")
export edges, new_indices, nodes, num_nodes, physical_neighbours, set_coords!, shift_coords!
export coords, edge_indices, edge_types, node_indices, node_types
export set_coords!, set_edge_indices!, set_edge_types!, set_node_indices!, set_node_types!
include("code_graph_functions.jl")
export combine, contract, contract_by_coords, fusion
include("contraction_functions.jl")
export plot_code, plot_operator
include("plotting_functions.jl")
export five_qubit_code, five_qubit_surface_code, steane_code
export random_code, random_stabilizer_state
include("examples/simple.jl")
export almost_rotated_surface_code, rotated_surface_code, surface_code
include("examples/surface.jl")
export code_to_Itensor, identity_coset, all_cosets
export physical_tensor, create_virtual_tensor
include("itensors_functions.jl")

# include submodules (not exported)
include("TNDecode.jl")
include("TNDistance.jl")
include("QecsimAdaptors.jl")

end


text:
@dict_readable struct RunInTerminalRequestArguments <: Outbound
    kind::Union{Missing,String}
    title::Union{Missing,String}
    cwd::String
    args::Vector{String}
    env::Union{Missing,Dict{String,Union{Missing,String}}}
end

@dict_readable struct InitializeRequestArguments <: Outbound
    clientID::Union{Missing,String}
    clientName::Union{Missing,String}
    adapterID::String
    locale::Union{Missing,String}
    linesStartAt1::Union{Missing,Bool}
    columnsStartAt1::Union{Missing,Bool}
    pathFormat::Union{Missing,String}
    supportsVariableType::Union{Missing,Bool}
    supportsVariablePaging::Union{Missing,Bool}
    supportsRunInTerminalRequest::Union{Missing,Bool}
    supportsMemoryReferences::Union{Missing,Bool}
end

@dict_readable struct ConfigurationDoneArguments <: Outbound
end

@dict_readable struct LaunchArguments <: Outbound
    noDebug::Union{Missing,Bool}
    __restart::Any
end

@dict_readable struct AttachArguments <: Outbound
    __restart::Any
end

@dict_readable struct RestartArguments <: Outbound
end

@dict_readable struct DisconnectArguments <: Outbound
    restart::Union{Missing,Bool}
    terminateDebuggee::Union{Missing,Bool}
end

@dict_readable struct TerminateArguments <: Outbound
    restart::Union{Missing,Bool}
end

@dict_readable struct SetBreakpointsArguments <: Outbound
    source::Source
    breakpoints::Union{Missing,Vector{SourceBreakpoint}}
    lines::Union{Missing,Vector{Int}}
    sourceModified::Union{Missing,Bool}
end

@dict_readable struct SetFunctionBreakpointsArguments <: Outbound
    breakpoints::Vector{FunctionBreakpoint}
end

@dict_readable struct SetExceptionBreakpointsArguments <: Outbound
    filters::Vector{String}
    exceptionOptions::Union{Missing,Vector{ExceptionOptions}}
end

@dict_readable struct DataBreakpointInfoArguments <: Outbound
    variableReference::Union{Missing,Int}
    name::String
end

@dict_readable struct SetDataBreakpointsArguments <: Outbound
    breakpoints::Vector{DataBreakpoint}
end

@dict_readable struct ContinueArguments <: Outbound
    threadId::Int
end

@dict_readable struct NextArguments <: Outbound
    threadId::Int
end

@dict_readable struct StepInArguments <: Outbound
    threadId::Int
    targetId::Union{Missing,Int}
end

@dict_readable struct StepOutArguments <: Outbound
    threadId::Int
end

@dict_readable struct StepBackArguments <: Outbound
    threadId::Int
end

@dict_readable struct ReverseContinueArguments <: Outbound
    threadId::Int
end

@dict_readable struct RestartFrameArguments <: Outbound
    frameId::Int
end

@dict_readable struct GotoArguments <: Outbound
    threadId::Int
    targetId::Int
end

@dict_readable struct PauseArguments <: Outbound
    threadId::Int
end

@dict_readable struct StackTraceArguments <: Outbound
    threadId::Int
    startFrame::Union{Missing,Int}
    levels::Union{Missing,Int}
    format::Union{Missing,StackFrameFormat}
end

@dict_readable struct ScopesArguments <: Outbound
    frameId::Int
end

@dict_readable struct VariablesArguments <: Outbound
    variablesReference::Int
    filter::Union{Missing,String}
    start::Union{Missing,Int}
    count::Union{Missing,Int}
    format::Union{Missing,ValueFormat}
end

@dict_readable struct SetVariableArguments <: Outbound
    variablesReference::Int
    name::String
    value::String
    format::Union{Missing,ValueFormat}
end

@dict_readable struct SourceArguments <: Outbound
    source::Union{Missing,Source}
    sourceReference::Int
end

@dict_readable struct TerminateThreadsArguments <: Outbound
    threadIds::Union{Missing,Vector{Int}}
end

@dict_readable struct ModulesArguments <: Outbound
    startModule::Union{Missing,Int}
    moduleCount::Union{Missing,Int}
end

@dict_readable struct LoadedSourcesArguments <: Outbound
end

@dict_readable struct EvaluateArguments <: Outbound
    expression::String
    frameId::Union{Missing,Int}
    context::Union{Missing,String}
    format::Union{Missing,ValueFormat}
end

@dict_readable struct SetExpressionArguments <: Outbound
    expression::String
    value::String
    frameId::Union{Missing,Int}
    format::Union{Missing,ValueFormat}
end

@dict_readable struct StepInTargetsArguments <: Outbound
    frameId::Int
end

@dict_readable struct GotoTargetsArguments <: Outbound
    source::Source
    line::Int
    column::Union{Missing,Int}
end

@dict_readable struct CompletionsArguments <: Outbound
    frameId::Union{Missing,Int}
    text::String
    column::Int
    line::Union{Missing,Int}
end

@dict_readable struct ExceptionInfoArguments <: Outbound
    threadId::Int
end

@dict_readable struct ReadMemoryArguments <: Outbound
    memoryReference::String
    offset::Union{Missing,Int}
    count::Int
end

@dict_readable struct DisassembleArguments <: Outbound
    memoryReference::String
    offset::Union{Missing,Int}
    instructionOffset::Union{Missing,Int}
    instructionCount::Int
    resolveSymbols::Union{Missing,Bool}
end

@dict_readable struct BreakpointLocationsArguments <: Outbound
    source::Source
    line::Int
    column::Union{Missing,Int}
    endLine::Union{Missing,Int}
    endColumn::Union{Missing,Int}
end

@dict_readable struct BreakpointLocationsResponseArguments <: Outbound
    breakpoints::Vector{BreakpointLocation}
end


text:
module KeyTuples

import DataFrames: DataFrame
import NullableArrays: NullableArray

export KeyTuple, keytuples_to_df

"""
Tuple-like type with keys for accessing elements. Valid keys are
symbols. Use the constructor
```julia
KeyTuple(:a => 1, :b => 2)
```
to create objects.
"""
immutable KeyTuple{K,T <: Tuple} <: Associative
    values::T
    function KeyTuple(values::T)
        isa(K, Tuple{Vararg{Symbol}}) ||
            throw(ArgumentError("Keys need to be a tuple of symbols."))
        length(K) == length(values) ||
            throw(ArgumentError("Length mismatch between keys and values."))
        new{K,T}(values)
    end
end

KeyTuple{T <: Tuple}(keys::Tuple{Vararg{Symbol}}, values::T) =
    KeyTuple{keys,T}(values)

function KeyTuple(key_value_pairs::Pair...)
    keys = tuple(map(first, key_value_pairs)...)
    values = tuple(map(last, key_value_pairs)...)
    KeyTuple{keys,typeof(values)}(values)
end

function Base.show{K,T}(io::IO, x::KeyTuple{K,T})
    print(io, "KeyTuple")
    show(io, tuple([Pair(key, value) for (key, value) in zip(K, x.values)]...))
end

Base.keys{K,T}(x::KeyTuple{K,T}) = K

Base.values(x::KeyTuple) = x.values

Base.getindex{K,T}(x::KeyTuple{K,T}, index) = getindex(x.values, index)

function Base.getindex{K,T}(x::KeyTuple{K,T}, key::Symbol)
    i = findfirst(K, key)
    i == 0 && throw(BoundsError("no key $(key)"))
    x.values[i]
end

Base.length{K,T}(::KeyTuple{K,T}) = length(K)

Base.start(x::KeyTuple) = 1

Base.done{K,T}(x::KeyTuple{K,T}, iter) = iter > length(K)

Base.next{K,T}(x::KeyTuple{K,T}, iter) = (Pair(K[iter], x.values[iter]), iter + 1)

Base.endof(x::KeyTuple) = length(x)

Base.first{K,T}(x::KeyTuple{K,T}) = Pair(K[1], x.values[1])

Base.last{K,T}(x::KeyTuple{K,T}) = Pair(K[end], x.values[end])

function Base.:(==){Kx,Tx,Ky,Ty}(x::KeyTuple{Kx,Tx}, y::KeyTuple{Ky,Ty})
    (Kx == Ky) && (x.values == y.values)
end

Base.hash{K,T}(x::KeyTuple{K,T}, h::UInt) = hash(x.values, hash(K, hash(:KeyTuple, h)))

"""
Internal helper function for conversion to DataFrame. Return a vector
of `key => column` pairs.
"""
function _keytuples_to_df{K,T}(vector::Vector{KeyTuple{K,T}})
    (Pair(key, map(x->x.values[index], vector)) for (index,key) in enumerate(K))
end

_keytuples_to_df(itr) = _keytuples_to_df(collect(itr)) # FIXME use iterators

"""
Convert an iterable that returns conformable `KeyTuple`s to a
`DataFrame` with the given columns.
"""
keytuples_to_df(itr) = DataFrame(; _keytuples_to_df(itr)...)

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

#=
Möller, T. & Trumbore, B., 1997.
(https://www.tandfonline.com/doi/abs/10.1080/10867651.1997.10487468)
=#
function intersection(f, r::Ray{3,T}, t::Triangle{3,T}) where {T}
  vs = vertices(t)
  o = origin(r)
  d = direction(r)

  e₁ = vs[3] - vs[1]
  e₂ = vs[2] - vs[1]
  p = d × e₂
  det = e₁ ⋅ p

  # keep det > 0, modify T accordingly
  if det > atol(T)
    τ = o - vs[1]
  else
    τ = vs[1] - o
    det = -det
  end

  if det < atol(T)
    # This ray is parallel to the plane of the triangle.
    return @IT NoIntersection nothing f
  end

  # calculate u parameter and test bounds
  u = τ ⋅ p
  if u < -atol(T) || u > det
    return @IT NoIntersection nothing f
  end

  q = τ × e₁

  # calculate v parameter and test bounds
  v = d ⋅ q
  if v < -atol(T) || u + v > det
    return @IT NoIntersection nothing f
  end

  λ = (e₂ ⋅ q) * (one(T) / det)

  if λ < -atol(T)
    return @IT NoIntersection nothing f
  end

  λ = clamp(λ, zero(T), typemax(T))

  return @IT IntersectingRayTriangle r(λ) f
end


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAM3.png")

# -----------------------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->sin(11t),t->250,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM4.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->sin(11t),t->250,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM5.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->sin(11t),t->250,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAM6.png")

#------------------------------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->cos(t)+exp(-t),t->100,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM7.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->cos(t)+exp(-t),t->100,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0)
png(path*"IS_exAM8.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->cos(t)+exp(-t),t->100,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=-1.0:0.0005:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAM9.png")


text:
using GLMakie, Statistics, Printf

"""
visualize(states::AbstractArray; statenames = string.(1:length(states)), quantiles = (0.1, 0.99), aspect = (1,1,1), resolution = (1920, 1080), statistics = false, title = "Field = ")
# Description
Visualize 3D states
# Arguments
- `states`: Array{Array{Float64,3},1}. An array of arrays containing different fields
# Keyword Arguments
- `statenames`: Array{String,1}. An array of stringnames
- `aspect`: Tuple{Int64,Int64,Float64}. Determines aspect ratio of box for volumes
- `resolution`: Resolution of preliminary makie window
- `statistics`: boolean. toggle for displaying statistics
# Return
- `scene`: Scene. A preliminary scene object for manipulation
"""
function visualize(
    states::AbstractArray;
    statenames = string.(1:length(states)),
    units = ["" for i in eachindex(states)],
    aspect = (1, 1, 1),
    resolution = (1920, 1080),
    statistics = false,
    title = "Field = ",
    bins = 300,
)
    # Create scene
    scene, layout = layoutscene(resolution = resolution)
    lscene = layout[2:4, 2:4] = LScene(scene)
    width = round(Int, resolution[1] / 4) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    if statistics
        llscene =
            layout[4, 1] = Axis(
                scene,
                xlabel = @lift(statenames[$statenode] * units[$statenode]),
                xlabelcolor = :black,
                ylabel = "pdf",
                ylabelcolor = :black,
                xlabelsize = 40,
                ylabelsize = 40,
                xticklabelsize = 25,
                yticklabelsize = 25,
                xtickcolor = :black,
                ytickcolor = :black,
                xticklabelcolor = :black,
                yticklabelcolor = :black,
            )
        layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)
    end

    # x,y,z are for determining the aspect ratio of the box
    if (typeof(aspect) <: Tuple) & (length(aspect) == 3)
        x, y, z = aspect
    else
        x, y, z = size(states[1])
    end

    # Clim sliders
    upperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    clims = @lift((
        quantile($state[:], $lowerclim_node),
        quantile($state[:], $upperclim_node),
    ))
    cmap_rgb = @lift(to_colormap($colornode))
    titlename = @lift(title * $statename) # use padding and appropriate centering

    # Statistics
    if statistics
        histogram_node = @lift(histogram($state, bins = bins))
        xs = @lift($histogram_node[1])
        ys = @lift($histogram_node[2])
        pdf = GLMakie.AbstractPlotting.barplot!(
            llscene,
            xs,
            ys,
            color = :red,
            strokecolor = :red,
            strokewidth = 1,
        )
        @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
        @lift(GLMakie.AbstractPlotting.ylims!(
            llscene,
            extrema($histogram_node[2]),
        ))
        vlines!(
            llscene,
            @lift($clims[1]),
            color = :black,
            linewidth = width / 100,
        )
        vlines!(
            llscene,
            @lift($clims[2]),
            color = :black,
            linewidth = width / 100,
        )
    end

    # Volume Plot
    volume!(
        lscene,
        0..x,
        0..y,
        0..z,
        state,
        camera = cam3d!,
        colormap = cmap_rgb,
        colorrange = clims,
    )
    # Camera
    cam = cameracontrols(scene.children[1])
    eyeposition = Float32[2, 2, 1.3]
    lookat = Float32[0.82, 0.82, 0.1]
    # Title
    supertitle =
        layout[1, 2:4] = Label(scene, titlename, textsize = 50, color = :black)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end
    lowerclim_string = @lift(
        "lower clim quantile = " *
        @sprintf("%0.2f", $lowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[1])
    )
    upperclim_string = @lift(
        "upper clim quantile = " *
        @sprintf("%0.2f", $upperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[2])
    )
    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )
    layout[1, 1] = Label(scene, "Menu", width = width, textsize = 50)

    # Modify Axis
    axis = scene.children[1][OldAxis]
    # axis[:names][:axisnames] = ("↓ Zonal [m] ", "Meriodonal [m]↓ ", "Depth [m]↓ ")
    axis[:names][:axisnames] = ("↓", "↓ ", "↓ ")
    axis[:names][:align] =
        ((:left, :center), (:right, :center), (:right, :center))
    # need to adjust size of ticks first and then size of axis names
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)
    # axis[:ticks][:ranges_labels].val # current axis labels
    xticks = collect(range(-0, aspect[1], length = 2))
    yticks = collect(range(-0, aspect[2], length = 6))
    zticks = collect(range(-0, aspect[3], length = 2))
    ticks = (xticks, yticks, zticks)
    axis[:ticks][:ranges] = ticks
    xtickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    xtickslabels[end] = "1e6"
    ytickslabels = ["", "south", "", "", "north", ""]
    ztickslabels = [@sprintf("%0.1f", (xtick)) for xtick in xticks]
    labels = (xtickslabels, ytickslabels, ztickslabels)
    axis[:ticks][:labels] = labels

    display(scene)
    # Change the default camera position after the fact
    # note that these change dynamically as the plot is manipulated
    return scene
end

"""
histogram(array; bins = 100)
# Description
return arrays for plotting histogram
"""
function histogram(
    array;
    bins = minimum([100, length(array)]),
    normalize = true,
)
    tmp = zeros(bins)
    down, up = extrema(array)
    down, up = down == up ? (down - 1, up + 1) : (down, up) # edge case
    bucket = collect(range(down, up, length = bins + 1))
    normalization = normalize ? length(array) : 1
    for i in eachindex(array)
        # normalize then multiply by bins
        val = (array[i] - down) / (up - down) * bins
        ind = ceil(Int, val)
        # handle edge cases
        ind = maximum([ind, 1])
        ind = minimum([ind, bins])
        tmp[ind] += 1 / normalization
    end
    return (bucket[2:end] + bucket[1:(end - 1)]) .* 0.5, tmp
end

# 2D visualization
function visualize(
    states::Array{Array{S, 2}, 1};
    statenames = string.(1:length(states)),
    units = ["" for i in eachindex(states)],
    aspect = (1, 1, 1),
    resolution = (2412, 1158),
    title = "Zonal and Temporal Average of ",
    xlims = (0, 1),
    ylims = (0, 1),
    bins = 300,
) where {S}
    # Create scene
    scene, layout = layoutscene(resolution = resolution)
    lscene =
        layout[2:4, 2:4] = Axis(
            scene,
            xlabel = "South to North [m]",
            xlabelcolor = :black,
            ylabel = "Depth [m]",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = 25,
            yticklabelsize = 25,
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
            titlesize = 50,
        )
    width = round(Int, resolution[1] / 4) # make menu 1/4 of preliminary resolution

    # Create choices and nodes
    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    interpolationlabels = ["contour", "heatmap"]
    interpolationchoices = [true, false]
    interpolationnode = Node(interpolationchoices[1])

    # Statistics
    llscene =
        layout[4, 1] = Axis(
            scene,
            xlabel = @lift(statenames[$statenode] * " " * units[$statenode]),
            xlabelcolor = :black,
            ylabel = "pdf",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = 25,
            yticklabelsize = 25,
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
        )
    layout[3, 1] = Label(scene, "Statistics", width = width, textsize = 50)

    # Clim sliders
    upperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    #ylims = @lift(range($lowerval, $upperval, length = $))
    # Lift Nodes
    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    unit = @lift(units[$statenode])
    oclims = @lift((
        quantile($state[:], $lowerclim_node),
        quantile($state[:], $upperclim_node),
    ))
    cmap_rgb = @lift(
        $oclims[1] < $oclims[2] ? to_colormap($colornode) :
        reverse(to_colormap($colornode))
    )
    clims = @lift(
        $oclims[1] != $oclims[2] ? (minimum($oclims), maximum($oclims)) :
        (minimum($oclims) - 1, maximum($oclims) + 1)
    )
    xlims = Array(range(xlims[1], xlims[2], length = 4)) #collect(range(xlims[1], xlims[2], length = size(states[1])[1]))
    ylims = Array(range(ylims[1], ylims[2], length = 4)) #@lift(collect(range($lowerval], $upperval, length = size($state)[2])))
    # newrange = @lift(range($lowerval, $upperval, length = 4))
    # lscene.yticks = @lift(Array($newrange))
    titlename = @lift(title * $statename) # use padding and appropriate centering
    layout[1, 2:4] = Label(scene, titlename, textsize = 50)
    # heatmap
    heatmap1 = heatmap!(
        lscene,
        xlims,
        ylims,
        state,
        interpolate = interpolationnode,
        colormap = cmap_rgb,
        colorrange = clims,
    )

    # statistics
    histogram_node = @lift(histogram($state, bins = bins))
    xs = @lift($histogram_node[1])
    ys = @lift($histogram_node[2])
    pdf = GLMakie.AbstractPlotting.barplot!(
        llscene,
        xs,
        ys,
        color = :red,
        strokecolor = :red,
        strokewidth = 1,
    )
    @lift(GLMakie.AbstractPlotting.xlims!(llscene, extrema($state)))
    @lift(GLMakie.AbstractPlotting.ylims!(llscene, extrema($histogram_node[2])))
    vlines!(llscene, @lift($clims[1]), color = :black, linewidth = width / 100)
    vlines!(llscene, @lift($clims[2]), color = :black, linewidth = width / 100)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end

    interpolationmenu =
        Menu(scene, options = zip(interpolationlabels, interpolationchoices))
    on(interpolationmenu.selection) do s
        interpolationnode[] = s
        heatmap1 = heatmap!(
            lscene,
            xlims,
            ylims,
            state,
            interpolate = s,
            colormap = cmap_rgb,
            colorrange = clims,
        )
    end

    newlabel = @lift($statename * " " * $unit)
    cbar = Colorbar(scene, heatmap1, label = newlabel)
    cbar.width = Relative(1 / 3)
    cbar.height = Relative(5 / 6)
    cbar.halign = :center
    cbar.flipaxisposition = true
    # cbar.labelpadding = -350
    cbar.labelsize = 50

    lowerclim_string = @lift(
        "clim quantile = " *
        @sprintf("%0.2f", $lowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[1])
    )
    upperclim_string = @lift(
        "clim quantile = " *
        @sprintf("%0.2f", $upperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[2])
    )

    # depends on makie version, vbox for old, vgrid for new
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(
            scene,
            "plotting options",
            width = width,
            textsize = 30,
            padding = (0, 0, 10, 0),
        ),
        interpolationmenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )

    layout[2:4, 5] = vgrid!(
        Label(
            scene,
            "Color Bar",
            width = width / 2,
            textsize = 50,
            padding = (25, 0, 0, 00),
        ),
        cbar,
    )
    layout[1, 1] = Label(scene, "Menu", width = width, textsize = 50)
    display(scene)
    return scene
end

function volumeslice(
    states::AbstractArray;
    statenames = string.(1:length(states)),
    units = ["" for i in eachindex(states)],
    aspect = (1, 1, 32 / 192),
    resolution = (2678, 1030),
    statistics = false,
    title = "Volume plot of ",
    bins = 300,
    statlabelsize = (20, 20),
)
    scene, layout = layoutscene(resolution = resolution)
    volumescene = layout[2:4, 2:4] = LScene(scene)
    menuwidth = round(Int, 350)
    layout[1, 1] = Label(scene, "Menu", width = menuwidth, textsize = 50)

    slice_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.0)
    slice_node = slice_slider.value

    directionindex = [1, 2, 3]
    directionnames = ["x-slice", "y-slice", "z-slice"]
    directionnode = Node(directionindex[1])

    stateindex = collect(1:length(states))
    statenode = Node(stateindex[1])

    layout[1, 2:4] =
        Label(scene, @lift(title * statenames[$statenode]), textsize = 50)

    colorchoices = [:balance, :thermal, :dense, :deep, :curl, :thermometer]
    colornode = Node(colorchoices[1])

    state = @lift(states[$statenode])
    statename = @lift(statenames[$statenode])
    unit = @lift(units[$statenode])
    nx = @lift(size($state)[1])
    ny = @lift(size($state)[2])
    nz = @lift(size($state)[3])
    nr = @lift([$nx, $ny, $nz])

    nslider = 100
    xrange = range(0.00, aspect[1], length = nslider)
    yrange = range(0.00, aspect[2], length = nslider)
    zrange = range(0.00, aspect[3], length = nslider)
    constx = collect(reshape(xrange, (nslider, 1, 1)))
    consty = collect(reshape(yrange, (1, nslider, 1)))
    constz = collect(reshape(zrange, (1, 1, nslider)))
    matx = zeros(nslider, nslider, nslider)
    maty = zeros(nslider, nslider, nslider)
    matz = zeros(nslider, nslider, nslider)
    matx .= constx
    maty .= consty
    matz .= constz
    sliceconst = [matx, maty, matz]
    planeslice = @lift(sliceconst[$directionnode])

    upperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    upperclim_node = upperclim_slider.value
    lowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    lowerclim_node = lowerclim_slider.value

    clims = @lift((
        quantile($state[:], $lowerclim_node),
        quantile($state[:], $upperclim_node),
    ))

    volume!(
        volumescene,
        0..aspect[1],
        0..aspect[2],
        0..aspect[3],
        state,
        overdraw = false,
        colorrange = clims,
        colormap = @lift(to_colormap($colornode)),
    )

    alpha_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.5)
    alphanode = alpha_slider.value

    slicecolormap = @lift(cgrad(:viridis, alpha = $alphanode))
    v = volume!(
        volumescene,
        0..aspect[1],
        0..aspect[2],
        0..aspect[3],
        planeslice,
        algorithm = :iso,
        isorange = 0.005,
        isovalue = @lift($slice_node * aspect[$directionnode]),
        transparency = true,
        overdraw = false,
        visible = true,
        colormap = slicecolormap,
        colorrange = [-1, 0],
    )

    # Volume histogram

    layout[3, 1] = Label(scene, "Statistics", textsize = 50)
    hscene =
        layout[4, 1] = Axis(
            scene,
            xlabel = @lift(statenames[$statenode] * " " * units[$statenode]),
            xlabelcolor = :black,
            ylabel = "pdf",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = statlabelsize[1],
            yticklabelsize = statlabelsize[2],
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
        )

    histogram_node = @lift(histogram($state, bins = bins))
    vxs = @lift($histogram_node[1])
    vys = @lift($histogram_node[2])
    pdf = GLMakie.AbstractPlotting.barplot!(
        hscene,
        vxs,
        vys,
        color = :red,
        strokecolor = :red,
        strokewidth = 1,
    )

    @lift(GLMakie.AbstractPlotting.xlims!(hscene, extrema($vxs)))
    @lift(GLMakie.AbstractPlotting.ylims!(hscene, extrema($vys)))
    vlines!(
        hscene,
        @lift($clims[1]),
        color = :black,
        linewidth = menuwidth / 100,
    )
    vlines!(
        hscene,
        @lift($clims[2]),
        color = :black,
        linewidth = menuwidth / 100,
    )

    # Slice
    sliceupperclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
    sliceupperclim_node = sliceupperclim_slider.value
    slicelowerclim_slider =
        Slider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
    slicelowerclim_node = slicelowerclim_slider.value

    slicexaxislabel = @lift(["y", "x", "x"][$directionnode])
    sliceyaxislabel = @lift(["z", "z", "y"][$directionnode])

    slicexaxis = @lift([[1, $ny], [1, $nx], [1, $nx]][$directionnode])
    sliceyaxis = @lift([[1, $nz], [1, $nz], [1, $ny]][$directionnode])

    slicescene =
        layout[2:4, 5:6] =
            Axis(scene, xlabel = slicexaxislabel, ylabel = sliceyaxislabel)

    sliced_state1 = @lift( $state[
        round(Int, 1 + $slice_node * (size($state)[1] - 1)),
        1:size($state)[2],
        1:size($state)[3],
    ])
    sliced_state2 = @lift( $state[
        1:size($state)[1],
        round(Int, 1 + $slice_node * (size($state)[2] - 1)),
        1:size($state)[3],
    ])
    sliced_state3 = @lift( $state[
        1:size($state)[1],
        1:size($state)[2],
        round(Int, 1 + $slice_node * (size($state)[3] - 1)),
    ])
    sliced_states = @lift([$sliced_state1, $sliced_state2, $sliced_state3])
    sliced_state = @lift($sliced_states[$directionnode])

    oclims = @lift((
        quantile($sliced_state[:], $slicelowerclim_node),
        quantile($sliced_state[:], $sliceupperclim_node),
    ))
    slicecolormapnode = @lift($oclims[1] < $oclims[2] ? $colornode : $colornode)
    sliceclims = @lift(
        $oclims[1] != $oclims[2] ? (minimum($oclims), maximum($oclims)) :
        (minimum($oclims) - 1, maximum($oclims) + 1)
    )

    heatmap1 = heatmap!(
        slicescene,
        slicexaxis,
        sliceyaxis,
        sliced_state,
        interpolate = true,
        colormap = slicecolormapnode,
        colorrange = sliceclims,
    )

    # Colorbar
    newlabel = @lift($statename * " " * $unit)
    cbar = Colorbar(scene, heatmap1, label = newlabel)
    cbar.width = Relative(1 / 3)
    # cbar.height = Relative(5/6)
    cbar.halign = :left
    # cbar.flipaxisposition = true
    # cbar.labelpadding = -250
    cbar.labelsize = 50

    @lift(GLMakie.AbstractPlotting.xlims!(slicescene, extrema($slicexaxis)))
    @lift(GLMakie.AbstractPlotting.ylims!(slicescene, extrema($sliceyaxis)))

    sliceindex = @lift([
        round(Int, 1 + $slice_node * ($nx - 1)),
        round(Int, 1 + $slice_node * ($ny - 1)),
        round(Int, 1 + $slice_node * ($nz - 1)),
    ][$directionnode])
    slicestring =
        @lift(directionnames[$directionnode] * " of " * statenames[$statenode])
    layout[1, 5:6] = Label(scene, slicestring, textsize = 50)

    axis = scene.children[1][OldAxis]
    axis[:names][:axisnames] = ("↓", "↓ ", "↓ ")
    axis[:names][:align] =
        ((:left, :center), (:right, :center), (:right, :center))
    axis[:names][:textsize] = (50.0, 50.0, 50.0)
    axis[:ticks][:textsize] = (00.0, 00.0, 00.0)

    # Menus
    statemenu = Menu(scene, options = zip(statenames, stateindex))
    on(statemenu.selection) do s
        statenode[] = s
    end

    colormenu = Menu(scene, options = zip(colorchoices, colorchoices))
    on(colormenu.selection) do s
        colornode[] = s
    end

    # Slice Statistics
    layout[1, 7] = Label(scene, "Slice Menu", width = menuwidth, textsize = 50)
    layout[3, 7] = Label(scene, "Slice Statistics", textsize = 50)
    hslicescene =
        layout[4, 7] = Axis(
            scene,
            xlabel = @lift(statenames[$statenode] * " " * units[$statenode]),
            xlabelcolor = :black,
            ylabel = "pdf",
            ylabelcolor = :black,
            xlabelsize = 40,
            ylabelsize = 40,
            xticklabelsize = statlabelsize[1],
            yticklabelsize = statlabelsize[2],
            xtickcolor = :black,
            ytickcolor = :black,
            xticklabelcolor = :black,
            yticklabelcolor = :black,
        )

    slicehistogram_node = @lift(histogram($sliced_state, bins = bins))
    xs = @lift($slicehistogram_node[1])
    ys = @lift($slicehistogram_node[2])
    pdf = GLMakie.AbstractPlotting.barplot!(
        hslicescene,
        xs,
        ys,
        color = :blue,
        strokecolor = :blue,
        strokewidth = 1,
    )

    @lift(GLMakie.AbstractPlotting.xlims!(hslicescene, extrema($xs)))
    @lift(GLMakie.AbstractPlotting.ylims!(hslicescene, extrema($ys)))
    vlines!(
        hslicescene,
        @lift($sliceclims[1]),
        color = :black,
        linewidth = menuwidth / 100,
    )
    vlines!(
        hslicescene,
        @lift($sliceclims[2]),
        color = :black,
        linewidth = menuwidth / 100,
    )

    interpolationnames = ["contour", "heatmap"]
    interpolationchoices = [true, false]
    interpolationnode = Node(interpolationchoices[1])
    interpolationmenu =
        Menu(scene, options = zip(interpolationnames, interpolationchoices))

    on(interpolationmenu.selection) do s
        interpolationnode[] = s
        # hack
        heatmap!(
            slicescene,
            slicexaxis,
            sliceyaxis,
            sliced_state,
            interpolate = s,
            colormap = slicecolormapnode,
            colorrange = sliceclims,
        )
    end

    directionmenu = Menu(scene, options = zip(directionnames, directionindex))

    on(directionmenu.selection) do s
        directionnode[] = s
    end

    slicemenustring = @lift(
        directionnames[$directionnode] *
        " at index " *
        string(round(Int, 1 + $slice_node * ($nr[$directionnode] - 1)))
    )
    lowerclim_string = @lift(
        "quantile = " *
        @sprintf("%0.2f", $lowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[1])
    )
    upperclim_string = @lift(
        "quantile = " *
        @sprintf("%0.2f", $upperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $clims[2])
    )
    alphastring = @lift("Slice alpha = " * @sprintf("%0.2f", $alphanode))
    layout[2, 1] = vgrid!(
        Label(scene, "State", width = nothing),
        statemenu,
        Label(scene, "Color", width = nothing),
        colormenu,
        Label(scene, "Slice Direction", width = nothing),
        directionmenu,
        Label(scene, alphastring, width = nothing),
        alpha_slider,
        Label(scene, slicemenustring, width = nothing),
        slice_slider,
        Label(scene, lowerclim_string, width = nothing),
        lowerclim_slider,
        Label(scene, upperclim_string, width = nothing),
        upperclim_slider,
    )

    slicelowerclim_string = @lift(
        "quantile = " *
        @sprintf("%0.2f", $slicelowerclim_node) *
        ", value = " *
        @sprintf("%0.1e", $sliceclims[1])
    )
    sliceupperclim_string = @lift(
        "quantile = " *
        @sprintf("%0.2f", $sliceupperclim_node) *
        ", value = " *
        @sprintf("%0.1e", $sliceclims[2])
    )

    layout[2, 7] = vgrid!(
        Label(scene, "Contour Plot Type", width = nothing),
        interpolationmenu,
        Label(scene, slicelowerclim_string, width = nothing),
        slicelowerclim_slider,
        Label(scene, sliceupperclim_string, width = nothing),
        sliceupperclim_slider,
        cbar,
    )

    display(scene)
    return scene
end


text:
export marginalrule

@marginalrule Dirichlet(:out_a) (m_out::Dirichlet, m_a::PointMass) = begin
    return (out = prod(ProdAnalytical(), Dirichlet(mean(m_a)), m_out), a = m_a)
end


text:
module BenchMapCatEduction

using BenchmarkTools
using Transducers

const SUITE = BenchmarkGroup()

SUITE["xf-eduction"] =
    @benchmarkable foldl(+, 1:1000 |> MapCat(x -> 1:x |> Map(x -> x + 1)))
SUITE["xf-xf"] = @benchmarkable foldl(+, 1:1000 |> MapCat(x -> 1:x) |> Map(x -> x + 1))
SUITE["base"] = @benchmarkable sum(y + 1 for x in 1:1000 for y in 1:x)

end  # module
BenchMapCatEduction.SUITE


text:
# Entry for script
function main_FVGQ_1_joint(args=ARGS)
    d = parse_commandline_FVGQ_1(args)
    return estimate_FVGQ_1_joint((; d...)) # to named tuple
end

function estimate_FVGQ_1_joint(d)
    # Or move these into main package when loading?
    Turing.setadbackend(:zygote)
    HMCExamples.set_BLAS_threads()
    use_tensorboard = true # could add toggle later

    # load data relative to the current path
    data_path = joinpath(pkgdir(HMCExamples), d.data_path)
    df = Matrix(DataFrame(CSV.File(data_path)))
    z = [df[i, :] for i in 1:size(df, 1)]

    ϵ0 = Matrix(DataFrame(CSV.File(joinpath(pkgdir(HMCExamples), "data/epsilons_burnin_FVGQ20_1.csv");header=false)))

    # Create the perturbation and the turing models
    m = FirstOrderPerturbationModel(FVGQ20_1)
    #create H prior
    Hx = zeros(6, m.n_x)
    Hy = zeros(6, m.n_y)
    Hy[1, 19] = 1 # Π
    Hy[2, 20] = 1 # R
    Hy[3, 18] = 1 # dw, trend is μz
    Hy[4, 18] = 1 # dy, trend is μz
    Hy[6, 24] = 1 # μ-1
    params = (
        β = Gamma_tr(d.beta_prior[1],d.beta_prior[2]),
        h = Beta_tr(d.h_prior[1], d.h_prior[2]),
        κ = (d.kappa_prior[1], d.kappa_prior[2], d.kappa_prior[3], d.kappa_prior[4]),
        α = (d.alpha_prior[1], d.alpha_prior[2], d.alpha_prior[3], d.alpha_prior[4]),
        γy = (d.gamma_y_prior[1], d.gamma_y_prior[2], d.gamma_y_prior[3],d.gamma_y_prior[4]),
        γΠ = (d.gamma_Pi_prior[1], d.gamma_Pi_prior[2], d.gamma_Pi_prior[3], d.gamma_Pi_prior[4]),
        θp = Beta_tr(d.theta_p_prior[1], d.theta_p_prior[2]),
        χ = Beta_tr(d.chi_prior[1], d.chi_prior[2]),
        γR = Beta_tr(d.gamma_R_prior[1], d.gamma_R_prior[2]),
        Πbar = Gamma_tr(d.Pi_bar_prior[1], d.Pi_bar_prior[2]),
        ρd = Beta_tr(d.rho_d_prior[1], d.rho_d_prior[2]),
        ρφ = Beta_tr(d.rho_psi_prior[1], d.rho_psi_prior[2]),
        ρg = Beta_tr(d.rho_g_prior[1], d.rho_g_prior[2]),
        g_bar = Beta_tr(d.g_bar_prior[1], d.g_bar_prior[2]),
        σ_A = InvGamma_tr(d.sigma_A_prior[1], d.sigma_A_prior[2]),
        σ_d = InvGamma_tr(d.sigma_d_prior[1], d.sigma_d_prior[2]),
        σ_φ = InvGamma_tr(d.sigma_psi_prior[1], d.sigma_psi_prior[2]),
        σ_μ = InvGamma_tr(d.sigma_mu_prior[1], d.sigma_mu_prior[2]),
        σ_m = InvGamma_tr(d.sigma_m_prior[1], d.sigma_m_prior[2]),
        σ_g = InvGamma_tr(d.sigma_g_prior[1], d.sigma_g_prior[2]),
        Λμ = Gamma_tr(d.Lambda_mu_prior[1], d.Lambda_mu_prior[2]),
        ΛA = Gamma_tr(d.Lambda_A_prior[1], d.Lambda_A_prior[2]),
        Hx = Hx,
        Hy = Hy
    )

    turing_model = FVGQ20_joint(
        z, m, d.p_f, params, allocate_cache(m)
    )

    # Sampler
    name = "FQGV-joint-1-s$(d.num_samples)-seed$(d.seed)"
    include_vars = ["β_draw", "h", "κ", "α", "θp", "χ", "γR", "γy", "γΠ", "Πbar_draw", "ρd", "ρφ", "ρg", "g_bar", "σ_A", "σ_d", "σ_φ", "σ_μ", "σ_m", "σ_g", "Λμ", "ΛA"]  # variables to log
    callback = TensorBoardCallback(d.results_path; name, include=include_vars)
    num_adapts = convert(Int64, floor(d.num_samples * d.adapts_burnin_prop))

    Random.seed!(d.seed)
    @info "Generating $(d.num_samples) samples with $(num_adapts) adapts"
    alg = NUTS(num_adapts, d.target_acceptance_rate)

    chain = sample(
        turing_model,
        NUTS(num_adapts, d.target_acceptance_rate; d.max_depth),
        d.num_samples;
        init_params=[d.p,ϵ0],
        progress=true,
        save_state=true,
        callback,
    )

    # Calculate and save results into the logdir
    calculate_experiment_results(chain, callback.logger, include_vars, d.full_results)

    # Store parameters in log directory
    parameter_save_path = joinpath(callback.logger.logdir, "parameters.json")

    @info "Storing Parameters at $(parameter_save_path) "
    open(parameter_save_path, "w") do f
        write(f, JSON.json(d))
    end
end

function parse_commandline_FVGQ_1(args)
    s = ArgParseSettings(; fromfile_prefix_chars=['@'])

    # See the appropriate _defaults.txt file for the default vvalues.
    @add_arg_table! s begin
        "--data_path"
        help = "relative path to data from the root of the package"
        arg_type = String
        "--p"
        help = "Initialization of parameters"
        arg_type = Vector{Float64}
        "--p_f"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}

        "--alpha_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}
        "--kappa_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}
        "--gamma_y_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}
        "--gamma_Pi_prior"
        help = "Value of fixed parameters"
        arg_type = Vector{Float64}

        "--beta_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--h_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--theta_p_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--chi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--gamma_R_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--Pi_bar_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_d_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_psi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--rho_g_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--g_bar_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        "--sigma_A_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_d_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_psi_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_mu_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_m_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--sigma_g_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        "--Lambda_mu_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--Lambda_A_prior"
        help = "Parameters for the prior"
        arg_type = Vector{Float64}

        help = "Parameters for the prior"
        arg_type = Vector{Float64}
        "--num_samples"
        help = "samples to draw in chain"
        arg_type = Int64
        "--adapts_burnin_prop"
        help = "Proportion of Adaptations burned in"
        arg_type = Float64
        "--target_acceptance_rate"
        help = "Target acceptance rate for dual averaging for NUTS"
        arg_type = Float64
        "--max_depth"
        help = "Maximum depth for NUTS"
        arg_type = Int64
        "--seed"
        help = "Random number seed"
        arg_type = Int64
        "--results_path"
        arg_type = String
        help = "Location to store results and logs"

        "--full_results"
        arg_type = Bool
        help = "Save the complete set of figures and results for the chain"
    end

    args_with_default = vcat("@$(pkgdir(HMCExamples))/src/FVGQ_1_joint_defaults.txt", args)
    return parse_args(args_with_default, s; as_symbols=true)

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Coloring

import DataStructures

include("topological_sort.jl")

# indexed sparse set of integers
mutable struct IndexedSet
    nzidx::Vector{Int}
    empty::BitVector
    nnz::Int
end

IndexedSet(n::Integer) = IndexedSet(zeros(Int, n), trues(n), 0)

function Base.push!(v::IndexedSet, i::Integer)
    if v.empty[i]  # new index
        v.nzidx[v.nnz+=1] = i
        v.empty[i] = false
    end
    return
end

function Base.empty!(v::IndexedSet)
    for i in 1:v.nnz
        v.empty[v.nzidx[i]] = true
    end
    v.nnz = 0
    return v
end

Base.length(v::IndexedSet) = length(v.nzidx)

function Base.resize!(v::IndexedSet, n::Integer)
    if n > length(v)
        @assert v.nnz == 0 # only resize empty vector
        resize!(v.nzidx, n)
        resize!(v.empty, n)
        fill!(v.empty, true)
    end
    return
end

Base.collect(v::IndexedSet) = v.nzidx[1:v.nnz]

function Base.union!(v::IndexedSet, s)
    for x in s
        push!(v, x)
    end
    return
end

# compact storage for an undirected graph
# neighbors of vertex i start at adjlist[offsets[i]]
struct UndirectedGraph
    adjlist::Vector{Int}
    edgeindex::Vector{Int} # corresponding edge number, indexed as adjlist
    offsets::Vector{Int}
    edges::Vector{Tuple{Int,Int}}
end

_num_vertices(g::UndirectedGraph) = length(g.offsets) - 1

_num_edges(g::UndirectedGraph) = length(g.edges)

_num_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i+1] - g.offsets[i]

_start_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i]

function UndirectedGraph(I, J, nel)
    adjcount = zeros(Int, nel)
    n_edges = 0
    for k in 1:length(I)
        i = I[k]
        j = J[k]
        if i == j
            continue
        end
        n_edges += 1
        adjcount[i] += 1
        adjcount[j] += 1
    end
    offsets = Vector{Int}(undef, nel + 1)
    offsets[1] = 1
    for k in 1:nel
        offsets[k+1] = offsets[k] + adjcount[k]
    end
    fill!(adjcount, 0)
    edges = Vector{Tuple{Int,Int}}(undef, n_edges)
    adjlist = Vector{Int}(undef, offsets[nel+1] - 1)
    edgeindex = Vector{Int}(undef, length(adjlist))
    edge_count = 0
    for k in 1:length(I)
        i = I[k]
        j = J[k]
        if i == j
            continue
        end
        edge_count += 1
        adjlist[offsets[i]+adjcount[i]] = j
        edgeindex[offsets[i]+adjcount[i]] = edge_count
        adjcount[i] += 1
        adjlist[offsets[j]+adjcount[j]] = i
        edgeindex[offsets[j]+adjcount[j]] = edge_count
        adjcount[j] += 1
        edges[edge_count] = (i, j)
    end
    @assert edge_count == n_edges
    return UndirectedGraph(adjlist, edgeindex, offsets, edges)
end

struct _Edge
    index::Int
    source::Int
    target::Int
end

function _prevent_cycle(
    v,
    w,
    x,
    e_idx1,
    e_idx2,
    S,
    firstVisitToTree,
    forbiddenColors,
    color,
)
    er = DataStructures.find_root!(S, e_idx2)
    @inbounds first = firstVisitToTree[er]
    p = first.source # but this depends on the order?
    q = first.target
    @inbounds if p != v
        firstVisitToTree[er] = _Edge(e_idx1, v, w)
    elseif q != w
        forbiddenColors[color[x]] = v
    end
    return
end

function _grow_star(v, w, e_idx, firstNeighbor, color, S)
    @inbounds e = firstNeighbor[color[w]]
    p = e.source
    @inbounds if p != v
        firstNeighbor[color[w]] = _Edge(e_idx, v, w)
    else
        union!(S, e_idx, e.index)
    end
    return
end

function _merge_trees(eg, eg1, S)
    e1 = DataStructures.find_root!(S, eg)
    e2 = DataStructures.find_root!(S, eg1)
    if e1 != e2
        union!(S, eg, eg1)
    end
    return
end

"""
    acyclic_coloring(g::UndirectedGraph)

Implement the acyclic coloring algorithm of Gebremdehin, Tarafdar, Manne, and
Pothen, "New Acyclic and Star Coloring Algorithms with Application to Computing
Hessians." SIAM J. Sci. Comput. 2007.

Returns `Tuple{Vector{Int},Int}` giving the color index of each node in `g`, as
well as the total number of colors used.
"""
function acyclic_coloring(g::UndirectedGraph)
    if _num_edges(g) == 0
        return fill(1, _num_vertices(g)), 1
    end
    num_colors = 0
    forbiddenColors = Int[]
    firstNeighbor = _Edge[]
    firstVisitToTree = fill(_Edge(0, 0, 0), _num_edges(g))
    color = fill(0, _num_vertices(g))
    # disjoint set forest of edges in the graph
    S = DataStructures.IntDisjointSets(_num_edges(g))
    @inbounds for v in 1:_num_vertices(g)
        n_neighbor = _num_neighbors(v, g)
        start_neighbor = _start_neighbors(v, g)
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            if color[w] == 0
                continue
            end
            forbiddenColors[color[w]] = v
        end
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            if color[w] == 0
                continue
            end
            n_neighbor_w = _num_neighbors(w, g)
            start_neighbor_w = _start_neighbors(w, g)
            for k2 in 0:(n_neighbor_w-1)
                x = g.adjlist[start_neighbor_w+k2]
                e2_idx = g.edgeindex[start_neighbor_w+k2]
                if color[x] == 0
                    continue
                end
                if forbiddenColors[color[x]] != v
                    _prevent_cycle(
                        v,
                        w,
                        x,
                        e_idx,
                        e2_idx,
                        S,
                        firstVisitToTree,
                        forbiddenColors,
                        color,
                    )
                end
            end
        end
        # find feasible color
        found = false
        for k in 1:num_colors
            if forbiddenColors[k] != v
                color[v] = k
                found = true
                break
            end
        end
        if !found
            num_colors += 1
            push!(forbiddenColors, 0)
            push!(firstNeighbor, _Edge(0, 0, 0))
            color[v] = num_colors
        end
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            if color[w] == 0
                continue
            end
            _grow_star(v, w, e_idx, firstNeighbor, color, S)
        end
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            if color[w] == 0
                continue
            end
            n_neighbor_w = _num_neighbors(w, g)
            start_neighbor_w = _start_neighbors(w, g)
            for k2 in 0:(n_neighbor_w-1)
                x = g.adjlist[start_neighbor_w+k2]
                e2_idx = g.edgeindex[start_neighbor_w+k2]
                if color[x] == 0 || x == v
                    continue
                end
                if color[x] == color[v]
                    _merge_trees(e_idx, e2_idx, S)
                end
            end
        end
    end
    return color, num_colors
end

struct RecoveryInfo
    vertexmap::Vector{Vector{Int}}
    postorder::Vector{Vector{Int}}
    parents::Vector{Vector{Int}}
    color::Vector{Int}
    num_colors::Int
    nnz::Int # number of off-diagonal nonzeros
    local_indices::Vector{Int} # map back to global indices
end

function RecoveryInfo()
    return RecoveryInfo(
        Vector{Vector{Int}}(undef, 0),
        Vector{Vector{Int}}(undef, 0),
        Vector{Vector{Int}}(undef, 0),
        Vector{Int}(undef, 0),
        0,
        0,
        Vector{Int}(undef, 0),
    )
end

function recovery_preprocess(
    g::UndirectedGraph,
    color,
    num_colors,
    local_indices,
)
    # represent two-color subgraph as:
    # list of vertices (with map to global indices)
    # adjacency list in a single vector (with list of offsets)
    # linear index of pair of colors
    twocolorindex = zeros(Int32, num_colors, num_colors)
    seen_twocolors = 0
    # count of edges in each subgraph
    edge_count = Int[]
    for k in 1:length(g.edges)
        u, v = g.edges[k]
        i = min(color[u], color[v])
        j = max(color[u], color[v])
        if twocolorindex[i, j] == 0
            seen_twocolors += 1
            twocolorindex[i, j] = seen_twocolors
            push!(edge_count, 0)
        end
        idx = twocolorindex[i, j]
        edge_count[idx] += 1
    end
    # edges sorted by twocolor subgraph
    sorted_edges = Vector{Vector{Tuple{Int,Int}}}(undef, seen_twocolors)
    for idx in 1:seen_twocolors
        sorted_edges[idx] = Tuple{Int,Int}[]
        sizehint!(sorted_edges[idx], edge_count[idx])
    end
    for i in 1:length(g.edges)
        u, v = g.edges[i]
        i = min(color[u], color[v])
        j = max(color[u], color[v])
        idx = twocolorindex[i, j]
        push!(sorted_edges[idx], (u, v))
    end
    # list of unique vertices in each twocolor subgraph
    vertexmap = Vector{Vector{Int}}(undef, seen_twocolors)
    postorder = Vector{Vector{Int}}(undef, seen_twocolors)
    parents = Vector{Vector{Int}}(undef, seen_twocolors)
    # temporary lookup map from global index to subgraph index
    revmap = zeros(Int, _num_vertices(g))
    adjcount = zeros(Int, _num_vertices(g))
    cmap = zeros(Int, 0) # shared storage for DFS
    for idx in 1:seen_twocolors
        my_edges = sorted_edges[idx]
        vlist = Int[]
        # build up the vertex list and adjacency count
        for k in 1:length(my_edges)
            u, v = my_edges[k]
            # seen these vertices yet?
            if revmap[u] == 0
                push!(vlist, u)
                revmap[u] = length(vlist)
                adjcount[u] = 0
            end
            if revmap[v] == 0
                push!(vlist, v)
                revmap[v] = length(vlist)
                adjcount[v] = 0
            end
            adjcount[u] += 1
            adjcount[v] += 1
        end
        # set up offsets for adjlist
        offset = Vector{Int}(undef, length(vlist) + 1)
        offset[1] = 1
        for k in 1:length(vlist)
            offset[k+1] = offset[k] + adjcount[vlist[k]]
            adjcount[vlist[k]] = 0
        end
        # adjlist for node u in twocolor idx starts at
        # vec[offset[u]]
        # u has global index vlist[u]
        vec = Vector{Int}(undef, offset[length(vlist)+1] - 1)
        # now fill in
        for k in 1:length(my_edges)
            u, v = my_edges[k]
            u_rev = revmap[u] # indices in the subgraph
            v_rev = revmap[v]
            vec[offset[u_rev]+adjcount[u]] = v_rev
            vec[offset[v_rev]+adjcount[v]] = u_rev
            adjcount[u] += 1
            adjcount[v] += 1
        end
        resize!(cmap, length(vlist))
        order, parent =
            reverse_topological_sort_by_dfs(vec, offset, length(vlist), cmap)
        for k in 1:length(vlist)
            revmap[vlist[k]] = 0  # clear for reuse
        end
        postorder[idx] = order
        parents[idx] = parent
        vertexmap[idx] = vlist
    end
    return RecoveryInfo(
        vertexmap,
        postorder,
        parents,
        color,
        num_colors,
        _num_edges(g),
        local_indices,
    )
end

_normalize(i, j) = (j > i) ? (j, i) : (i, j)

function _indirect_recover_structure(rinfo::RecoveryInfo)
    N = length(rinfo.color)
    I = zeros(Int, rinfo.nnz + N)
    J = zeros(Int, rinfo.nnz + N)
    k = 0
    for i in 1:N
        k += 1
        I[k] = J[k] = i
    end
    for t in 1:length(rinfo.postorder)
        vmap = rinfo.vertexmap[t]
        order = rinfo.postorder[t]
        parent = rinfo.parents[t]
        for z in 1:length(order)
            v = order[z]
            p = parent[v]
            if p == 0
                continue
            end
            k += 1
            I[k], J[k] = _normalize(vmap[v], vmap[p])
        end
    end
    @assert k == rinfo.nnz + N
    return I, J
end

"""
    hessian_color_preprocess(
        edgelist,
        num_total_var,
        seen_idx = IndexedSet(0),
    )

`edgelist` contains the nonzeros in the Hessian, *including* nonzeros on the
diagonal.
"""
function hessian_color_preprocess(
    edgelist,
    num_total_var,
    seen_idx = IndexedSet(0),
)
    resize!(seen_idx, num_total_var)
    I, J = Int[], Int[]
    for (i, j) in edgelist
        push!(seen_idx, i)
        push!(seen_idx, j)
        push!(I, i)
        push!(J, j)
    end
    local_indices = sort!(seen_idx.nzidx[1:seen_idx.nnz])
    empty!(seen_idx)
    global_to_local_idx = seen_idx.nzidx # steal for storage
    for k in 1:length(local_indices)
        global_to_local_idx[local_indices[k]] = k
    end
    # only do the coloring on the local indices
    for k in 1:length(I)
        I[k] = global_to_local_idx[I[k]]
        J[k] = global_to_local_idx[J[k]]
    end
    g = UndirectedGraph(I, J, length(local_indices))
    color, num_colors = acyclic_coloring(g)
    @assert length(color) == _num_vertices(g)
    rinfo = recovery_preprocess(g, color, num_colors, local_indices)
    I, J = _indirect_recover_structure(rinfo)
    # convert back to global indices
    for k in 1:length(I)
        I[k] = local_indices[I[k]]
        J[k] = local_indices[J[k]]
    end
    return I, J, rinfo
end

"""
    seed_matrix(rinfo::RecoveryInfo)

Allocate a seed matrix for the Coloring.
"""
function seed_matrix(rinfo::RecoveryInfo)
    return Matrix{Float64}(undef, length(rinfo.local_indices), rinfo.num_colors)
end

"""
    prepare_seed_matrix!(R, rinfo::RecoveryInfo)
"""
function prepare_seed_matrix!(R, rinfo::RecoveryInfo)
    N = length(rinfo.color)
    @assert N == size(R, 1) == length(rinfo.local_indices)
    @assert size(R, 2) == rinfo.num_colors
    fill!(R, 0.0)
    for i in 1:N
        R[i, rinfo.color[i]] = 1
    end
    return
end

"""
    recover_from_matmat!(
        V::AbstractVector{T},
        R::AbstractMatrix{T},
        rinfo::RecoveryInfo,
        stored_values::AbstractVector{T},
    ) where {T}

recover the hessian values from the hessian-matrix solution
stored_values is a temporary vector of length >= length(rinfo.local_indices)
"""
function recover_from_matmat!(
    V::AbstractVector{T},
    R::AbstractMatrix{T},
    rinfo::RecoveryInfo,
    stored_values::AbstractVector{T},
) where {T}
    N = length(rinfo.color) # number of local variables
    @assert length(V) == rinfo.nnz + N
    @assert length(stored_values) >= length(rinfo.local_indices)
    k = 0
    for i in 1:N
        k += 1
        V[k] = R[i, rinfo.color[i]]
    end
    for t in 1:length(rinfo.vertexmap)
        vmap = rinfo.vertexmap[t]
        order = rinfo.postorder[t]
        parent = rinfo.parents[t]
        for z in 1:length(order)
            @inbounds stored_values[z] = 0.0
        end
        @inbounds for z in 1:length(order)
            v = order[z]
            p = parent[v]
            if p == 0
                continue
            end
            i, j = vmap[v], vmap[p]
            value = R[i, rinfo.color[j]] - stored_values[v]
            stored_values[p] += value
            k += 1
            V[k] = value
        end
    end
    @assert k == rinfo.nnz + N
    return
end

end  # module

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
function max_heapify!(a::Array, s::Int, i::Int)
	l = i * 2
	r = i * 2 + 1
	m = 0
	if l <= s && a[l] > a[i]
		m = l
	else
		m = i
	end

	if r <= s && a[r] > a[m]
		m = r
	end

	if m != i
		a[i], a[m] = a[m], a[i]
		max_heapify!(a, s, m)
	end
end

function make_heap!(a::Array)
	s = length(a)
	for i in fld(length(a), 2) : -1 : 1
		max_heapify!(a, s, i)
	end
end

function heapsort!(a::Array)
	make_heap!(a)
	s = length(a)
	for i in length(a) : -1 : 2
		a[1], a[s] = a[s], a[1]
		s -= 1
		max_heapify!(a, s, 1)
	end
end


text:
# From Julia 1.0's online docs. File countheads.jl available to all machines:

function count_heads(n)
    c::Int = 0
    for i = 1:n
        c += rand(Bool)
    end
    c
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
expr:
- Julia version: 1.3.1
- Author: Berk
- Date: 2020-11-04
Note: This is a dump file when BBFs used to take raw expressions instead of Expr(Function)'s
=#

"""
Generates a map for the auxiliary variables for Expr -> Function conversion.
"""
function generate_auxmap(ex::Expr, vars::Array{JuMP.VariableRef}, model::JuMP.Model)
    outers = get_outers(ex)
    outer_var_mapping = outers_to_vars(outers, model)
    dictmap = Dict()
    count = 1;
    for (key, value) in outer_var_mapping
        if value isa Array
            dictmap[key] = [:(aux[$i]) for i = count:count+length(value)]
            count += length(value)
        else
            count += 1
            dictmap[key] = :(aux[$count])
        end
    end
    new_expr = copy(ex);
    new_expr = substitute_expr(new_expr, dictmap)
    return
end

import MLStyle
import JuliaVariables

"""
    get_locals(ex::Expr)

Get variables in Expr that are locally defined.

get_locals borrowed from here:
https://discourse.julialang.org/t/determine-the-scope-of-symbols-in-an-expr/47993/2
"""
function get_locals(ex::Expr)
    vars = (JuliaVariables.solve_from_local ∘ JuliaVariables.simplify_ex)(ex).args[1].bounds
    map(x -> x.name, vars)
end

_get_outers(_) = Symbol[]
_get_outers(x::JuliaVariables.Var) = x.is_global ? [x.name] : Symbol[]
function _get_outers(ex::Expr)
    MLStyle.@match ex begin
        Expr(:(=), _, rhs) => _get_outers(rhs)
        Expr(:tuple, _..., Expr(:(=), _, rhs)) => _get_outers(rhs)
        Expr(_, args...) => mapreduce(_get_outers, vcat, args)
    end
end

"""
    get_outers(ex::Expr)

Gets the Symbols in the expression that are not local.

get_outers borrowed from here:
https://discourse.julialang.org/t/determine-the-scope-of-symbols-in-an-expr/47993/2
"""
get_outers(ex) = (unique! ∘ _get_outers ∘ JuliaVariables.solve_from_local ∘ JuliaVariables.simplify_ex)(ex)

# Dummy code generating functions from Exprs. Could come in useful.

constr_fn = let lhs = lhsexpr, rhs = rhsexpr, op = op, consts = consts
    if op in [:(=), :>]
        function (x)
            for (key, value) in union(x, consts, sets)
                eval(Meta.parse("$key = $value"))
            end
            lhs_evaled = eval(lhs)
            rhs_evaled = eval(rhs)
            return lhs_evaled-rhs_evaled
        end
    elseif op == :<
        function (x)
            for (key, value) in union(x, consts, sets)
                eval(Meta.parse("$key = $value"))
            end
            lhs_evaled = eval(lhs)
            rhs_evaled = eval(rhs)
            return rhs_evaled-lhs_evaled
        end
    else
        error("unexpected operator ", op)
    end

"""
    substitute_expr(expr::Expr, data::Dict)

Substitutes for the outer variables (symbols from the expr) using data.
"""
function substitute_expr(expr::Expr, data::Dict)
    new_expr = copy(expr)
    for (key, value) in data
        new_expr = substitute(new_expr, :($key) => value)
    end
    return eval(new_expr)
end

"""
    substitute_args is a helper function, by Oscar Dowson.
    See here for more details:
    https://discourse.julialang.org/t/procedural-nonlinear-constaint-generation-string-to-nlconstraint/34799/4
"""
substitute_args(ex, vars) = ex
substitute_args(ex::Symbol, vars) = get(vars, ex, ex)
function substitute_args(ex::Expr, vars)
    for (i, arg) in enumerate(ex.args)
        ex.args[i] = substitute_args(arg, vars)
    end
    return ex
end


text:
import SpecialFunctions: erfc
import IntervalArithmetic: Interval, mid

# Monkey-patch a few functions for Intervals
# ... this is far from proper and a bit specific for our use case here
# (that's why it's not contributed upstream).
# should be done e.g. by changing  the rounding mode ...
erfc(i::Interval) = Interval(prevfloat(erfc(i.lo)), nextfloat(erfc(i.hi)))

function compute_Glims_fast(lattice::AbstractMatrix{<:Interval}, args...; kwargs...)
    # This is done to avoid a call like ceil(Int, ::Interval)
    # in the above implementation of compute_fft_size,
    # where it is in general cases not clear, what to do.
    # In this case we just want a reasonable number for Gmax,
    # so replacing the intervals in the lattice with
    # their midpoints should be good.
    compute_Glims_fast(mid.(lattice), args...; kwargs...)
end
function compute_Glims_precise(::AbstractMatrix{<:Interval}, args...; kwargs...)
    error("fft_size_algorithm :precise not supported with intervals")
end

function _is_well_conditioned(A::AbstractArray{<:Interval}; kwargs...)
    # This check is used during the lattice setup, where it frequently fails with intervals
    # (because doing an SVD with intervals leads to a large overestimation of the rounding error)
    _is_well_conditioned(mid.(A); kwargs...)
end

function symmetry_operations(lattice::AbstractMatrix{<:Interval}, atoms, magnetic_moments=[];
                             tol_symmetry=max(1e-5, maximum(radius, lattice)))
    @assert tol_symmetry < 1e-2
    symmetry_operations(mid.(lattice), atoms, magnetic_moments; tol_symmetry)
end

function local_potential_fourier(el::ElementCohenBergstresser, q::T) where {T <: Interval}
    lor = round(q.lo, digits=5)
    hir = round(q.hi, digits=5)
    @assert iszero(round(lor - hir, digits=3))
    T(local_potential_fourier(el, mid(q)))
end


text:
function runge_kutta(y0::Float64, t0::Float64, h::Float64, steps::Int, func::AbstractString)
    result = Float64[]
    t_array = Float64[]
    expression = parse(func)
    global y = y0
    global t = t0
    push!(result, y0)
    push!(t_array, t0)

    for i = 2:(steps+1)
        tn = t
        kn1 = eval(expression)
        t = tn + h/2
        y = result[i-1] + kn1*h/2
        kn2 = eval(expression)
        t = tn + h/2
        y = result[i-1] + kn2*h/2
        kn3 = eval(expression)
        t = tn + h
        y = result[i-1] + kn3*h
        kn4 = eval(expression)
        push!(result,result[i-1] + ((kn1+2*kn2+2kn3+kn4)/6)*h)

        t = tn+h
        push!(t_array, t)
        y = result[i];
    end
    return result, t_array
end

result, t_array = runge_kutta(0.0, 0.0, 0.1, 20, "1+2t")
for i in result
	println(i)
end

max_stars_repo_name: CarbonDDR/al-go-rithms


text:
typealias DefaultRefType UInt32

## Pools

# V is always set to CategoricalValue{T}
# This workaround is needed since this type not defined yet
# See JuliaLang/julia#269
type CategoricalPool{T, R <: Integer, V}
    index::Vector{T}
    invindex::Dict{T, R}
    order::Vector{R}
    levels::Vector{T}
    valindex::Vector{V}
    ordered::Bool

    function CategoricalPool{T, R}(index::Vector{T},
                                   invindex::Dict{T, R},
                                   order::Vector{R},
                                   ordered::Bool)
        pool = new(index, invindex, order, index[order], V[], ordered)
        buildvalues!(pool)
        return pool
    end
end

## Values

immutable CategoricalValue{T, R <: Integer}
    level::R
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end

## Arrays

abstract AbstractCategoricalArray{T, N, R} <: AbstractArray{CategoricalValue{T, R}, N}

type CategoricalArray{T, N, R <: Integer} <: AbstractCategoricalArray{T, N, R}
    refs::Array{R, N}
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end
typealias CategoricalVector{T, R} CategoricalArray{T, 1, R}
typealias CategoricalMatrix{T, R} CategoricalArray{T, 2, R}

## Nullable Arrays

abstract AbstractNullableCategoricalArray{T, N, R} <: AbstractArray{Nullable{CategoricalValue{T, R}}, N}

type NullableCategoricalArray{T, N, R <: Integer} <: AbstractNullableCategoricalArray{T, N, R}
    refs::Array{R, N}
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end
typealias NullableCategoricalVector{T, R} NullableCategoricalArray{T, 1, R}
typealias NullableCategoricalMatrix{T, R} NullableCategoricalArray{T, 2, R}

## Type Aliases

typealias CatArray Union{CategoricalArray, NullableCategoricalArray}