text:
@inline g(contact::Contact, x) = contact.Nx[SA[1; 2; 3]]' * (x - contact.offset[SA[1; 2; 3]])
@inline g(contact::Contact, state::State, Δt) = g(contact, getx3(state, Δt))

@inline ∂g∂pos(contact::Contact) = contact.Nx
@inline ∂g∂pos(contact::Contact, state::State) = ∂g∂pos(contact)

@inline ∂g∂vel(contact::Contact, Δt) = contact.Nx * Δt
@inline ∂g∂vel(contact::Contact, state::State, Δt) = ∂g∂vel(contact, Δt)

@inline schurf(contact::Contact, φ, γ, s, μ) = ∂g∂pos(contact)' * (γ / s * φ - μ / s)
@inline function schurf(contact::Contact, state::State, γ, s, μ, Δt)
    φ = g(contact, getx3(state, Δt))
    return schurf(contact, φ, γ, s, μ)
end

@inline function schurD(contact::Contact, γ, s, Δt)
    Nx = ∂g∂pos(contact)
    Nv = ∂g∂vel(contact, Δt)

    return Nx' * γ / s * Nv
end
@inline schurD(contact::Contact, state::State, γ, s, Δt) = schurD(contact, γ, s, Δt)

@inline function calcFrictionForce!(mechanism::Mechanism{T}, friction::Contact, body::Body, γ) where T
    cf = friction.cf
    D = friction.D
    state = body.state

    d = state.d
    v = state.vsol[2]
    ω = state.ωsol[2]
    state.vsol[2] = szeros(T, 3)
    state.ωsol[2] = szeros(T, 3)
    dyn = dynamics(mechanism, body)
    state.vsol[2] = v
    ω = state.ωsol[2] = ω
    state.d = d

    b0 = D*dyn + friction.b # remove old friction force

    if norm(b0) > cf*γ
        friction.b = b0/norm(b0)*cf*γ
    else
        friction.b = b0
    end
    return
end


text:
using ITensors

let
  Ny = 6
  Nx = 12

  N = Nx*Ny

  sites = siteinds("S=1/2", N;
                   conserve_qns = true)

  lattice = square_lattice(Nx, Ny; yperiodic = false)

  ampo = AutoMPO()
  for b in lattice
    ampo .+= 0.5, "S+", b.s1, "S-", b.s2
    ampo .+= 0.5, "S-", b.s1, "S+", b.s2
    ampo .+=      "Sz", b.s1, "Sz", b.s2
  end
  H = MPO(ampo,sites)

  state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  # Initialize wavefunction to a random MPS
  # of bond-dimension 10 with same quantum
  # numbers as `state`
  psi0 = randomMPS(sites,state,20)

  sweeps = Sweeps(10)
  setmaxdim!(sweeps,20,60,100,100,200,400,800)
  setcutoff!(sweeps,1E-8)
  @show sweeps

  energy,psi = dmrg(H,psi0,sweeps)

  return
end

ext: jl


text:
# include("approximators/approximators.jl")
# include("CommandGenerator.jl")
include("costs.jl")
include("partitionTrainTest.jl")
include("ned2enu.jl")


text:
abstract type AbstractSimulation end
struct DefaultSimulation end

simulate(model::Model, args...; kwargs...) = simulate(DefaultSimulation(), model, args...; kwargs...)
function simulate(::DefaultSimulation, model::Model, param, u0)
end


text:
export dump_gate, yaotofile, yaotoscript
export yaofromfile, yaofromstring, @yao_str, @yaoscript
using MLStyle: @match

include("dump.jl")
include("load.jl")
include("address_manipulate.jl")
include("optimise.jl")

print_blocktree() = print_subtypetree(AbstractBlock)


text:
module Calc

include("Add.jl")
include("Multiply.jl")
include("Subtract.jl")
include("Divide.jl")

export add, subtract, divide, multiply

end


text:
# initialise evotree
function init_evotree(params::EvoTypes{T,U,S},
    X::AbstractMatrix, Y::AbstractVector) where {T,U,S}

    K = 1
    levels = nothing
    X = convert(Matrix{T}, X)

    if typeof(params.loss) == Poisson
        Y = T.(Y)
        μ = fill(log(mean(Y)), 1)
    elseif typeof(params.loss) == Softmax
        if eltype(Y) <: CategoricalValue
            levels = CategoricalArrays.levels(Y)
            K = length(levels)
            μ = zeros(T, K)
            Y = UInt32.(CategoricalArrays.levelcode.(Y))
        else
            levels = sort(unique(Y))
            yc = CategoricalVector(Y, levels=levels)
            K = length(levels)
            μ = zeros(T, K)
            Y = UInt32.(CategoricalArrays.levelcode.(yc))
        end
    elseif typeof(params.loss) == Gaussian
        K = 2
        Y = T.(Y)
        μ = [mean(Y), log(std(Y))]
    else
        Y = T.(Y)
        μ = [mean(Y)]
    end

    # initialize preds
    X_size = size(X)
    pred = zeros(T, K, X_size[1])
    @inbounds for i in 1:X_size[1]
        pred[:,i] .= μ
    end

    bias = Tree(μ)
    evotree = GBTree([bias], params, Metric(), K, levels)

    # initialize gradients and weights
    δ𝑤 = ones(T, 2 * K + 1, X_size[1])

    # binarize data into quantiles
    edges = get_edges(X, params.nbins)
    X_bin = binarize(X, edges)

    𝑖_ = UInt32.(collect(1:X_size[1]))
    𝑗_ = UInt32.(collect(1:X_size[2]))
    𝑗 = zeros(eltype(𝑗_), ceil(Int, params.colsample * X_size[2]))

    # initializde histograms
    nodes = [TrainNode(X_size[2], params.nbins, K, T) for n in 1:2^params.max_depth - 1]
    nodes[1].𝑖 = zeros(eltype(𝑖_), ceil(Int, params.rowsample * X_size[1]))
    out = zeros(UInt32, length(nodes[1].𝑖))
    left = zeros(UInt32, length(nodes[1].𝑖))
    right = zeros(UInt32, length(nodes[1].𝑖))

    cache = (params = deepcopy(params),
        X = X, Y = Y, K = K,
        nodes = nodes,
        pred = pred,
        𝑖_ = 𝑖_, 𝑗_ = 𝑗_, 𝑗 = 𝑗,
        out = out, left = left, right = right,
        δ𝑤 = δ𝑤,
        edges = edges,
        X_bin = X_bin)

    cache.params.nrounds = 0

    return evotree, cache
end

function grow_evotree!(evotree::GBTree{T}, cache) where {T,S}

    # initialize from cache
    params = evotree.params
    δnrounds = params.nrounds - cache.params.nrounds

    # loop over nrounds
    for i in 1:δnrounds
        # select random rows and cols
        sample!(params.rng, cache.𝑖_, cache.nodes[1].𝑖, replace=false, ordered=true)
        sample!(params.rng, cache.𝑗_, cache.𝑗, replace=false, ordered=true)

        # build a new tree
        update_grads!(params.loss, cache.δ𝑤, cache.pred, cache.Y, params.α)
        # assign a root and grow tree
        tree = Tree(params.max_depth, evotree.K, zero(T))
        grow_tree!(tree, cache.nodes, params, cache.δ𝑤, cache.edges, cache.𝑗, cache.out, cache.left, cache.right, cache.X_bin, cache.K)
        push!(evotree.trees, tree)
        predict!(params.loss, cache.pred, tree, cache.X, cache.K)

    end # end of nrounds
    cache.params.nrounds = params.nrounds
    return nothing
end

# grow a single tree
function grow_tree!(
    tree::Tree{T},
    nodes::Vector{TrainNode{T}},
    params::EvoTypes{T,U,S},
    δ𝑤::Matrix{T},
    edges,
    𝑗, out, left, right,
    X_bin::AbstractMatrix, K) where {T,U,S}

    # reset nodes
    @threads for n in eachindex(nodes)
        [nodes[n].h[j] .= 0 for j in 𝑗]
        nodes[n].∑ .= 0
        nodes[n].gain = 0
        fill!(nodes[n].gains, -Inf)
    end

    # reset
    n_next = [1]
    n_current = copy(n_next)
    depth = 1

    # initialize summary stats
    nodes[1].∑ .= vec(sum(δ𝑤[:, nodes[1].𝑖], dims=2))
    nodes[1].gain = get_gain(params.loss, nodes[1].∑, params.λ, K)
    # grow while there are remaining active nodes
    while length(n_current) > 0 && depth <= params.max_depth
        offset = 0 # identifies breakpoint for each node set within a depth

        if depth < params.max_depth
            for n_id ∈ 1:length(n_current)
                n = n_current[n_id]
                if n_id % 2 == 0
                    if n % 2 == 0
                        nodes[n].h .= nodes[n >> 1].h .- nodes[n + 1].h
                    else
                        nodes[n].h .= nodes[n >> 1].h .- nodes[n - 1].h
                    end
                else
                    update_hist!(params.loss, nodes[n].h, δ𝑤, X_bin, nodes[n].𝑖, 𝑗, K)
                end
            end
        end

        for n ∈ sort(n_current)
            if depth == params.max_depth || nodes[n].∑[end] <= params.min_weight
                pred_leaf_cpu!(params.loss, tree.pred, n, nodes[n].∑, params, K, δ𝑤, nodes[n].𝑖)
            else
                # histogram subtraction
                update_gains!(params.loss, nodes[n], 𝑗, params, K)
                best = findmax(nodes[n].gains)
                if best[2][1] != params.nbins && best[1] > nodes[n].gain + params.γ
                    tree.gain[n] = best[1] - nodes[n].gain
                    tree.cond_bin[n] = best[2][1]
                    tree.feat[n] = best[2][2]
                    tree.cond_float[n] = edges[tree.feat[n]][tree.cond_bin[n]]
                end
                tree.split[n] = tree.cond_bin[n] != 0
                if !tree.split[n]
                    pred_leaf_cpu!(params.loss, tree.pred, n, nodes[n].∑, params, K, δ𝑤, nodes[n].𝑖)
                    popfirst!(n_next)
                else
                    # println("typeof(nodes[n].𝑖): ", typeof(nodes[n].𝑖))
                    _left, _right = split_set_threads!(out, left, right, nodes[n].𝑖, X_bin, tree.feat[n], tree.cond_bin[n], offset)
                    nodes[n << 1].𝑖, nodes[n << 1 + 1].𝑖 = _left, _right
                    offset += length(nodes[n].𝑖)
                    update_childs_∑!(params.loss, nodes, n, best[2][1], best[2][2], K)
                    nodes[n << 1].gain = get_gain(params.loss, nodes[n << 1].∑, params.λ, K)
                    nodes[n << 1 + 1].gain = get_gain(params.loss, nodes[n << 1 + 1].∑, params.λ, K)

                    if length(_right) >= length(_left)
                        push!(n_next, n << 1)
                        push!(n_next, n << 1 + 1)
                    else
                        push!(n_next, n << 1 + 1)
                        push!(n_next, n << 1)
                    end
                    popfirst!(n_next)
                    # println("n_next split post: ", n, " | ", n_next)
                end
            end
        end
        n_current = copy(n_next)
        depth += 1
    end # end of loop over active ids for a given depth
    return nothing
end

"""
    fit_evotree(params, X_train, Y_train;
        X_eval=nothing, Y_eval=nothing,
        early_stopping_rounds=9999,
        print_every_n=9999,
        verbosity=1)

Main training function. Performorms model fitting given configuration `params`, `X_train`, `Y_train` input data.

# Arguments

- `params::EvoTypes`: configuration info providing hyper-paramters. `EvoTypes` comprises EvoTreeRegressor, EvoTreeClassifier, EvoTreeCount or EvoTreeGaussian
- `X_train::Matrix`: training data of size `[#observations, #features]`.
- `Y_train::Vector`: vector of train targets of length `#observations`.

# Keyword arguments

- `X_eval::Matrix`: training data of size `[#observations, #features]`.
- `Y_eval::Vector`: vector of train targets of length `#observations`.
- `early_stopping_rounds::Integer`: number of consecutive rounds without metric improvement after which fitting in stopped.
- `print_every_n`: sets at which frequency logging info should be printed.
- `verbosity`: set to 1 to print logging info during training.
"""
function fit_evotree(params, X_train, Y_train;
    X_eval=nothing, Y_eval=nothing,
    early_stopping_rounds=9999,
    print_every_n=9999,
    verbosity=1)

    # initialize metric
    iter_since_best = 0
    if params.metric != :none
        metric_track = Metric()
        metric_best = Metric()
    end

    nrounds_max = params.nrounds
    params.nrounds = 0

    if params.device == "gpu"
        model, cache = init_evotree_gpu(params, X_train, Y_train)
        if params.metric != :none && !isnothing(X_eval)
            X_eval = CuArray(eltype(cache.Y).(X_eval))
            Y_eval = CuArray(eltype(cache.Y).(Y_eval))
            pred_eval = predict(params.loss, model.trees[1], X_eval, model.K)
            eval_vec = CUDA.zeros(eltype(cache.pred), size(Y_eval, 1))
        elseif params.metric != :none
            eval_vec = CUDA.zeros(eltype(cache.pred), size(Y_train, 1))
        end
    else
        model, cache = init_evotree(params, X_train, Y_train)
        if params.metric != :none && !isnothing(X_eval)
            pred_eval = predict(params.loss, model.trees[1], X_eval, model.K)
            Y_eval = convert.(eltype(cache.Y), Y_eval)
        end
    end

    while model.params.nrounds < nrounds_max && iter_since_best < early_stopping_rounds
        model.params.nrounds += 1
        grow_evotree!(model, cache)
        # callback function
        if params.metric != :none
            if X_eval !== nothing
                predict!(params.loss, pred_eval, model.trees[model.params.nrounds + 1], X_eval, model.K)
                if params.device == "gpu"
                    metric_track.metric = eval_metric(Val{params.metric}(), eval_vec, pred_eval, Y_eval, params.α)
                else
                    metric_track.metric = eval_metric(Val{params.metric}(), pred_eval, Y_eval, params.α)
                end
            else
                if params.device == "gpu"
                    # println("mean(pred_eval): ", mean(cache.pred))
                    metric_track.metric = eval_metric(Val{params.metric}(), eval_vec, cache.pred, cache.Y, params.α)
                else
                    metric_track.metric = eval_metric(Val{params.metric}(), cache.pred, cache.Y, params.α)
                end
            end
            if metric_track.metric < metric_best.metric
                metric_best.metric = metric_track.metric
                metric_best.iter =  model.params.nrounds
                iter_since_best = 0
            else
                iter_since_best += 1
            end
            if mod(model.params.nrounds, print_every_n) == 0 && verbosity > 0
                display(string("iter:", model.params.nrounds, ", eval: ", metric_track.metric))
            end
        end # end of callback
    end
    if params.metric != :none
        model.metric.iter = metric_best.iter
        model.metric.metric = metric_best.metric
    end
    params.nrounds = nrounds_max
    return model
end


text:
# Implementation of Atomistic API with NBodySimulator.jl

include("nbodysimulator/nbodysimulator_base.jl")
include("nbodysimulator/nbodysimulator_simulator.jl")
include("nbodysimulator/nbodysimulator_result.jl")


text:
module TestBijection

include("preamble.jl")

@testset begin
    l = (@lens _.y[2]) ∘ converting(fromfield = y -> y - 1, tofield = x -> x + 1)
    obj = (x=0, y=(1, 2, 3))
    @test get(obj, l) == 1
    @test set(obj, l, 10).y[2] == 11
end

end  # module


text:
using Conda
using Compat
using PyCall

backend = "theano"

try
    keras = pyimport("keras")
    if VersionNumber(keras[:__version__]) >= v"2.0.2"
        Compat.@info("Using Keras $(keras[:__version__]) -> $(keras[:__path__])")
        global backend = keras[:backend][:backend]()
    else
        Compat.@error("Invalid Keras version ($(keras[:__version__]))")
    end
catch _
    Compat.@info("No valid (>=2.0.2) install of keras found.")
    Compat.@info("Installing Keras via Conda.jl...")
    Conda.add_channel("conda-forge")
    Conda.add("keras==2.0.2")

    if Compat.Sys.iswindows()
        Compat.@warn("Tensorflow only supports python 3.5 on windows. Using the included theano backend.")
    else
        Compat.@info("Installing Tensorflow (1.0.0) via Conda.jl")
        Conda.add("tensorflow==1.0.0")
        global backend = "tensorflow"
    end
end

keras_path = joinpath(homedir(), ".keras")
config_path = joinpath(keras_path, "keras.json")
default_settings = """
{
    "image_dim_ordering": "tf",
    "epsilon": 1e-07,
    "floatx": "float32",
    "backend": "$backend"
}
"""

if !ispath(config_path)
    if !ispath(keras_path)
        mkdir(keras_path)
    end

    Compat.@info("Writing default config to $config_path")

    open(config_path, "w+") do fstream
        write(fstream, default_settings)
    end
end


text:
struct Polygon{PointType <: AbstractPoint}
    point_array::Vector{PointType} # ccw order is assumed, first point != last point
end

function area(polygon)
    L = length(polygon.point_array)
    if L <= 2
        return 0.0
    end
    v1 = polygon.point_array[1]
    a = zero(typeof(v1.x))
    for kk = 1:L
        ii = kk
        jj = util.next_index(ii)

        v1 = polygon.point_array[ii]
        v2 = polygon.point_array[jj]
        a_12 = cross2(v1.x, v1.y, v2.x, v2.y)
        a = a + a_12
    end
    return a * 0.5
end

function create_polygon(x, y)
    # ccw order is assumed
    point_array = map(Point, x, y)
    return Polygon(point_array)
end

function create_polygon(::Type{PointType}, x, y) where {PointType <: AbstractPoint}
    # ccw order is assumed
    point_array = map(PointType, x, y)
    return Polygon(point_array)
end

function n_line_segments(poly)
    # closed polygon is assumed
    return length(poly.point_array)
end

# function resample_line_segments(poly, N_intermediate_points)
#
# end

function discretize_for_plot(poly::Polygon, closed_curve::Bool = true)
    return discretize_for_plot(poly.point_array, closed_curve)
end

function is_diagonal_helper(polygon, a_idx, b_idx)
    L = length(polygon.point_array)
    a = polygon.point_array[a_idx]
    b = polygon.point_array[b_idx]
    for kk = 1:L
        ii = kk
        jj = util.next_index(ii, L)

        v1 = polygon.point_array[ii]
        v2 = polygon.point_array[jj]

        if (a_idx != ii) && (a_idx != jj) &&
            (b_idx != ii) && (b_idx != jj) &&
            segment_intersection(a, b, v1, v2)
            return false
        end
    end
    return true
end

function is_in_cone(polygon, a_idx::Int, b_idx::Int)::Bool
    # a is the base point
    # b is the other point
    L = length(polygon.point_array)
    a_1_idx = util.prev_index(a_idx, L)
    a_3_idx = util.next_index(a_idx, L)

    b = polygon.point_array[b_idx]
    a = polygon.point_array[a_idx]
    a_1 = polygon.point_array[a_1_idx]
    a_3 = polygon.point_array[a_3_idx]

    if is_left_on(a, a_3, a_1)
        return is_left_on(a, b, a_1) && is_left_on(b, a, a_3)
    else
        return !( is_left_on(a, b, a_3) && is_left_on(b, a, a_1) )
    end
end

function is_diagonal(polygon, a_idx, b_idx)
    return  is_in_cone(polygon, a_idx, b_idx) &&
            is_in_cone(polygon, b_idx, a_idx) &&
            is_diagonal_helper(polygon, a_idx, b_idx)
end

function __init_ear(polygon)
    L = length(polygon.point_array)
    bool = falses(L)
    for kk = 1:L
        aa = util.prev_index(kk, L)
        bb = kk
        cc = util.next_index(kk, L)
        bool[kk] = is_diagonal(polygon, aa, cc)
    end
    return bool
end

mutable struct PolygonTriangle
    index_1::Int
    index_2::Int
    index_3::Int
end

function triangulate(polygon)
    L = length(polygon.point_array)
    if L <= 3
        return [PolygonTriangle(1, 2, 3)]
    end
    ears = __init_ear(polygon)
    if L == 4
        if all(ears)
            t1 = PolygonTriangle(1, 2, 3)
            t2 = PolygonTriangle(3, 4, 1)
            return [t1, t2]
        else
            # concave case
            i1 = 0
            i2 = 0
            for kk = 1:4
                if ears[kk] == false && i1 == 0
                    i1 = kk
                elseif ears[kk] == false && i2 == 0
                    i2 = kk
                end
            end
            if i2 != i1 + 2
                error("Something went wrong")
            end
            k1 = i1
            k2 = util.next_index(k1, 4)
            k3 = util.next_index(k2, 4)
            k4 = util.next_index(k3, 4)
            t1 = PolygonTriangle(k1, k2, k3)
            t2 = PolygonTriangle(k3, k4, k1)
            return [t1, t2]
        end

    end # L 4

    poly_list = datastructs.array_to_linked_list(polygon.point_array)
    ears_list = datastructs.array_to_linked_list(Vector(ears))
    idx_list = datastructs.array_to_linked_list(collect(1:length(ears)))
    poly_tri_array = Vector{PolygonTriangle}(undef, L-2)
    ptr = 1
    n = L
    kk = 1
    first_node = poly_list.node
    while n > 3
        v2 = poly_list.node
        v2next = v2.next
        ear_node  = ears_list.node
        idx_node = idx_list.node
        while true
            if ear_node.data
                v1 = v2.prev
                v3 = v2.next
                v0 = v1.prev
                v4 = v2.next

                # v1 - v3 is a diagonal
                t = PolygonTriangle(
                    idx_node.prev.data,
                    idx_node.data,
                    idx_node.next.data)
                poly_tri_array[ptr] = t
                ptr = ptr + 1

                datastructs.remove_element!(idx_list, idx_node)
                datastructs.remove_element!(ears_list, ear_node)
                datastructs.remove_element!(poly_list, v2)
                n -= 1
                first_node = v3
                break
            end
            v2 = v2next
            if first_node == v2
                break
            end
        end
        return poly_tri_array
    end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function find_close(pos, roots, eps)
    pos = [pos[1], abs(pos[2])] # Root only has positive conjugate
    closest = Root(Point2f0(Inf, Inf), false, false)
    for root in roots
        if sum(abs2, root.pos .- pos) < sum(abs2, closest.pos .- pos)
            closest = root
        end
    end
    if all(abs.(closest.pos .- pos) .< eps)
        return closest
    end
end

""" Root only contain the conjugate with positive imaginary part
"""
mutable struct Root
    pos::Point2f0
    pole::Bool
    selected::Bool
end

function get_poles(roots)
    a = Point2f0[]
    for node in roots
        if node.pole
            push!(a, node.pos)
            if node.pos[2] != 0
                push!(a, node.pos .* [1, -1])
            end
        end
    end
    return a
end

function to_points(roots)
    a = Point2f0[]
    for node in roots
        push!(a, node.pos)
        if node.pos[2] != 0
            push!(a, node.pos .* [1, -1])
        end
    end
    return a
end

function get_zeros(roots)
    a = Point2f0[]
    for node in roots
        if !node.pole
            push!(a, node.pos)
            if node.pos[2] != 0
                push!(a, node.pos .* [1, -1])
            end
        end
    end
    return a
end

function get_selected(roots)
    for i in eachindex(roots)
        if roots[i].selected
            if roots[i].pos[2] == 0
                return i, Point2f0[roots[i].pos], roots[i].pole ? '+' : 'o'
            else
                return i, Point2f0[roots[i].pos, [1, -1] .* roots[i].pos], roots[i].pole ? '+' : 'o'
            end
        end
    end
    return 0, Point2f0[], '+'
end

function unselect_all!(roots)
    for root in roots[]
        root.selected = false
    end
end

max_stars_repo_name: albheim/PoleExploration

max_issues_repo_name: albheim/PoleExploration

max_forks_repo_name: albheim/PoleExploration


text:
module Acquisitions
using AxisArrays, Images
export Acquisition, SIMAcquisition, times, FluorophoreMarkers

using ..FluorophoreMarkers

abstract type AbstractAcquisition end

include("./acquisition.jl")
include("./simacquisition.jl")

end


text:
using Dates
using LinearAlgebra
using SparseArrays
include("../src/wildart.jl")

# using Logging
# debug_logger = ConsoleLogger(stderr, Logging.Debug)
# global_logger(debug_logger)

# X = [1 2 2 7;
#      1 2 6 1;
#      4 7 1 10]
# Y = copy(X)
# snf(Y)[3]
# using SmithNormalForm
# smith(Y)

include("../src/boundary_map_Alex.jl")
bigY = bound_matrix_4(9,4,4)
result = snf(bigY)

# 170141142824549539861647102733501356222 < typemax(Int128)
# @time for i in 1:10000
#     rand(BigInt) ÷ rand(BigInt)
# end

# M = rand(0:9,3,3)

# convert.(Real, M)

# x = -9223372036854775808; typeof(x)
# y = Int128(x); typeof(y)
# div(x,-1)
# div(y,-1)

# div(92233720368547758080,-1)

# typemax(BigInt)
# -9223372036854775808
# typemax(-9223372036854775808)

# convert.(Int128,rand(0:9, 4,4))


text:
function test_arctan_all()
    print(rpad("Math.Series.arctan[h](1/x)...", 50, ' '))

    arctan = Brainstorm._Math._Series._Arctan.arctan
    arctanh = Brainstorm._Math._Series._Arctan.arctanh

    setprecision(400_000) do
        @test all([arctan(d, x)  - trunc(BigInt,  atan(1//big(x))*(big(10)^d)) for x in (2, 3, 5, 17), d in (1, 10, 100, 1_000, 10_000, 100_000)] .== 0)
        @test all([arctanh(d, x) - trunc(BigInt, atanh(1//big(x))*(big(10)^d)) for x in (2, 3, 5, 17), d in (1, 10, 100, 1_000, 10_000, 100_000)] .== 0)
    end

    println("PASS")
end


text:
using KhepriTikZ
using Test

@testset "KhepriTikZ.jl" begin
    # Write your tests here.
end


text:
struct ParamInfo
    name::Symbol
    typ::Type
end

abstract type Node end
abstract type ExprNode <: Node end
abstract type ValueNode <: Node end

# utility

function replace_input_symbols(symbols::Set{Symbol}, value, trace::Symbol)
    value
end

function replace_input_symbols(symbols::Set{Symbol}, expr::Expr, trace::Symbol)
    args = [replace_input_symbols(symbols, arg, trace) for arg in expr.args]
    Expr(expr.head, args...)
end

function replace_input_symbols(symbols::Set{Symbol}, symbol::Symbol, trace::Symbol)
    if symbol in symbols
        return Expr(:(.), trace, QuoteNode(value_field(symbol)))
    else
        return Expr(:(.), :Main, QuoteNode(symbol))
    end
end

# statement nodes

struct ArgsChangeNode <: ExprNode
    output::ValueNode
end

function Base.print(io::IO, node::ArgsChangeNode)
    write(io, "@argschange")
end

parents(::ArgsChangeNode) = []

struct AddrDistNode{T} <: ExprNode
    input_nodes::Vector{ValueNode}
    #output::Union{ValueNode,Nothing} # TODO always have an output nod
    output::ValueNode # TODO always have an output nod
    dist::Distribution{T}
    address::Symbol
    line::LineNumberNode
end

function Base.print(io::IO, node::AddrDistNode)
    write(io, "@addr $(node.dist) $(node.address)")
end

parents(node::AddrDistNode) = node.input_nodes

struct AddrGeneratorNode{T,U} <: ExprNode
    input_nodes::Vector{ValueNode}
    #output::Union{ValueNode,Nothing} # TODO always have an output node
    output::ValueNode # TODO always have an output node
    gen::Generator{T,U}
    address::Symbol
    change_node::ValueNode
    line::LineNumberNode
end

function Base.print(io::IO, node::AddrGeneratorNode)
    write(io, "@addr $(node.gen) $(node.address)")
end

parents(node::AddrGeneratorNode) = vcat([node.change_node], node.input_nodes)

struct AddrChangeNode <: ExprNode
    address::Symbol
    addr_node::Union{AddrDistNode,AddrGeneratorNode}
    output::ValueNode
end

function Base.print(io::IO, node::AddrChangeNode)
    write(io, "@change $(node.address)")
end

# so that the @change will be placed after its @addr node
parents(node::AddrChangeNode) = [node.addr_node]

struct JuliaNode <: ExprNode
    input_nodes::Vector{ValueNode}
    output::ValueNode
    expr_or_value::Any
    line::Union{LineNumberNode,Nothing}
end

function Base.print(io::IO, node::JuliaNode)
    write(io, "Julia $(gensym())")
end

parents(node::JuliaNode) = node.input_nodes

function expr_read_from_trace(node::JuliaNode, trace::Symbol)
    input_symbols = Set{Symbol}([node.name for node in node.input_nodes])
    replace_input_symbols(input_symbols, node.expr_or_value, trace)
end

# value nodes
struct ArgumentValueNode <: ValueNode
    name::Symbol
    typ::Type
end

function Base.print(io::IO, node::ArgumentValueNode)
    write(io, "Arg $(node.name)")
end

parents(node::ArgumentValueNode) = []
get_type(node::ArgumentValueNode) = node.typ

mutable struct ExprValueNode <: ValueNode
    name::Symbol
    typ::Type
    finished::Bool
    source::ExprNode
    function ExprValueNode(name, typ)
        new(name, typ, false)
    end
end

function finish!(node::ExprValueNode, source::ExprNode)
    @assert !node.finished
    node.finished = true
    node.source = source
end

function Base.print(io::IO, node::ExprValueNode)
    write(io, "$(node.name)  $(node.typ)")
end

parents(node::ExprValueNode) = [node.source]
get_type(node::ExprValueNode) = node.typ

struct ParamValueNode <: ValueNode
    param::ParamInfo
    name::Symbol
end

ParamValueNode(param::ParamInfo) = ParamValueNode(param, param.name)

function Base.print(io::IO, node::ParamValueNode)
    write(io, "Param $(node.param.name)  $(node.param.typ)")
end

parents(node::ParamValueNode) = []
get_type(node::ParamValueNode) = node.param.typ

mutable struct BasicBlockIR
    output_ad::Bool
    args_ad::Tuple
    arg_nodes::Vector{ArgumentValueNode}
    params::Vector{ParamInfo}
    value_nodes::Dict{Symbol, ValueNode}
    addr_dist_nodes::Dict{Symbol, AddrDistNode}
    addr_gen_nodes::Dict{Symbol, AddrGeneratorNode}
    all_nodes::Set{Node}
    output_node::Union{ValueNode,Nothing}
    retchange_node::Union{ValueNode,Nothing}
    args_change_node::Union{ArgsChangeNode,Nothing}
    addr_change_nodes::Set{AddrChangeNode}
    generator_input_change_nodes::Set{ValueNode}
    finished::Bool

    # nodes that depend on incremental values (@argchange() and @change())
    # that aren't part of the actual function. used to check that there is no
    # data flow from incremental values to the actual function's observable
    # behavior (trace or return value).
    incremental_nodes::Set{ValueNode}

    expr_nodes_sorted::Vector{ExprNode}

    function BasicBlockIR(output_ad::Bool, args_ad::Tuple)
        arg_nodes = Vector{ArgumentValueNode}()
        params = Vector{ParamInfo}()
        value_nodes = Dict{Symbol,ValueNode}()
        addr_dist_nodes = Dict{Symbol, AddrDistNode}()
        addr_gen_nodes = Dict{Symbol, AddrGeneratorNode}()
        all_nodes = Set{Node}()
        output_node = nothing
        retchange_node = nothing
        args_change_node = nothing
        addr_change_nodes = Set{AddrChangeNode}()
        generator_input_change_nodes = Set{ValueNode}()
        finished = false
        incremental_nodes = Set{ValueNode}()
        new(output_ad, args_ad,
            arg_nodes, params, value_nodes, addr_dist_nodes,
                     addr_gen_nodes, all_nodes, output_node, retchange_node,
                     args_change_node, addr_change_nodes,
                     generator_input_change_nodes,
                     finished, incremental_nodes)
    end
end

function toposort_visit!(node, permanent_marked, temporary_marked,
                         unmarked, reverse_list)
    if node in permanent_marked
        return
    end
    if node in temporary_marked
        error("cycle found")
    end
    push!(temporary_marked, node)
    for parent in parents(node)
        toposort_visit!(parent, permanent_marked, temporary_marked, unmarked, reverse_list)
    end
    push!(permanent_marked, node)
    push!(reverse_list, node)
end

function toposort_expr_nodes(nodes)
    permanent_marked = Set{Node}()
    temporary_marked = Set{Node}()
    unmarked = Set{Node}(nodes)
    reverse_list = Vector{Node}()
    while !isempty(unmarked)
        node = pop!(unmarked)
        toposort_visit!(node, permanent_marked, temporary_marked, unmarked, reverse_list)
    end
    priorities = Dict(node => i for (i, node) in enumerate(reverse_list))
    expr_nodes = collect(
        Iterators.filter((node) -> isa(node, ExprNode), keys(priorities)))
    convert(Vector{ExprNode}, sort(expr_nodes, by=(node) -> priorities[node]))
end

function finish!(ir::BasicBlockIR)
    ir.finished = true
    ir.expr_nodes_sorted = toposort_expr_nodes(ir.all_nodes)
    if ir.output_node === nothing && ir.output_ad
        error("Compiled Gen function marked for output @ad but there is no return statement")
    end
end

# NOTE: currently we don't properly support list comprehensions and closures
# and let (anything that defines a new environment)
function add_argument!(ir::BasicBlockIR, name::Symbol, typ)
    @assert !ir.finished
    node = ArgumentValueNode(name, typ)
    push!(ir.arg_nodes, node)
    push!(ir.all_nodes, node)
    ir.value_nodes[name] = node
end

function add_param!(ir::BasicBlockIR, name::Symbol, typ)
    @assert !ir.finished
    if !(typ <: AbstractFloat) || (typ <: Array{T} where {T <: AbstractFloat})
        error("@param type $typ is not <: AbstractFloat or <: Array{T} where {T <: AbstractFloat}")
    end
    param_info = ParamInfo(name, typ)
    push!(ir.params, param_info)
    node = ParamValueNode(param_info)
    if haskey(ir.value_nodes, name)
        error("Name $name already used")
    end
    ir.value_nodes[name] = node
    push!(ir.all_nodes, node)
end

function resolve_symbols_set(ir::BasicBlockIR, symbol::Symbol)
    symbols = Set{Symbol}()
    if haskey(ir.value_nodes, symbol)
        push!(symbols, symbol)
    end
    symbols
end

resolve_symbols_set(ir::BasicBlockIR, expr) = Set{Symbol}()

function resolve_symbols_set(ir::BasicBlockIR, expr::Expr)
    symbols = Set{Symbol}()
    if expr.head == :(.)
        union!(symbols, resolve_symbols_set(ir, expr.args[1]))
    elseif expr.head != :line
        for operand in expr.args
            union!(symbols, resolve_symbols_set(ir, operand))
        end
    end
    symbols
end

function _get_input_node!(ir::BasicBlockIR, expr, typ)
    name = gensym()
    value_node = ExprValueNode(name, typ)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, value_node)
    input_nodes = map((sym) -> ir.value_nodes[sym],
                      collect(resolve_symbols_set(ir, expr)))
    expr_node = JuliaNode(input_nodes, value_node, expr, nothing)
    finish!(value_node, expr_node)
    push!(ir.all_nodes, expr_node)
    if any([node in ir.incremental_nodes for node in input_nodes])
        push!(ir.incremental_nodes, value_node)
    end
    value_node
end

function _get_input_node!(ir::BasicBlockIR, name::Symbol, typ)
    if haskey(ir.value_nodes, name)
        node = ir.value_nodes[name]
        #@assert get_type(node) == typ # TODO
        node
    else
        # this is a symbol defined in the top-level environment
        value_node = ExprValueNode(name, typ)
        expr_node = JuliaNode(ValueNode[], value_node, name, nothing)
        ir.value_nodes[name] = value_node
        push!(ir.all_nodes, value_node)
        finish!(value_node, expr_node)
        push!(ir.all_nodes, expr_node)
        value_node
    end
end

function incremental_dependency_error(addr)
    error("@addr argument cannot depend on @change or @argschange statements (address $addr)")
end

function add_addr!(ir::BasicBlockIR, addr::Symbol, line::LineNumberNode,
                   dist::Distribution{T}, args::Vector,
                   typ::Type, name::Symbol) where {T}
    @assert !ir.finished
    types = get_static_argument_types(dist)
    input_nodes = ValueNode[
        _get_input_node!(ir, expr, typ) for (expr, typ) in zip(args, types)]
    if any([node in ir.incremental_nodes for node in input_nodes])
        incremental_dependency_error(addr)
    end
    value_node = ExprValueNode(name, typ)
    expr_node = AddrDistNode(input_nodes, value_node, dist, addr, line)
    finish!(value_node, expr_node)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, expr_node)
    ir.addr_dist_nodes[addr] = expr_node
    push!(ir.all_nodes, value_node)
    nothing
end

function add_addr!(ir::BasicBlockIR, addr::Symbol, line::LineNumberNode,
                   dist::Distribution{T}, args::Vector) where {T}
    add_addr!(ir, addr, line, dist, args, get_return_type(dist), gensym())
end

function add_addr!(ir::BasicBlockIR, addr::Symbol, line::LineNumberNode,
                   gen::Generator{T,U}, args::Vector,
                   typ::Type, name::Symbol, change_expr) where {T,U}
    @assert !ir.finished
    types = get_static_argument_types(gen)
    # NOTE: these could be abstract types....
    input_nodes = ValueNode[
        _get_input_node!(ir, expr, typ) for (expr, typ) in zip(args, types)]
    if any([node in ir.incremental_nodes for node in input_nodes])
        incremental_dependency_error(addr)
    end
    value_node = ExprValueNode(name, typ)
    # TODO make it not Some{}
    change_node = _get_input_node!(ir, change_expr, Union{Some{get_change_type(gen)}, Nothing})
    push!(ir.generator_input_change_nodes, change_node)
    expr_node = AddrGeneratorNode(input_nodes, value_node, gen, addr, change_node, line)
    finish!(value_node, expr_node)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, expr_node)
    ir.addr_gen_nodes[addr] = expr_node
    push!(ir.all_nodes, value_node)
    nothing
end

function add_addr!(ir::BasicBlockIR, addr::Symbol, line::LineNumberNode,
                   gen::Generator{T,U}, args::Vector, change_expr) where {T,U}
    typ = get_return_type(gen)
    add_addr!(ir, addr, line, gen, args, typ, gensym(), change_expr)
end

function add_julia!(ir::BasicBlockIR, expr, typ, name::Symbol, line::LineNumberNode)
    @assert !ir.finished
    value_node = ExprValueNode(name, typ)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, value_node)
    input_nodes = map((sym) -> ir.value_nodes[sym],
                      collect(resolve_symbols_set(ir, expr)))
    expr_node = JuliaNode(input_nodes, value_node, expr, line)
    finish!(value_node, expr_node)
    push!(ir.all_nodes, expr_node)
    if any([node in ir.incremental_nodes for node in input_nodes])
        push!(ir.incremental_nodes, value_node)
    end
    value_node
end

function add_argschange!(ir::BasicBlockIR, typ, name::Symbol)
    @assert !ir.finished
    if ir.args_change_node !== nothing
        error("@argschange can only be called once")
    end
    value_node = ExprValueNode(name, typ)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, value_node)
    expr_node = ArgsChangeNode(value_node)
    finish!(value_node, expr_node)
    push!(ir.all_nodes, expr_node)
    push!(ir.incremental_nodes, value_node)
    ir.args_change_node = expr_node
    nothing
end

function add_change!(ir::BasicBlockIR, addr::Symbol, typ, name::Symbol)
    @assert !ir.finished
    value_node = ExprValueNode(name, typ)
    ir.value_nodes[name] = value_node
    push!(ir.all_nodes, value_node)
    if haskey(ir.addr_dist_nodes, addr)
        addr_node = ir.addr_dist_nodes[addr]
    elseif haskey(ir.addr_gen_nodes, addr)
        addr_node = ir.addr_gen_nodes[addr]
    else
        error("@change must occur after an @addr. Address: $addr ")
    end
    expr_node = AddrChangeNode(addr, addr_node, value_node)
    finish!(value_node, expr_node)
    push!(ir.all_nodes, expr_node)
    push!(ir.incremental_nodes, value_node)
    push!(ir.addr_change_nodes, expr_node)
    nothing
end

function set_return!(ir::BasicBlockIR, name::Symbol)
    if ir.output_node !== nothing
        error("Basic block can only have one return statement, found a second: $name")
    end
    value_node = ir.value_nodes[name]
    if value_node in ir.incremental_nodes
        error("Return value cannot depend on @change or @argschange statements")
    end
    ir.output_node = value_node
end

function set_return!(ir::BasicBlockIR, expr::Expr, typ::Type)
    if ir.output_node !== nothing
        error("Basic block can only have one return statement, found a second: $name")
    end
    ir.output_node = _get_input_node!(ir, expr, typ)
end

function set_retchange!(ir::BasicBlockIR, name::Symbol)
    if ir.retchange_node !== nothing
        error("Basic block can only have one @retchange statement, found a second: $name")
    end
    value_node = ir.value_nodes[name]
    ir.retchange_node = value_node
end

function set_retchange!(ir::BasicBlockIR, expr::Expr, typ::Type)
    if ir.retchange_node !== nothing
        error("Basic block can only have one @retchange statement, found a second: $name")
    end
    ir.retchange_node = _get_input_node!(ir, expr, typ)
end

# some helper functions

has_output(node::ExprNode) = true
has_output(node::Union{AddrDistNode,AddrGeneratorNode}) = node.output !== nothing

function get_value_info(node::Union{AddrDistNode,AddrGeneratorNode})
    value_node::ValueNode = node.output
    (value_node.typ, value_field(value_node))
end

function get_value_info(node::ExprNode)
    value_node::ValueNode = node.output
    (value_node.typ, value_field(value_node))
end

function get_args(trace::Symbol, node::ExprNode)
    map(input_node -> Expr(:(.), trace, QuoteNode(value_field(input_node))), node.input_nodes)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
######### multiple distributions  ###########
using Distributions

VECTOR_SIZE = 1000
v = ones(VECTOR_SIZE)

function bench(ex::Expr)  # ex = :(Weibull(1, 1))
      # to get valid intial values, use the mean of the distribution
      exactMean = eval( :(mean($ex)) )
      # some distribs (Cauchy) have no defined mean
      exactMean = isfinite(exactMean) ? exactMean : 1.0

      # exactStd = eval( :(std($ex)) )
      # exactStd = isfinite(exactStd) ? exactStd : 1.0

	mex = Expr(:block, :(y = x * v ; y ~ $ex))

	m = model(mex, x=exactMean)
	mg = model(mex, gradient=true, x=exactMean)

	name = "$ex on vector of $VECTOR_SIZE"

	[    benchmark( ()-> m.eval(m.init),     "loglik eval", name, 100) ;
	     benchmark( () -> mg.evalallg(mg.init), "loglik and gradient eval", name, 100) ;
	     benchmark( () -> run(m * RWM(0.1), steps=100), "100 RWM steps", name, 10) ]
end

distribs = [:(Normal(1, 1)),
            :(Normal(3, 12)),
            :(Weibull(1, 1)),
            :(Weibull(3, 1)),
            :(Uniform(0, 2)),
            :(TDist(2.2)),
            :(TDist(4)),
            :(Beta(1,2)),
            :(Beta(3,2)),
            :(Gamma(1,2)),
            :(Gamma(3,0.2)),
            :(Cauchy(0,1)),
            :(Cauchy(-1,0.2)),
            :(Exponential(3)),
            :(Exponential(0.2)),
            :(LogNormal(-1, 1)),
            :(LogNormal(2, 0.1))]

res = mapreduce(bench, vcat, distribs)


text:
# # Tutorial 🐧

# This is a gentle and lighthearted tutorial on how to use tools from AlgebraOfGraphics,
# using as example dataset a collection of measurements on penguins[^1]. See
# the [Palmer penguins website](https://allisonhorst.github.io/palmerpenguins/index.html)
# for more information.
#
# [^1]: Gorman KB, Williams TD, Fraser WR (2014) Ecological Sexual Dimorphism and Environmental Variability within a Community of Antarctic Penguins (Genus Pygoscelis). PLoS ONE 9(3): e90081. [DOI](https://doi.org/10.1371/journal.pone.0090081)

# To follow along this tutorial, you will need to install a few packages.
# All the required packages can be installed with the following command.
#
# ```julia
# julia> import Pkg; Pkg.add(["AlgebraOfGraphics", "CairoMakie", "DataFrames", "LIBSVM", "PalmerPenguins"])
# ```
#
# After the above command completes, we are ready to go.

using PalmerPenguins, DataFrames

penguins = dropmissing(DataFrame(PalmerPenguins.load()))
first(penguins, 6)

# ## Frequency plots
#
# Let us start by getting a rough idea of how the data is distributed.
#
# !!! note
#     Due to julia's compilation model, the first plot may take a while to appear.

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

axis = (width = 225, height = 225)
penguin_frequency = data(penguins) * frequency() * mapping(:species)

draw(penguin_frequency; axis)

# ### Small intermezzo: saving the plot
#
# If you are working in an interactive enviroment with inline plotting support,
# such VSCode or Pluto.jl, the above should have displayed a bar plot.
# If you are working directly in the console, you can simply save the plot and
# inspect it in the file explorer.
#
# ```julia
# fg = draw(penguin_frequency; axis)
# save("figure.png", fg, px_per_unit = 3) # save high-resolution png
# ```
#
# ### Styling by categorical variables
#
# Next, let us see whether the distribution is the same across islands.

plt = penguin_frequency * mapping(color = :island)
draw(plt; axis)

# Ups! The bars are in the same spot and are hiding each other. We need to specify
# how we want to fix this. Bars can either `dodge` each other, or be `stack`ed on top
# of each other.

plt = penguin_frequency * mapping(color = :island, dodge = :island)
draw(plt; axis)

# This is our first finding. `Adelie` is the only species of penguins that can be
# found on all three islands. To be able to see both which species is more numerous
# and how different species are distributed across islands in a unique plot,
# we could have used `stack`.

plt = penguin_frequency * mapping(color = :island, stack = :island)
draw(plt; axis)

# ## Correlating two variables
#
# Now that we have understood the distribution of these three penguin species, we can
# start analyzing their features.

penguin_bill = data(penguins) * mapping(:bill_length_mm, :bill_depth_mm)
draw(penguin_bill; axis)

# We would actually prefer to visualize these measures in centimeters, and to have
# cleaner axes labels. As we want this setting to be preserved in all of our `bill`
# visualizations, let us save it in the variable `penguin_bill`, to be reused in
# subsequent plots.

penguin_bill = data(penguins) * mapping(
    :bill_length_mm => (t -> t / 10) => "bill length (cm)",
    :bill_depth_mm => (t -> t / 10) => "bill depth (cm)",
)
draw(penguin_bill; axis)

# Much better! Note the parentheses around the function `t -> t / 10`. They are
# necessary to specify that the function maps `t` to `t / 10`, and not to
# `t / 10 => "bill length (cm)"`.

# There does not seem to be a strong correlation between the two dimensions, which
# is odd. Maybe dividing the data by species will help.

plt = penguin_bill * mapping(color = :species)
draw(plt; axis)

# Ha! Within each species, penguins with a longer bill also have a deeper bill.
# We can confirm that with a linear regression

plt = penguin_bill * linear() * mapping(color = :species)
draw(plt; axis)

# This unfortunately no longer shows our data!
# We can use `+` to plot both things on top of each other:

plt = penguin_bill * linear() * mapping(color = :species) + penguin_bill * mapping(color = :species)
draw(plt; axis)

# Note that the above expression seems a bit redundant, as we wrote the same thing twice.
# We can "factor it out" as follows

plt = penguin_bill * (linear() + mapping()) * mapping(color = :species)
draw(plt; axis)

# where `mapping()` is a neutral multiplicative element.
# Of course, the above could be refactored as

layers = linear() + mapping()
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# We could actually take advantage of the spare `mapping()` and use it to pass some
# extra info to the scatter, while still using all the species members to compute
# the linear fit.

layers = linear() + mapping(marker = :sex)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# This plot is getting a little bit crowded. We could instead show female and
# male penguins in separate subplots.

layers = linear() + mapping(col = :sex)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# See how both plots show the same fit, because the `sex` mapping is not applied
# to `linear()`. The following on the other hand produces a separate fit for
# males and females:

layers = linear() + mapping()
plt = penguin_bill * layers * mapping(color = :species, col = :sex)
draw(plt; axis)

# ## Smooth density plots
#
# An alternative approach to understanding how two variables interact is to consider
# their joint probability density distribution (pdf).

using AlgebraOfGraphics: density
plt = penguin_bill * density(npoints=50) * mapping(col = :species)
draw(plt; axis)

# The default colormap is multi-hue, but it is possible to pass single-hue colormaps as well:

plt *= visual(colormap = :grayC)
draw(plt; axis)

# A `Heatmap` (the default visualization for a 2D density) is a bit unfortunate if
# we want to mark species by color. In that case, one can use `visual` to change
# the default visualization and, optionally, fine tune some arguments.
# In this case, a `Wireframe` with thin lines looks quite nice. (Note that, for the
# time being, we must specify explicitly that we require a 3D axis.)

axis = (type = Axis3, width = 300, height = 300)
layer = density() * visual(Wireframe, linewidth=0.05)
plt = penguin_bill * layer * mapping(color = :species)
draw(plt; axis)

# Of course, a more traditional approach would be to use a `Contour` plot instead:

axis = (width = 225, height = 225)
layer = density() * visual(Contour)
plt = penguin_bill * layer * mapping(color = :species)
draw(plt; axis)

# The data and the linear fit can also be added back to the plot:

layers = density() * visual(Contour) + linear() + mapping()
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# In the case of many layers (contour, density and scatter) it is important to think
# about balance. In the above plot, the markers are quite heavy and can obscure the linear
# fit and the contour lines.
# We can lighten the markers using alpha transparency.

layers = density() * visual(Contour) + linear() + visual(alpha = 0.5)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# ## Correlating three variables
#
# We are now mostly up to speed with `bill` size, but we have not considered how
# it relates to other penguin features, such as their weight.
# For that, a possible approach is to use a continuous color
# on a gradient to denote weight and different marker shapes to denote species.
# Here we use `group` to split the data for the linear regression without adding
# any additional style.

body_mass = :body_mass_g => (t -> t / 1000) => "body mass (kg)"
layers = linear() * mapping(group = :species) + mapping(color = body_mass, marker = :species)
plt = penguin_bill * layers
draw(plt; axis)

# Naturally, within each species, heavier penguins have bigger bills, but perhaps
# counter-intuitively the species with the shallowest bills features the heaviest penguins.
# We could also try and see the interplay of these three variables in a 3D plot.

axis = (type = Axis3, width = 300, height = 300)
plt = penguin_bill * mapping(body_mass, color = :species)
draw(plt; axis)

#

plt = penguin_bill * mapping(body_mass, color = :species, layout = :sex)
draw(plt; axis)

# Note that static 3D plot can be misleading, as they only show one projection
# of 3D data. They are mostly useful when shown interactively.
#
# ## Machine Learning
#
# Finally, let us use Machine Learning techniques to build an automated penguin classifier!
#
# We would like to investigate whether it is possible to predict the species of a penguin
# based on its bill size. To do so, we will use a standard classifier technique
# called [Support-Vector Machine](https://en.wikipedia.org/wiki/Support-vector_machine).
#
# The strategy is quite simple. We split the data into training and testing
# subdatasets. We then train our classifier on the training dataset and use it to
# make predictions on the whole data. We then add the new columns obtained this way
# to the dataset and visually inspect how well the classifier performed in both
# training and testing.

using LIBSVM, Random

## use approximately 80% of penguins for training
Random.seed!(1234) # for reproducibility
N = nrow(penguins)
train = fill(false, N)
perm = randperm(N)
train_idxs = perm[1:floor(Int, 0.8N)]
train[train_idxs] .= true
nothing # hide

## fit model on training data and make predictions on the whole dataset
X = hcat(penguins.bill_length_mm, penguins.bill_depth_mm)
y = penguins.species
model = SVC() # Support-Vector Machine Classifier
fit!(model, X[train, :], y[train])
ŷ = predict(model, X)

## incorporate relevant information in the dataset
penguins.train = train
penguins.predicted_species = ŷ
nothing #hide

# Now, we have all the columns we need to evaluate how well our classifier performed.

axis = (width = 225, height = 225)
dataset =:train => renamer(true => "training", false => "testing") => "Dataset"
accuracy = (:species, :predicted_species) => isequal => "accuracy"
plt = data(penguins) *
    expectation() *
    mapping(:species, accuracy) *
    mapping(col = dataset)
draw(plt; axis)

# That is a bit hard to read, as all values are very close to `1`.
# Let us visualize the error rate instead.

error_rate = (:species, :predicted_species) => !isequal => "error rate"
plt = data(penguins) *
    expectation() *
    mapping(:species, error_rate) *
    mapping(col = dataset)
draw(plt; axis)

# So, mostly our classifier is doing quite well, but there are some mistakes,
# especially among `Chinstrap` penguins. Using *at the same time* the `species` and
# `predicted_species` mappings on different attributes, we can see which penguins
# are problematic.

prediction = :predicted_species => "predicted species"
datalayer = mapping(color = prediction, row = :species, col = dataset)
plt = penguin_bill * datalayer
draw(plt; axis)

# Um, some of the penguins are indeed being misclassified... Let us try to understand why
# by adding an extra layer, which describes the density of the distributions of the three
# species.

pdflayer = density() * visual(Contour, colormap=Reverse(:grays)) * mapping(group = :species)
layers = pdflayer + datalayer
plt = penguin_bill * layers
draw(plt; axis)

# We can conclude that the classifier is doing a reasonable job:
# it is mostly making mistakes on outlier penguins.


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

using Markdown
using InteractiveUtils

# ╔═╡ 4c4bfc78-7224-11ec-211d-f1e60824d35c
begin
	using PlutoUI
	using Symbolics, SymbolicUtils
	using LaTeXStrings,Latexify
	# using Bijections  # could use to specify replacement of variables/names
end

# ╔═╡ b9b5c2e0-537d-4a72-adb9-8474dc8acf98
md"""
# Julia-functions with keyword arguments from symbolic expressions

### (I'm grateful for any advice or improvement suggestions!!!)

## Aim
In this notebook, the (not yet achieved) goal is to eliminate the danger of messing up when copying equations between publications (or LaTeX manuscripts) and Julia-code.

Such mistakes can lead to graphs or simulation results in a paper no longer matching reported equations. This can make papers difficult to understand, reproduce and build upon their results. Parameters used for plots or simulations/optimizations are sometimes not completely specified or reported incorrectly, which can even make it difficult to say if one's attempt to reproduce the results succeeded.

My current idea of how to work towards eliminating this problem is this:

- Use a `Pluto.jl` notebook instead of a paper, or at least give it as a supplement.
- Enter the equations using `Symbolics.jl` (or `ModelingToolkit.jl`, etc.) and render them as LaTeX. This makes it much easier to read and check them then bracketed Julia code.
- Manipulate the equations and expressions to plug in numbers, run simulations or optimizations, and finally plot results.

However, there are several problems on the way towards this goal. Also, being a beginner at using the Julia SciML ecosystem, it's hard to follow best-practices, especially given the fast pace of development, changes and lack of accessible documentation. The *this I focus on right now* are

- dealing with weird variable names (e.g. with superscripts)
- plugging in numbers for variables, plotting equations and preparing for simulation/optimization, in a way that is simple, convenient, stable and helps catch errors

"""

# ╔═╡ 42f6afaa-cf0d-46dd-b1d9-ebaeafc9e5db
md"""
## Contents
- Render equation with complicated variable names
- substitute an equation that is split up into sub-expressions back into one equation
- Plug in numbers into equation using `substitute`
- Compile Julia function (array input) from the expression (using `build_function`)
- Compile Julia function (keyword-argument input) from the expression
"""

# ╔═╡ b580a960-7780-417e-ae0f-3ac59f4a4a30
md"""
## Use non-standard variable names

Suppose you are dealing with symbols that do not map uniquely to Julia code, for example:
"""

# ╔═╡ 0c2701f4-8e93-4fbc-b18c-52e9d1c213f8
begin
	proper_vars = @variables x a b var"{c^{+}_{\alpha}}"

	c_plus_alpha = var"{c^{+}_{\alpha}}" # abbreviation for latex-symbol

	proper_vars
end

# ╔═╡ 6e521fb4-190c-4b68-890b-d39b22537471
md"""
_I'm not sure if it's a good idea to represent such symbols with this hack, i.e., abusing non-standard Julia variable names. This seems to mostly work but leads to problems when the variable is programatically converted to a `Base.Symbol`._

Now imagine you copy an equation with such symbols from a publication and want to keep your names and symbols as close to the original as possible. A toy-example equation:
"""

# ╔═╡ bf6bac04-1957-4aa9-976e-efd76b41632b
begin
	true_eqn =
		x ~ a^2 + log2((a * b^2 * c_plus_alpha ) / sin(2a * b^2)^2) * sqrt(2a * b^2)
	latexify(true_eqn)
end

# ╔═╡ ea8f56c9-11d1-48e3-ae36-678fab1ad4bd
md"## Reconstruct a split-up equation"

# ╔═╡ 8e96dcbb-5320-4222-9b83-aeedce9ed816
md"Now suppose that the equation (due to being more complicated) is presented as split up into parts (using some helper variables `h_1`, `h_2`). These have to be 'substituted back' into each other to get back the full equation. *For using the full equation, it may not be neccessary (or even a good idea) to reassemble the full symbolic expression. Are there better ways?*"

# ╔═╡ 79537fac-84f7-4bc1-8816-7d0ee433408e
begin
	# auxilliary variables
	helper_vars = @variables h_1 h_2

	split_up_eqns = [
		x ~ a^2 + log2(h_1) * sqrt(h_2),
		h_1 ~ h_2 * c_plus_alpha / 2 / sin(h_2)^2,
		h_2 ~ 2a * b^2
	]
	latexify(split_up_eqns)
end

# ╔═╡ 15d8d20e-a209-4941-9d95-dc43a1cde690
md"""
I'm not exactly sure what is the _best_ way of getting back the original equation, but this dirty hack works:
"""

# ╔═╡ 745aa374-ade1-4474-99fd-f51869405a00
"""
**NOTE: Beware, this is just a crappy prototype!**

Substitutes expressions inside `target_expr` as specified by equations `source_eqns`.
For each equation in `source_eqns`, substitutes all occurences of its `lhs` by its `rhs` until none of the `lhs` are left.

**If there are cyclic dependencies among `source_eqns`, the output is undefined!**
"""
function substitute_all(target_expr, source_eqns,)
    all_eqns = [target_expr, source_eqns...]

	# ------ some consistency checks ----
	# left hand sides of `source_eqns` are unique
    eliminated_expressions = map(x->x.lhs, source_eqns)
    length(Set(eliminated_expressions)) == length(source_eqns) || error(
		"Redundant equations! Eliminated expressions must be unique. Nonunique:
		$eliminated_expressions")

	# left hand sides of `source_eqns` have no variables in common
	intersecting_lhs_vars = intersect(
		Symbolics.get_variables.(eliminated_expressions)...)
	isempty(intersecting_lhs_vars) || error(
		"Output depends on order of source equations due to lhs dependencies: lhs's
		$eliminated_expressions with intersecting variables $intersecting_lhs_vars")

	# ------- actual 'work' -------------
    for i in 2:length(all_eqns)
        all_eqns = substitute.(all_eqns, (Dict(all_eqns[i].lhs => all_eqns[i].rhs),))
    end

    return all_eqns[1]
end

# ╔═╡ 92801458-b302-4e89-8434-4e6bb4c0135e
begin # a basic test:
	recovered_eqn = substitute_all(split_up_eqns[1], split_up_eqns[2:end],)
	recovered_eqn == true_eqn || error("Failure: equations don't match!")
end

# ╔═╡ f2357d35-0c51-4b04-a491-3893ac35f5a4
md"""
## Plug in numbers into the equation

this is the easiest way to plug in numbers into the equation that I have seen so far.
"""

# ╔═╡ fcdd4f4a-59d2-4027-8c30-f56d6782b7ef
begin
	values_to_plug_in = Dict([x=>1e9,
			a=>2.0, b=>3.0, c_plus_alpha=>4.0])
	example_equation_result = substitute(recovered_eqn.rhs, values_to_plug_in)
end

# ╔═╡ 8fdb9251-1865-4e5e-9fac-b888c340e066
md"""
- **Note:** `x=>1e9` doesn't do anything! `substitute` does not (nor should it) check if all substitutions actually take place! This is a point where one might want something to help catch errors.

- **Note:**  `c_plus_alpha` is an alias for var"{c^{+}_{\alpha}}", which is how hat symbol thinks it should LaTeX-render and transform into `Symbol`. At this point, using the alias works.
"""

# ╔═╡ f91d951f-2704-4244-a102-c1fc419337cb
md"""
## Compile equation into a Julia function

Suppose you want to plug in lots of different numbers for your variables and look at the restults, maybe 'plot the equation', which may require to broadcast over arrays, etc. My fear is that the `substitute` way would become too slow.

One way to do it is to create a keyword-argument Julia function from the expression.
(This might not be a *good* way to do it. Maybe there are better ways to achieve such tasks?)

`Symbolics.jl` provides `build_function`, which can generate a function that takes its arguments (numbers we want to plug in) as an array.

To get keyword-arguments, one can use `SymbolicUtils.jl`'s `Func` datatype, as well as `toexpr`. However, it seems that these keyword-arguments *must* be provided with default values (I guess that's not so bad), and *must* be named according to `Symbol(variable)`, which interferes with my idea to have weird variable names such as
```julia
@variables var"stupidstuff"
abbreviation = var"stupidstuff" # continue to use abbreviation
```

*Do you have an idea how to deal with that? Of course, using meta-programming, one could build the function manually. However, that's quite a bit of work and hard to get completely right.*
"""

# ╔═╡ 85718649-6050-4cd6-9c6d-83023a6e6030
md"""
### Compile into function that takes arguments as a `Vector`
"""

# ╔═╡ 76c64481-d2af-40b8-bdda-ea2ebc05e8c8
begin
	func_expr = build_function(recovered_eqn, proper_vars)
	Base.remove_linenums!(func_expr) # why does this not remove the linenums?
	func_expr
end

# ╔═╡ 756f6398-2357-4602-8269-d27d60e32b2b
md"""
### Compile into function that takes keyword arguments
"""

# ╔═╡ 27adf4ac-d11a-4e57-b8d9-d69a57f9e045
md"for some reason, `SymbolicUtils.Code` has its own type for saying `x => 1.0`, called `Assignment`. Lets transform the values-to-plug-in above into defaults for the keyword-argument function.
"

# ╔═╡ 105d0d41-6680-48ec-a16c-2c719ab06234
begin
	default_values = Dict([x=>1e9,
				a=>2.0, b=>3.0, c_plus_alpha=>4.0])

	kwargs_spec = SymbolicUtils.Code.Assignment.(
		keys(default_values), values(default_values)
	)
end

# ╔═╡ 26454574-c62a-4851-9211-c56e4da9ae8c
md"Build the keyword-argument function:"

# ╔═╡ eea19aaa-a295-4ae5-88fd-e2730710fbb9
begin
	func = Symbolics.Func(
		[], # positional arguments
		kwargs_spec,
		recovered_eqn.rhs
	)
	kwarg_expr = SymbolicUtils.Code.toexpr(func)
end

# ╔═╡ 30ffb100-2331-494e-9433-09bf96c701fa
md"""
**Note:** we get a keyword-argument named `var"{c^{+}_{\\alpha}}"`, which I don't know how to customize.
"""

# ╔═╡ 5a4b516e-c7d4-4bad-8793-d4c63db23f9b
begin # use kwarg function
	kwargfunction = eval(kwarg_expr)

	kwargfunction() == example_equation_result || error("Mismatched results")

	kwargfunction(a = 42.), kwargfunction(a = 42., var"{c^{+}_{\alpha}}"=1.2e5)
end

# ╔═╡ cd16d92b-c605-4b37-ad42-1125c931d2c6
md"Maybe plot something?"

# ╔═╡ a2464c22-3904-44ab-a00b-59570cdf7b7d
begin
	a_range = range(0.1,stop=10.,length=1000)

	# note: keyword arguments cannot (I think) be broadcast with `func.()` notation
	eqn_values = broadcast(a_range) do aval
		kwargfunction(a=aval)
	end
end

# ╔═╡ a5c6de65-baeb-49eb-a28a-ac1109806794
# begin # compiling Plots.jl can take a few minutes. Uncomment if you dare!
# 	using Plots
# 	plot(a_range, eqn_values, label=recovered_eqn.lhs)
# end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

[compat]
LaTeXStrings = "~1.3.0"
Latexify = "~0.15.9"
PlutoUI = "~0.7.29"
SymbolicUtils = "~0.19.1"
Symbolics = "~4.2.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgCheck]]
git-tree-sha1 = "dedbbb2ddb876f899585c4ec4433265e3017215a"
uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
version = "2.1.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "1ee88c4c76caa995a885dc2f22a5d548dfbbc0ba"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.2.2"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.AutoHashEquals]]
git-tree-sha1 = "45bb6705d93be619b81451bb2006b7ee5d4e4453"
uuid = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
version = "0.2.0"

[[deps.BangBang]]
deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
git-tree-sha1 = "95831c49cf801756a922e50641361e3b4476a782"
uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
version = "0.3.33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Baselet]]
git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
version = "0.1.1"

[[deps.Bijections]]
git-tree-sha1 = "705e7822597b432ebe152baa844b49f8026df090"
uuid = "e2ed5e7c-b2de-5872-ae92-c73ca462fb04"
version = "0.1.3"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "926870acb6cbcf029396f2f2de030282b6bc1941"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.4"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[deps.CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CompositeTypes]]
git-tree-sha1 = "d5b014b216dc891e81fea299638e4c10c657b582"
uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657"
version = "0.1.2"

[[deps.CompositionsBase]]
git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
version = "0.1.1"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DefineSingletons]]
git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
version = "0.1.2"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.DiffRules]]
deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "9bc5dac3c8b6706b58ad5ce24cffd9861f07c94f"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.9.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "6a8dc9f82e5ce28279b6e3e2cea9421154f5bd0d"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.37"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.DomainSets]]
deps = ["CompositeTypes", "IntervalSets", "LinearAlgebra", "StaticArrays", "Statistics"]
git-tree-sha1 = "5f5f0b750ac576bcf2ab1d7782959894b304923e"
uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
version = "0.5.9"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DynamicPolynomials]]
deps = ["DataStructures", "Future", "LinearAlgebra", "MultivariatePolynomials", "MutableArithmetics", "Pkg", "Reexport", "Test"]
git-tree-sha1 = "585de0d658506cf0fe5808026edff662bef5bf03"
uuid = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
version = "0.4.1"

[[deps.EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "3fe985505b4b667e1ae303c9ca64d181f09d5c05"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.3"

[[deps.ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.InitialValues]]
git-tree-sha1 = "40c555f961d7ccf86d8ccd150b9eef379cbfa0a3"
uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
version = "0.3.0"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.LabelledArrays]]
deps = ["ArrayInterface", "ChainRulesCore", "LinearAlgebra", "MacroTools", "StaticArrays"]
git-tree-sha1 = "41158dee1d434944570b02547d404e075da15690"
uuid = "2ee39098-c373-598a-b85f-a56591580800"
version = "1.7.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a8f4f279b6fa3c3c4f1adadd78a621b13a506bce"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.9"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Metatheory]]
deps = ["AutoHashEquals", "DataStructures", "Dates", "DocStringExtensions", "Parameters", "Reexport", "TermInterface", "ThreadsX", "TimerOutputs"]
git-tree-sha1 = "0886d229caaa09e9f56bcf1991470bd49758a69f"
uuid = "e9d8d322-4543-424a-9be4-0cc815abe26c"
version = "1.3.3"

[[deps.MicroCollections]]
deps = ["BangBang", "Setfield"]
git-tree-sha1 = "4f65bdbbe93475f6ff9ea6969b21532f88d359be"
uuid = "128add7d-3638-4c79-886c-908ea0c25c34"
version = "0.1.1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.MultivariatePolynomials]]
deps = ["DataStructures", "LinearAlgebra", "MutableArithmetics"]
git-tree-sha1 = "fa6ce8c91445e7cd54de662064090b14b1089a6d"
uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
version = "0.4.2"

[[deps.MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "73deac2cbae0820f43971fad6c08f6c4f2784ff2"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.2"

[[deps.NaNMath]]
git-tree-sha1 = "f755f36b19a5116bb580de457cda0c140153f283"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.6"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.5"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "7711172ace7c40dc8449b7aed9d2d6f1cf56a5bd"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.29"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "6b96eb51a22af7e927d9618eaaf135a3520f8e2f"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.24.0"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Referenceables]]
deps = ["Adapt"]
git-tree-sha1 = "e681d3bfa49cd46c3c161505caddf20f0e62aaa9"
uuid = "42d2dcc6-99eb-4e98-b66c-637b7d73030e"
version = "0.1.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.2.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.RuntimeGeneratedFunctions]]
deps = ["ExprTools", "SHA", "Serialization"]
git-tree-sha1 = "cdc1e4278e91a6ad530770ebb327f9ed83cf10c4"
uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
version = "0.5.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.SciMLBase]]
deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
git-tree-sha1 = "68c2ee0351edd3c39f67a017d2d1d12f939126dc"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "1.25.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "0afd9e6c623e379f593da01f20590bacc26d1d14"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.8.1"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.0.0"

[[deps.SplittablesBase]]
deps = ["Setfield", "Test"]
git-tree-sha1 = "39c9f91521de844bad65049efd4f9223e7ed43f9"
uuid = "171d559e-b47b-412a-8079-5efa626c420e"
version = "0.1.14"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "7f5a513baec6f122401abfc8e9c074fdac54f6c1"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.4.1"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "88a559da57529581472320892576a486fa2377b9"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.3.1"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "d88665adc9bcf45903013af0982e2fd05ae3d0a6"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "51383f2d367eb3b444c961d485c565e4c0cf4ba0"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.14"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "bedb3e17cc1d94ce0e6e66d3afa47157978ba404"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.14"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.SymbolicUtils]]
deps = ["AbstractTrees", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "Metatheory", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TermInterface", "TimerOutputs"]
git-tree-sha1 = "d8a229358904e51e90a2aa4b48d7ce1f461c7dfe"
uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
version = "0.19.1"

[[deps.Symbolics]]
deps = ["ArrayInterface", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "Metatheory", "NaNMath", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "TermInterface", "TreeViews"]
git-tree-sha1 = "3f354e1fac72482a0ba2ad3d7f9a08e8acc367a5"
uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
version = "4.2.3"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.TermInterface]]
git-tree-sha1 = "7aa601f12708243987b88d1b453541a75e3d8c7a"
uuid = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
version = "0.2.3"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.ThreadsX]]
deps = ["ArgCheck", "BangBang", "ConstructionBase", "InitialValues", "MicroCollections", "Referenceables", "Setfield", "SplittablesBase", "Transducers"]
git-tree-sha1 = "6dad289fe5fc1d8e907fa855135f85fb03c8fa7a"
uuid = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
version = "0.1.9"

[[deps.TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "a5aed757f65c8a1c64503bc4035f704d24c749bf"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.14"

[[deps.Transducers]]
deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
git-tree-sha1 = "3f0945b47207a41946baee6d1385e4ca738c25f7"
uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
version = "0.4.68"

[[deps.TreeViews]]
deps = ["Test"]
git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
version = "0.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─b9b5c2e0-537d-4a72-adb9-8474dc8acf98
# ╟─42f6afaa-cf0d-46dd-b1d9-ebaeafc9e5db
# ╠═4c4bfc78-7224-11ec-211d-f1e60824d35c
# ╟─b580a960-7780-417e-ae0f-3ac59f4a4a30
# ╠═0c2701f4-8e93-4fbc-b18c-52e9d1c213f8
# ╟─6e521fb4-190c-4b68-890b-d39b22537471
# ╟─bf6bac04-1957-4aa9-976e-efd76b41632b
# ╟─ea8f56c9-11d1-48e3-ae36-678fab1ad4bd
# ╟─8e96dcbb-5320-4222-9b83-aeedce9ed816
# ╟─79537fac-84f7-4bc1-8816-7d0ee433408e
# ╟─15d8d20e-a209-4941-9d95-dc43a1cde690
# ╟─745aa374-ade1-4474-99fd-f51869405a00
# ╠═92801458-b302-4e89-8434-4e6bb4c0135e
# ╟─f2357d35-0c51-4b04-a491-3893ac35f5a4
# ╠═fcdd4f4a-59d2-4027-8c30-f56d6782b7ef
# ╟─8fdb9251-1865-4e5e-9fac-b888c340e066
# ╟─f91d951f-2704-4244-a102-c1fc419337cb
# ╟─85718649-6050-4cd6-9c6d-83023a6e6030
# ╠═76c64481-d2af-40b8-bdda-ea2ebc05e8c8
# ╟─756f6398-2357-4602-8269-d27d60e32b2b
# ╟─27adf4ac-d11a-4e57-b8d9-d69a57f9e045
# ╠═105d0d41-6680-48ec-a16c-2c719ab06234
# ╟─26454574-c62a-4851-9211-c56e4da9ae8c
# ╠═eea19aaa-a295-4ae5-88fd-e2730710fbb9
# ╟─30ffb100-2331-494e-9433-09bf96c701fa
# ╠═5a4b516e-c7d4-4bad-8793-d4c63db23f9b
# ╟─cd16d92b-c605-4b37-ad42-1125c931d2c6
# ╠═a2464c22-3904-44ab-a00b-59570cdf7b7d
# ╠═a5c6de65-baeb-49eb-a28a-ac1109806794
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
struct NotInField <: Exception end
Base.showerror(io::IO, e::NotInField) = print(io, "𝑛 is not within secp256k1 field")

struct NotOnCurve <: Exception end
Base.showerror(io::IO, e::NotOnCurve) = print(io, "Point is not on curve")

abstract type SignatureException <: Exception end

struct PrefixError <: SignatureException end
Base.showerror(io::IO, e::PrefixError) = print(io, "Incorrect signature prefix")

struct LengthError <: SignatureException end
Base.showerror(io::IO, e::LengthError) = print(io, "Encoded length does not match available data")


text:
using LimitOrderBook
using Test

@testset "LimitOrderBook.jl" begin
    # Write your tests here.
end


text:
function calc_Vandermonde_int_int(x,ni,nj)

    V=zeros(BigFloat, ni, nj)
    for i=1:ni
        for j=1:nj
            V[i,j] = (x[i]^BigInt(j+1))/BigInt(j+1)/BigInt(j)
        end
    end

    return V

end


text:
function burkey_cantrell()
    e_diabat = zeros(Float64, Ms)
    h0 = zeros(Float64, Ms+1, Ms+1)
    vm = zeros(Float64, Ms+1, Ms+1)
    gauss = zeros(Float64, Int(Ms/2), Int(Ms/2))
    for j in 1:(Int(Ms/2)-1)
        gauss[j, j+1] = Float64(j)/sqrt(4.0*Float64(j)^2 - 1.0)
        gauss[j+1, j] = gauss[j, j+1]
    end

    eigvals_gauss = eigvals(gauss)
    eigvecs_gauss = eigvecs(gauss)

    for j in 2:(Int(Ms/2) + 1)
        h0[j, j] = delta_E/4.0*eigvals_gauss[j-1] - delta_E/4.0
        h0[j+Int(Ms/2), j+Int(Ms/2)] = delta_E/4.0*eigvals_gauss[j-1] + delta_E/4.0
        e_diabat[j-1] = h0[j, j]
        e_diabat[j + Int(Ms/2) - 1] = h0[j+Int(Ms/2), j+Int(Ms/2)]
    end

    vm[1, 2:Int(Ms/2) + 1] = @. sqrt(delta_E/4.0 * (2.0*eigvecs_gauss[1, :]^2))
    vm[2:Int(Ms/2) + 1, 1] = @. sqrt(delta_E/4.0 * (2.0*eigvecs_gauss[1, :]^2))
    vm[1, (2 + Int(Ms/2)):(Ms + 1)] = @. sqrt(delta_E/4.0 * (2.0*eigvecs_gauss[1, :]^2))
    vm[(2 + Int(Ms/2)):(Ms + 1), 1] = @. sqrt(delta_E/4.0 * (2.0*eigvecs_gauss[1, :]^2))

    return h0, e_diabat, vm
end


text:
VERSION >= v"0.4" #&& __precompile__()

module GMT

using Compat
import Compat.String

export
	GMT_Call_Module,
	GMT_Create_Args,
	GMT_Create_Session,
	GMT_Create_Options,
	GMT_Create_Cmd,
	GMT_Create_Data,
	GMT_Destroy_Options,
	GMT_Destroy_Args,
	GMT_Destroy_Session,
	GMT_Encode_ID,
	GMT_Get_ID,
	GMT_Get_Default,
	GMT_Option,
	GMT_Read_Data,
	GMT_Register_IO,
	GMT_Get_Data,
	GMT_Retrieve_Data,
	GMT_Destroy_Data,
	GMT_Message,
	GMT_Call_Module,
	GMT_Write_Data,
	GMT_Init_VirtualFile, GMT_Open_VirtualFile, GMT_Close_VirtualFile, GMT_Read_VirtualFile,
	GMT_Get_Matrix, GMT_Put_Matrix, GMT_Get_Vector, GMT_Put_Vector,
	GMT_IS_DATASET, GMT_IS_TEXTSET, GMT_IS_GRID, GMT_IS_LINE,
	GMT_IS_CPT, GMT_IS_IMAGE, GMT_IS_VECTOR, GMT_IS_MATRIX,
	GMT_IS_COORD, GMT_IS_POINT,	GMT_IS_MATRIX, GMT_IS_SURFACE,
	GMT_DATASET, GMT_IS_PLP,
	GMT_GRID, GMT_MATRIX, GMT_PALETTE, GMT_PS,
	GMT_UNIVECTOR,
	GMT_IN, GMT_OUT, GMT_OK,
	GMT_IS_FILE, GMT_IS_STREAM,	GMT_IS_FDESC,
	GMT_IS_DUPLICATE, GMT_IS_REFERENCE,
	GMT_VIA_NONE, GMT_VIA_VECTOR, GMT_VIA_MATRIX, GMT_VIA_OUTPUT,
	GMT_IS_DUPLICATE_VIA_VECTOR, GMT_IS_REFERENCE_VIA_VECTOR,
	GMT_IS_DUPLICATE_VIA_MATRIX, GMT_IS_REFERENCE_VIA_MATRIX,
	GMT_MODULE_EXIST, GMT_MODULE_PURPOSE, GMT_MODULE_OPT, GMT_MODULE_CMD,
	GMT_GRID_DATA_ONLY, GMT_GRID_HEADER_ONLY, GMT_GRID_ALL,
	GMT_GRID_ALL, GMT_READ_NORMAL, GMT_WRITE_SET, GMT_Report,
	array_container,
	GMT_Encode_Options, GMT_Expand_Option, GMT_grid_flip_vertical,
	gmt_core_module_info,
	gmt, grid_type,
	GMT_RESOURCE, GMTAPI_CTRL, GMTAPI_DATA_OBJECT,
	NULL, CHAR

include("libgmt_h.jl")
include("libgmt.jl")
include("gmt_main.jl")

end # module


text:
"""
    BetaFourParameters
* Parametrization a, c, \\alpha, \\beta
* Score
* Fisher Information
* `time_varying_params` map.
* Default link

Right now when estimating a BetaFourParameters model is recomended to provide fixed parameters a and c, dont
with the code:

```julia
gas_beta_4 = Model([1, 2, 11, 12], [1, 2, 11, 12], BetaFourParameters, 0.0; time_varying_params=[3])
gas_beta_4.ω[1] = minimum(y) - 0.1*std(y) # parameter a
gas_beta_4.ω[2] = maximum(y) + 0.1*std(y) # parameter c
```

This code is a simple and not very accurate heuristic on how to estimate the a and c parameters. If you have estimated
using maximum likelihood it will be probably better

The fact occurs because if in the optimization the gradient makes `c < maximum(y)` or `a > minimum(y)` then it
is impossible that y comes from the BetaFourParameters distribution leading to DomainErrors
"""
BetaFourParameters

function score!(score_til::Matrix{T}, y::T, ::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = -(param[t, 3] - 1)/(y - param[t, 1]) + (param[t, 4] + param[t, 3] - 1)/(param[t, 2] - param[t, 1])
    score_til[t, 2] = (param[t, 4] - 1)/(param[t, 2] - y) - (param[t, 4] + param[t, 3] - 1)/(param[t, 2] - param[t, 1])
    score_til[t, 3] = log(y - param[t, 1]) - log(param[t, 2] - param[t, 1]) +
                        digamma(param[t, 3] + param[t, 4]) - digamma(param[t, 3])
    score_til[t, 4] = log(param[t, 2] - y) - log(param[t, 2] - param[t, 1]) +
                        digamma(param[t, 3] + param[t, 4]) - digamma(param[t, 4])
    return
end

function log_likelihood(::Type{BetaFourParameters}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = 0.0
    for t in 1:n
        loglik += (param[t, 3] - 1)*log(y[t] - param[t, 1]) + (param[t, 4] - 1) * log(param[t, 2] - y[t]) -
                  (param[t, 3] + param[t, 4] - 1) * log(param[t, 2] - param[t, 1]) - logbeta(param[t, 3], param[t, 4])
    end
    return -loglik
end

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

# utils
function update_dist(::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T
    small_threshold!(param[:, 3:4], SMALL_NUM, t)
    beta = Beta(param[t, 3], param[t, 4])
    return LocationScale(param[t, 1], param[t, 2] - param[t, 1], beta)
end

function params_sdm(d::BetaFourParameters)
    return (d.μ, d.σ + d.μ, Distributions.params(d.ρ)...)
end

function num_params(::Type{BetaFourParameters})
    return 4
end

# d = Beta(1, 1)
# d2 = LocationScale(-1, 3, d)
# maximum(d2)
# minimum(d2)


text:
# test against Farmer Toda MATLAB values

hex2num_me(s) = reinterpret(Float64, parse(UInt64, s, base = 16))

DFT99 = dropdims(hex2num_me.(["c0340cd923b15b89"   "c02e1345b58a094e"   "c0240cd923b15b89"   "c0140cd923b15b88"   "0000000000000000"   "40140cd923b15b88"   "40240cd923b15b8a"   "402e1345b58a094e"   "40340cd923b15b89"]); dims=1)
PFT99 = hex2num_me.(
      [ "3feebdff4a767113"   "3fa3c61f93aef8f1"   "3f464066ed540f5d"   "3edd4507c34018f2"   "3dded173035783dd"   "3c49f92e9f402b9c"   "3a21859dbed5564e"   "3762ec61a6257787"   "34105bcf2cf90500";
        "3d99649e3e85a1e9"   "3fef2e561795e507"   "3f96a88e64e8b5cd"   "3f64ba2e346747b5"   "3f4e8d7202f75e32"   "3ecfa9d74ae8221e"   "3de7124bc92f5c91"   "3c97a393256534a1"   "3ae10768d4a1c15f";
        "3c1af68f17efc8f0"   "3dc8b97cb927921b"   "3fef97050c545fe9"   "3f814097bd0adae2"   "3f5c24ab3f85ffdf"   "3f64e95cfe5cc002"   "3f1f902982d4999c"   "3e883178a38482b3"   "3da2d5af6cc1ee2f";
        "3c65cee5b29700db"   "3d2b083ffb64d194"   "3e3ef0b49f17ce6a"   "3fefe89723f18a15"   "3f3e2eb32f62d91c"   "3f12d4fab284def1"   "3f41eede5c549b17"   "3f529bdeda0b2764"   "3f450acbc309a532";
        "3f5444b0bd52b3ea"   "3ee6ebb33b5859ca"   "3e9906756359ae1f"   "3edcc643e8f3b048"   "3fefeb6f21308b9c"   "3edcc643e93f9d2c"   "3e99067563ddbec7"   "3ee6ebb33c0dc9f9"   "3f5444b0be28a035";
        "3f450acbd076ef71"   "3f529bdec9630975"   "3f41eede3d641bd9"   "3f12d4fa8fa25180"   "3f3e2eb308db0eb2"   "3fefe897240426f9"   "3e3ef0b4dc7bb8a6"   "3d2b0840799bcb9d"   "3c65cee6624c25a4";
        "3da2d5af6d4add66"   "3e883178a3f70f30"   "3f1f9029832ab837"   "3f64e95cfe76ec17"   "3f5c24ab3f8e8e0a"   "3f814097bd08e239"   "3fef97050c5446a5"   "3dc8b97cb9375469"   "3c1af68f182082f7";
        "3ae10768d5312c4c"   "3c97a39325f2c4a0"   "3de7124bc98ad89f"   "3ecfa9d74b3299b9"   "3f4e8d72031a74db"   "3f64ba2e346ea424"   "3f96a88e64e812c8"   "3fef2e561795d9ea"   "3d99649e3e90b030";
        "34105bd6cc1238fd"   "3762ec6842abd059"   "3a2185a21e0a922f"   "3c49f932f0f81e51"   "3dded1761626dc19"   "3edd450938876835"   "3f4640674ba35ba6"   "3fa3c61f936a7347"   "3feebdff4a626af6"]
      )

# setup
ρ = 0.99
n = 9
L = 2
σ = 1.

# set up state space
σlr = sqrt(1.0 / (1.0-ρ^2))
nσ = sqrt(n-1)
y = range(-σlr*nσ, stop=σlr*nσ, length=n)
@test y ≈ DFT99

# output
P, JN, Λ, numMoments, approxErr = MarkovTransitionMatrices.discreteNormalApprox(y, y, (x::Real, st::Real) -> x - ρ*st, 2)
@test maximum(abs.(P .- PFT99)) < 1e-7
@test maximum(abs.(0.5 * (P .- PFT99) ./ (P .+ PFT99))) < 1e-5

#