text:
export SellmeierModel, sellmeier_model

"""
    SellmeierModel{N}

Sellmeier model representing the Sellmeier equation ε(λ) = 1 + ∑ᵢ [Bᵢ λ² / (λ² - Cᵢ)], where
the sum is for `1 ≤ i ≤ N`.

# Fields
- `str::SVector{N,Float64}`: `str[i]` is the strength of the `i`th term in the Sellmeier equation; corresponds to Bᵢ.
- `λres::SVector{N,Float64}`: `λres[i]` is the resonance wavelength of the `i`th term in the Sellmeier equation; corresponds to √Cᵢ.
"""
struct SellmeierModel{N}
    str::SFloat{N}  # strength of terms; typically written str in equation
    λres::SFloat{N}  # resonance wavelength of terms; typically written √λres in equation

    SellmeierModel{N}(str::SReal{N}, λres::SReal{N}) where {N} = (p = sortperm(λres); new(str[p], λres[p]))
end

SellmeierModel(str::AbsVecReal, λres::AbsVecReal) = (N = length(str); SellmeierModel{N}(SFloat{N}(str), SFloat{N}(λres)))

"""
    Base.length(sm::SellmeierModel)

Return the number of terms in the sum of the Sellmeier equation represented by the given
`SellmeierModel`.
"""
Base.length(sm::SellmeierModel{N}) where {N} = N

"""
    (::SellmeierModel)(λ::Real)

Evaluate the Sellmeier equation represented by the given `SellmeierModel` at a wavelngth `λ`
in units of meters.
"""
(sm::SellmeierModel)(λ::Real) = sellmeier_model(λ, sm.str, sm.λres)

# Returns the relative permittivity ε = n² at a wavelength λ following the Sellmeier model.
# The fitting equation is ε (= n²) rather than n, because ε is sometimes evaluated negative
# during the optimization procedure; n = √ε is not defined for real negative ε.
sellmeier_model(λ::Real, str::AbsVecReal, λres::AbsVecReal) =
    1.0 + mapreduce((strᵢ,λresᵢ) -> strᵢ / (1.0 - (λresᵢ/λ)^2), +, str, λres)

ext: jl


text:
using Test
using RootedTrees

@testset "Create tree manually" begin
    #==== Test topology:

                c1
               /
        *p ―⟶ x ―⟶ c2
        |      \
        s      c3

    * `p` is the root.

    Or as a Newick string:
    ((c1,c2,c3)x,s)p;
    ====================#
    tree = createtree(6)

    setlabel!.(getnodes(tree), ("p", "x", "c1", "c2", "c3", "s"))

    p, x, c1, c2, c3, s = getnodes(tree)

    link!(p => x)
    link!(p => s)
    link!(x => c1)
    link!(x => c2)
    link!(x => c3)

    tree.root = p

    @test newick(tree) == "((c1,c2,c3)x,s)p;"

end

@testset "Read Newick" begin
    # Example strings from Wikipedia
    test_strings = [
        "(,,(,));",                              ## no nodes are named
        "(A,B,(C,D));",                          ## leaf nodes are named
        "(A,B,(C,D)E)F;",                        ## all nodes are named
        "(:0.1,:0.2,(:0.3,:0.4):0.5);",         ## all but root node have a distance to parent
        "(:0.1,:0.2,(:0.3,:0.4):0.5):0.0;" ,     ## all have a distance to parent
        "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);",      ## distances and leaf names (popular)
        "(A:0.1,
        B:0.2 , (C:0.3,D:0.4) : 0.5 ) ;",      ## same as above, with whitespace
        "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;",    ## distances and all names
        "((B:0.2,(C:0.3,D:0.4)E:0.5)F:0.1)A;"   ## a tree rooted on a leaf node (rare)
        ]

    for nwkstr in test_strings
        tree = parsenewick(nwkstr)
        @test newick(tree) == filter(! isspace, nwkstr)
    end

    test_string =
    "((\"B\":0.2,(\"C\":0.3,\"D\":0.4)\"E\":0.5)\"F\":0.1)A;"   ## with quote marks

    tree = parsenewick(test_string)
    @test newick(tree) == "((B:0.2,(C:0.3,D:0.4)E:0.5)F:0.1)A;"

end

@testset "Traversals" begin
    tree = createtree(6)

    setlabel!.(getnodes(tree), ("p", "x", "c1", "c2", "c3", "s"))

    p, x, c1, c2, c3, s = getnodes(tree)

    tree.root = p

    link!(p => x)
    link!(p => s)
    link!(x => c1)
    link!(x => c2)
    link!(x => c3)

    @test preorder(tree) == preorder(p) == [p, x, c1, c2, c3, s]
    @test preorder(x) == [x, c1, c2, c3]
    @test postorder(tree) == postorder(p) == [c1, c2, c3, x, s, p]
    @test postorder(x) == [c1, c2, c3, x]
end

@testset "Reroot" begin
    tree = createtree(6)

    setlabel!.(getnodes(tree), ("p", "x", "c1", "c2", "c3", "s"))

    p, x, c1, c2, c3, s = getnodes(tree)

    link!(p => x)
    link!(p => s)
    link!(x => c1)
    link!(x => c2)
    link!(x => c3)

    tree.root = p

    reroot!(tree, x)
    @test isnothing(parent(x))
    @test parent(p) ≡ x
    @test p ∈ children(x)
    @test c1 ∈ children(x)
    @test c2 ∈ children(x)
    @test s ∈ children(p)
end

function findnode(tree; label, which = :first)
    if which == :first
        ids = findfirst(x -> getlabel(x) == label, getnodes(tree))
    elseif which == :all
        ids = findall(x -> getlabel(x) == label, getnodes(tree))
    else
        throw(ArgumentError("invalid `which` argument"))
    end

    return isnothing(ids) ? nothing : getnode(tree, ids)
end

@testset "Reroot" begin
    tree = parsenewick("((A:4.0,B:4.0)F:2.0,(C:4.0,(D:8.0,E:8.0)G:4.0)H:2.0)I:1.0;")
    #=
    RNode{Dict}: #1 — "I"
    RNode{Dict}: #2 — "F"
    RNode{Dict}: #3 — "A"
    RNode{Dict}: #4 — "B"
    RNode{Dict}: #5 — "H"
    RNode{Dict}: #6 — "C"
    RNode{Dict}: #7 — "G"
    RNode{Dict}: #8 — "D"
    RNode{Dict}: #9 — "E"
    =#
    c = tree.nodes[6]
    d = tree.nodes[8]
    e = tree.nodes[9]
    f = tree.nodes[2]
    g = tree.nodes[7]
    h = tree.nodes[5]
    i = tree.nodes[1]

    reroot!(tree, g)

    @test Set(collect(children(g))) == Set([d, e, h])
    @test Set(collect(children(h))) == Set([i, c])
    @test Set(collect(children(i))) == Set([f])
    @test isnothing(parent(g))
    @test parent(h) ≡ g
    @test parent(i) ≡ h

    @test brlength(g) == 1.0
    @test brlength(h) == 4.0
    @test brlength(i) == 2.0
    @test brlength(f) == 2.0
end

@testset "Ancestors" begin
    #=
    Make sure that node IDs are in preorder
    =#
    tree16 = symmetric_tree(Nothing, Nothing, 16)

    @testset "path_to_root" begin
        c = getnode(tree16, 12)
        E = path_to_root(c)

        @test length(E) == 4
        @test all(getid.(E) .== [11, 10, 2, 1])
    end

    @testset "findmrca" begin
        #TODO: Add test for auto method selection
        #TODO: Fix method 2!
        @test 1 == findmrca(tree16, [5, 31], method = 1)
        @test 1 == findmrca(tree16, [5, 31], method = 2)
        @test 2 == findmrca(tree16, [5, 16], method = 1)
        @test_broken 2 == findmrca(tree16, [5, 16], method = 2)
        @test 25 == findmrca(tree16, [27, 31], method = 1)
        @test_broken 25 == findmrca(tree16, [27, 31], method = 2)
        @test 10 == findmrca(tree16, [12, 16], method = 1)
        @test_broken 10 == findmrca(tree16, [12, 16], method = 2)
        @test 19 == findmrca(tree16, [20, 21], method = 1)
        @test_broken 19 == findmrca(tree16, [20, 21], method = 2)

        # More than 2 targets
        @test 17 == findmrca(tree16, [20, 21, 31], method = 1)
        @test_broken 17 == findmrca(tree16, [20, 21, 31], method = 2)

        # Inner nodes
        @test 17 == findmrca(tree16, [19, 25, 29], method = 1)
        @test_broken 17 == findmrca(tree16, [19, 25, 29], method = 2)

    #     #=
    #     Need to decide on the behaviour of findmrca when one of the nodes is an ancestor
    #     of another (including the root!).
    #     Also: can a node be its own ancestor ?
    #     =#
    #     # Any node with the root
    #     @test_skip 0 == findmrca(tree16, [1, 25], method = 1)
    #     @test_skip 0 == findmrca(tree16, [1, 25], method = 2)
    end

    #! Need to define findmrca behaviour for all combinations of nodes
    # @testset "mrca_matrix" begin
    #     A = mrca_matrix(tree16)

    #     for i in 2:31
    #         for j in 1:(i-1)
    #             @test_skip A[j, i] == findmrca(tree16, [i, j])
    #         end
    #     end
    # end
end

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
module TestMOIwrapper

using CPLEX
using MathOptInterface
using Test

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_runtests()
    model = MOI.Bridges.full_bridge_optimizer(CPLEX.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    # Turn off presolve reductions so CPLEX will generate infeasibility
    # certificates.
    MOI.set(model, MOI.RawOptimizerAttribute("CPX_PARAM_REDUCE"), 0)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(atol = 1e-3, rtol = 1e-3),
        exclude = String[
            # CPLEX doesn't support nonconvex QCPs
            "test_quadratic_nonconvex_",
            # TODO(odow): fixed in MOI 0.10.1
            "test_solve_SOS2_add_and_delete",
            "test_objective_get_ObjectiveFunction_ScalarAffineFunction",
            "test_conic_SecondOrderCone_negative_initial_bound",
            "test_conic_SecondOrderCone_negative_post_bound",
            "test_conic_SecondOrderCone_negative_post_bound_ii",
            "test_conic_SecondOrderCone_no_initial_bound",
            "test_conic_SecondOrderCone_nonnegative_initial_bound",
            "test_conic_SecondOrderCone_nonnegative_post_bound",
        ],
    )
    return
end

"""
Test setting CPXPARAM_OptimalityTarget because it changes the problem type.
"""
function test_CPXPARAM_OptimalityTarget()
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.Interval(1.0, 4.0))
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        MOI.ScalarQuadraticFunction(
            [MOI.ScalarQuadraticTerm(2.0, x, x)],
            MOI.ScalarAffineTerm{Float64}[],
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(
        model,
        MOI.RawOptimizerAttribute("CPXPARAM_OptimalityTarget"),
        CPX_OPTIMALITYTARGET_OPTIMALGLOBAL,
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 16.0 atol = 1e-6
    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 4.0 atol = 1e-6
    return
end

function test_user_provided_env()
    env = CPLEX.Env()
    model_1 = CPLEX.Optimizer(env)
    @test model_1.env === env
    model_2 = CPLEX.Optimizer(env)
    @test model_2.env === env
    # Check that finalizer doesn't touch env when manually provided.
    finalize(model_1)
    @test env.ptr != C_NULL
    return
end

function test_automatic_env()
    model_1 = CPLEX.Optimizer()
    model_2 = CPLEX.Optimizer()
    @test model_1.env.ptr !== model_2.env.ptr
    return
end

function test_user_provided_env_empty()
    env = CPLEX.Env()
    model = CPLEX.Optimizer(env)
    @test model.env === env
    @test env.ptr != C_NULL
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr != C_NULL
    return
end

function test_automatic_env_empty()
    model = CPLEX.Optimizer()
    env = model.env
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr != C_NULL
    return
end

function test_manual_env()
    env = CPLEX.Env()
    model = CPLEX.Optimizer(env)
    finalize(env)
    @test env.finalize_called
    finalize(model)
    @test env.ptr == C_NULL
    return
end

function test_ZeroOne_NONE()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 0.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 1.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -CPLEX.CPX_INFBOUND
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == CPLEX.CPX_INFBOUND
    return
end

function test_ZeroOne_LESS_THAN()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.LessThan(2.0))
    c = MOI.add_constraint(model, x, MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 0.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -CPLEX.CPX_INFBOUND
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    return
end

function test_ZeroOne_GREATER_THAN()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(-2.0))
    c = MOI.add_constraint(model, x, MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 1.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == CPLEX.CPX_INFBOUND
    return
end

function test_ZeroOne_INTERVAL()
    model = CPLEX.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.Interval(-2.0, 2.0))
    c = MOI.add_constraint(model, x, MOI.ZeroOne())
    tmp = Ref{Cdouble}()
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    MOI.delete(model, c)
    CPLEX.CPXgetlb(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == -2.0
    CPLEX.CPXgetub(model.env, model.lp, tmp, 0, 0)
    @test tmp[] == 2.0
    return
end

function test_fake_status()
    model = CPLEX.Optimizer()
    model.ret_optimize = CPLEX.CPXERR_NO_MEMORY
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.MEMORY_LIMIT
    @test MOI.get(model, MOI.RawStatusString()) ==
          "CPLEX Error  1001: Out of memory.\n"
    return
end

end  # module TestMOIwrapper

TestMOIwrapper.runtests()


text:
using Pkg; Pkg.activate("../")
using WebAssemblyText

wat2wasm="/home/andy/programs/wabt/bin/wat2wasm"
wasmvalidate="/home/andy/programs/wabt/bin/wasm-validate"
wasmopt="/home/andy/programs/binaryen/bin/wasm-opt"

stripextension(filename) = filename[1:end-3]

filenames = readdir("./jl")
names = stripextension.(filenames)

for name in names
  printstyled("Compiling jl/$name.jl -> wat/$name.wat -> wasm/$name.wasm\n"; color=:yellow)
  outfile = "wat/$name.wat"

  wat = jl2wat("jl/$name.jl")
  open(outfile, "w") do io
    write(io, wat)
  end
  printstyled("...wat", color=:green)

  run(`$wat2wasm wat/$name.wat -o wasm/$name.wasm`)
  printstyled("...wasm", color=:green)
  run(`$wasmvalidate wasm/$name.wasm`)
  printstyled("...validated", color=:green)
  run(`$wasmopt --enable-multivalue -O4 wasm/$name.wasm -o wasm/$name.wasm`)
  printstyled("...optimized\n", color=:green)
end


text:
abstract type GenerateBibliography <: Builder.DocumentPipeline end
Selectors.order(::Type{GenerateBibliography}) = 5.0
Selectors.matcher(::Type{GenerateBibliography}) = true
Selectors.strict(::Type{GenerateBibliography}) = false
function Selectors.runner(::Type{GenerateBibliography}, doc::Documents.Document)
    Utilities.log(doc, "generating bibliograpies.")
    generate_bibliographies(doc)
end

function generate_bibliographies(doc::Documents.Document)
    for (src, page) in doc.internal.pages
        empty!(page.globals.meta)
        for element in page.elements
            generate_bibliography(element, page.mapping[element], page, doc)
        end
    end
end

generate_bibliography(element, node, page, doc) = nothing
function generate_bibliography(element, bib::BibliographyNode, page, doc)
    Selectors.dispatch(CitationStyle, bib.style, doc, bib)

    list = Markdown.List(1)
    for key in bib.keys
        push!(list.items, bib.database[key][1])
    end

    page.mapping[element] = list
    return nothing
end


text:
# Test supports_in_set functions
@testset "supports_in_set" begin
    # supports_in_set (IntervalSet)
    @testset "IntervalSet" begin
        set = IntervalSet(0, 1)
        @test supports_in_set(0, set)
        @test !supports_in_set(-1, set)
        @test !supports_in_set(2, set)
    end
    # supports_in_bounds (UnivariateDistribution set)
    @testset "Univariate Distribution" begin
        set = DistributionSet(Uniform())
        @test supports_in_set(0, set)
        @test !supports_in_set(-1, set)
        @test !supports_in_set(2, set)
    end
    # supports_in_set (Fallback)
    @testset "Fallback" begin
        @test supports_in_set(0, BadSet())
    end
end

# Test lower bound functions
@testset "Lower Bound" begin
    set1 = IntervalSet(0, 1)
    set2 = DistributionSet(Uniform())
    set3 = BadSet()
    # JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(set1)
        @test has_lower_bound(set2)
        @test !has_lower_bound(set3)
    end
    # JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(set1) == 0.0
        @test lower_bound(set2) == 0.0
        @test_throws ErrorException lower_bound(set3)
    end
    # JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test set_lower_bound(set1, 0.5) == IntervalSet(0.5, 1)
        @test_throws ErrorException set_lower_bound(set2, 2)
        @test_throws ErrorException set_lower_bound(set3, 2)
    end
end

# Test upper bound functions
@testset "Upper Bound" begin
    set1 = IntervalSet(0, 1)
    set2 = DistributionSet(Uniform())
    set3 = BadSet()
    # JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(set1)
        @test has_upper_bound(set2)
        @test !has_upper_bound(set3)
    end
    # JuMP.upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(set1) == 1.0
        @test upper_bound(set2) == 1.0
        @test_throws ErrorException upper_bound(set3)
    end
    # JuMP.set_upper_bound
    @testset "JuMP.set_upper_bound" begin
        @test set_upper_bound(set1, 0.5) == IntervalSet(0, 0.5)
        @test_throws ErrorException set_upper_bound(set2, 2)
        @test_throws ErrorException set_upper_bound(set3, 2)
    end
end

# Test support generation
@testset "generate_support_values" begin
    @testset "IntervalSet" begin
        set = IntervalSet(0., 1.)
        @test generate_support_values(set, num_supports = 10, sig_fig = 3) isa Vector{<:Number}
        @test generate_support_values(set, num_supports = 10, sig_fig = 3)[2] == 0.111
        @test generate_support_values(set, num_supports = 10, sig_fig = 3)[2] != 1/11
        @test length(generate_support_values(set, num_supports = 10, sig_fig = 3)) == 10
    end
    @testset "DistributionSet" begin
        dist1 = Normal(0., 1.)
        dist2 = MvNormal([0.; 0.], [1. 0.; 0. 2.])
        set1 = DistributionSet(dist1)
        set2 = DistributionSet(dist2)
        @test generate_support_values(set1, num_supports = 10) isa Vector{<:Number}
        @test generate_support_values(set2, num_supports = 10) isa Array{<:Number, 2}
        @test length(generate_support_values(set1, num_supports = 10)) == 10
        @test size(generate_support_values(set2, num_supports = 10)) == (2, 10)
    end
    @testset "Fallback" begin
        @test_throws ErrorException generate_support_values(BadSet())
    end
end


text:
ENV["JULIA_PARDISO"] = "/Users/gong/lib/Pardiso/libpardiso600-MACOS-X86-64.dylib"
ENV["PYTHON"] = string(ENV["HOME"], "/anaconda3/bin/python3")
#ENV["PYTHON"] = string(ENV["HOME"], "/.julia/conda/3/bin/python3")
#ENV["CONDA_JL_HOME"] =

ocean_julia = ENV["HOME"] * "/GitHub/ocean_julia"
if (ocean_julia in LOAD_PATH) == false
    push!(LOAD_PATH, ocean_julia);
end

ocean_julia_npp = ENV["HOME"] * "/GitHub/ocean_julia/NPP"
if (ocean_julia_npp in LOAD_PATH) == false
    push!(LOAD_PATH, ocean_julia_npp);
end

if (pwd() in LOAD_PATH) == false
    push!(LOAD_PATH, pwd());
end

using Pkg


text:
"""
    struct Jet{T, N}

Represents the truncated (N-1)-th order Taylor series

    f(a) + (x-a)f'(a) + 1/2(x-a)^2f''(a) + ...

Coefficients are stored in unscaled form.
For a jet `j`, several operations are supported:

1. Indexing `j[i]` returns fᵢ
2. Evaluation j(x) semantically evaluates the truncated taylor series at
   `x`. However, evaluation is restricted to be precisely at `a` - the
   additional information in the taylor series is only available through
   derivatives. Mathematically this corresponds to an infinitessimal ball
   around `a`.
"""
struct Jet{T, N}
    a::T
    f₀::T
    fₙ::NTuple{N, T}
end

function ChainRulesCore.rrule(::typeof(Base.getproperty), j::Jet, s)
    error("Raw getproperty not allowed in AD code")
end

function Base.:+(j1::Jet{T, N}, j2::Jet{T, N}) where {T, N}
    @assert j1.a === j2.a
    Jet{T, N}(j1.a, j1.f₀ + j2.f₀, map(+, j1.fₙ, j2.fₙ))
end

function Base.:+(j::Jet{T, N}, x::T) where {T, N}
    Jet{T, N}(j.a, j.f₀+x, j.fₙ)
end

struct One; end

function Base.:+(j::Jet, x::One)
    j + one(j[0])
end

function ChainRulesCore.rrule(::typeof(+), j::Jet, x::One)
    j + x, Δ->(NoTangent(), One(), ZeroTangent())
end

function Base.zero(j::Jet{T, N}) where {T, N}
    let z = zero(j[0])
        Jet{T, N}(j.a, z,
            ntuple(_->z, N))
    end
end
function ChainRulesCore.rrule(::typeof(Base.zero), j::Jet)
    zero(j), Δ->(NoTangent(), ZeroTangent())
end

function Base.getindex(j::Jet{T, N}, i::Integer) where {T, N}
    (0 <= i <= N) || throw(BoundsError(j, i))
    i == 0 && return j.f₀
    @inbounds j.fₙ[i]
end

function deriv(j::Jet{T, N}) where {T, N}
    Jet{T, N-1}(j.a, j.fₙ[1], Base.tail(j.fₙ))
end

function integrate(j::Jet{T, N}) where {T, N}
    Jet{T, N+1}(j.a, zero(j.f₀), tuple(j.f₀, j.fₙ...))
end

deriv(::NoTangent) = NoTangent()
integrate(::NoTangent) = NoTangent()

antideriv(j, Δ) = integrate(zero(j) + Δ)

function ChainRulesCore.rrule(::typeof(deriv), j::Jet)
    deriv(j), Δ->(NoTangent(), antideriv(j, Δ))
end

function ChainRulesCore.rrule(::typeof(integrate), j::Jet)
    integrate(j), Δ->(NoTangent(), deriv(Δ))
end

function Base.show(io::IO, j::Jet)
    print(io, join(map(enumerate(j.fₙ)) do (n, x)
        n -= 1
        buf = IOBuffer()
        if n != 0
            if n != 1
                print(buf, "1/",factorial(n))
            end
            print(buf, "(x - ", j.a, ")")
            if n != 1
                print(buf, "^", n)
            end
        end
        if x < 0
            print(buf, "(", x, ")")
        else
            print(buf, x)
        end
        String(take!(buf))
    end, " + "))
end

function domain_check(j::Jet, x)
    if j.a !== x
        throw(DomainError("Evaluation is only valid at a"))
    end
end
domain_check(j::Jet, x::ATB) = domain_check(j, primal(x))

function ChainRulesCore.rrule(::typeof(domain_check), j::Jet, x)
    domain_check(j, x), Δ->(ZeroTangent(), ZeroTangent(), ZeroTangent())
end

function (j::Jet)(x)
    domain_check(j, x)
    j[0]
end

function tderiv(bb::TaylorBundle{N}, i) where {N}
    i == 0 && return bb
    N-i == 0 && return bb.coeffs[i]
    TaylorBundle{N-i}(bb.coeffs[i], bb.coeffs[i+1:end])
end

function aldot(a::TaylorBundle{N}, b::TaylorBundle{M}) where {M,N}
    a.primal * b.primal +
        sum(1:min(M,N)) do i
            a.coeffs[i] * b.coeffs[i]
        end
end
aldot(a, b::TaylorBundle) = a * b.primal

*ₐₗ(a, b) = a * b
function *ₐₗ(a::Tangent{T}, b::TaylorBundle) where {N,B,T<:TaylorBundle{N,B}}
    # TODO: More general implementation of this
    bb = TaylorBundle{N}(a.primal, a.coeffs)
    cs = ntuple(1 + length(bb.coeffs)) do i
        aldot(tderiv(bb, i-1), b)
    end
    Tangent{T}(;primal = cs[1], coeffs = cs[2:end])
end

function ChainRulesCore.rrule(j::Jet, x)
    z = j(x)
    z, let djx = deriv(j)(x)
        function (Δ)
            (NoTangent(), Δ *ₐₗ djx)
        end
    end
end

function ChainRulesCore.rrule(::typeof(map), ::typeof(*), a, b)
    map(*, a, b), Δ->(NoTangent(), NoTangent(), map(*, Δ, b), map(*, a, Δ))
end

ChainRulesCore.rrule(::typeof(map), ::typeof(integrate), js::Array{<:Jet}) =
    map(integrate, js), Δ->(NoTangent(), NoTangent(), map(deriv, Δ))

struct derivBack
    js
end
(d::derivBack)(Δ::Union{ZeroTangent, NoTangent}) = (NoTangent(), NoTangent(), Δ)
(d::derivBack)(Δ::Array) = (NoTangent(), NoTangent(), broadcast(antideriv, d.js, Δ))

ChainRulesCore.rrule(::typeof(map), ::typeof(deriv), js::Array{<:Jet}) =
    map(deriv, js), derivBack(js)

ChainRulesCore.rrule(::typeof(map), ::typeof(antideriv), js::Array{<:Jet}, Δ) =
    map(antideriv, js, Δ), Δ->(NoTangent(), NoTangent(), map(deriv, Δ), One())

function mapev(js::Array{<:Jet}, xs::AbstractArray)
    map((j,x)->j(x), js, xs)
end

function ChainRulesCore.rrule(::typeof(mapev), js::Array{<:Jet}, xs::AbstractArray)
    mapev(js, xs), let djs=map(deriv, js)
        Δ->(NoTangent(), NoTangent(), map(*, unthunk(Δ), mapev(djs, xs)))
    end
end

function (∂⃖ₙ::∂⃖{N})(::typeof(map), f, a::Array) where {N}
    @assert Base.issingletontype(typeof(f))
    js = map(a) do x
        ∂f = ∂☆{N}()(ZeroBundle{N}(f),
                     TaylorBundle{N}(x,
                       (one(x), (zero(x) for i = 1:(N-1))...,)))
        @assert isa(∂f, TaylorBundle) || isa(∂f, TangentBundle{1})
        Jet{typeof(x), N}(x, ∂f.primal,
            isa(∂f, TangentBundle) ? ∂f.partials : ∂f.coeffs)
    end
    ∂⃖ₙ(mapev, js, a)
end

# Exempt these from the definition above
for f in (integrate, deriv, antideriv)
    function (∂⃖ₙ::∂⃖{N})(m::typeof(map), f::typeof(f), a::Array) where {N}
        invoke(∂⃖ₙ, Tuple{Any, Vararg{Any}}, m, f, a)
    end
end

"""
    jet_taylor_ev(::Val{}, jet, taylor)

Generates a closed form arithmetic expression for the N-th component
of the action of a 1d jet (of order at least N) on a maximally symmetric
(i.e. taylor) tangent bundle element. In particular, if we represent both
the `jet` and the `taylor` tangent bundle element by their associated canonical
taylor series:

    j = j₀ + j₁ (x - a) + j₂ 1/2 (x - a)^2 + ... + jₙ 1/n! (x - a)^n
    t = t₀ + t₁ (x - t₀) + t₂ 1/2 (x - t₀)^2 + ... + tₙ 1/n! (x - t₀)^n

then the action of evaluating `j` on `t`, is some other taylor series

    t′ = a + t′₁ (x - a) + t′₂ 1/2 (x - a)^2 + ... + t′ₙ 1/n! (x - a)^n

The t′ᵢ can be found by explicitly plugging in `t` for every `x` and expanding
out, dropping terms of orders that are higher. This computes closed form
expressions for the t′ᵢ that are hopefully easier on the compiler.
"""
@generated function jet_taylor_ev(::Val{N}, jet, taylor) where {N}
    elements = Any[nothing for _ = 1:N]
    for part in partitions(N)
        coeff = (factorial(N) ÷ (prod(factorial, values(countmap(part))) *
            prod(factorial, part)))
        factor = mapreduce((a,b)->:($a*$b), part) do i
            :(taylor[$i])
        end
        summand = coeff == 1 ? factor : :($coeff*$factor)
        let e = elements[length(part)]
            elements[length(part)] = e === nothing ? summand : :($summand + $e)
        end
    end
    return Expr(:call, +, map(enumerate(elements)) do (i, e)
        :(jet[$i]*$e)
    end...)
end

@generated function (j::Jet{T, N} where T)(x::TaylorBundle{M}) where {N, M}
    O = min(M,N)
    quote
        domain_check(j, x.primal)
        coeffs = x.coeffs
        TaylorBundle{$O}(j[0],
            ($((:(jet_taylor_ev(Val{$i}(), coeffs, j)) for i = 1:O)...),))
    end
end

function (j::Jet{T, 1} where T)(x::TangentBundle{1})
    domain_check(j, x.primal)
    coeffs = x.partials
    TangentBundle{1}(j[0], (jet_taylor_ev(Val{1}(), coeffs, j),))
end

function (j::Jet{T, N} where T)(x::TangentBundle{N, M}) where {N, M}
    error("TODO")
end


text:
# author: Tommaso Sartor

include("./riccati.jl")
using .riccati

# main script
length(ARGS) != 2 && error("Usage: julia test_riccati.jl <nmass> <nrep>")

nmass = parse(Int, ARGS[1])
nrep = parse(Int, ARGS[2])

nx = 2*nmass
nu = nmass
N = 5

# data

# mass spring system
Ts = 0.5

Ac = zeros(nx, nx)
for ii in 1:nmass
	Ac[ii, nmass+ii] = 1.0
end
for ii in 1:nmass
	Ac[nmass+ii, ii] = -2.0
end
for ii in 1:nmass-1
	Ac[nmass+ii+1, ii] = 1.0
end
for ii in 1:nmass-1
	Ac[nmass+ii, ii+1] = 1.0
end

Bc = zeros(nx, nu)
for ii in 1:nu
	Bc[nmass+ii, ii] = 1.0
end

#display(Ac)
#println()
#display(Bc)
#println()

MM = [ Ts*Ac Ts*Bc; zeros(nu, nx+nu) ]

#display(MM)
#println()

#MM = exp( MM )
MM = randn(nu+nx, nu+nx)

#display(MM)
#println()

A = MM[1:nx,1:nx]
B = MM[1:nx, nx+1:end]

#display(A)
#println()
#display(B)
#println()

Q = zeros(nx, nx)
for ii in 1:nx
	Q[ii, ii] = 1.0
end

R = zeros(nu, nu)
for ii in 1:nu
	R[ii, ii] = 2.0
end

#display(Q)
#println()
#display(R)
#println()

x0 = zeros(nx, 1)
x0[1] = 3.5
x0[2] = 3.5

#display(x0)
#println()

# work matrices

BAt = [transpose(B); transpose(A)]
#display(BAt)
#println()
RSQ = [R zeros(nu,nx); zeros(nx,nu) Q]

BAtL = zeros(nu+nx, nx)
L = zeros(nu+nx, nu+nx, N)
LN = zeros(nx, nx)
M = zeros(nu+nx, nu+nx)

# riccati recursion, square root algorithm

for rep in 1:nrep
	riccati_trf(N, nx, nu, BAt, RSQ, L, LN, BAtL, M)
end

time_start = time_ns()

for rep in 1:nrep
	riccati_trf(N, nx, nu, BAt, RSQ, L, LN, BAtL, M)
end

time_end = time_ns()

println((time_end-time_start)/nrep*1e-9)

#display(L)
#println()
#display(LN)
#println()

max_stars_repo_name: kokizzu/prometeo

max_issues_repo_name: kokizzu/prometeo

max_forks_repo_name: kokizzu/prometeo


text:
struct Result{T<:QuObject,A}
    times::Vector{Float64}
    states::Vector{T}
    evals::Matrix{ComplexF64}
    #probs::Vector{Matrix{Float64}}
    solver::A
end

function lsolve(L::Liouvillian,ψ₀::Ket,tspan,e_ops,alg;kwargs...)
    dimsmatch(L,ψ₀)
    prob = ODEProblem(L,complex(full(ψ₀)),tspan)
    sol  = solve(prob,alg;dense=false,abstol=1E-10,reltol=1E-8,kwargs...)
    states = Ket.(sol.u,[dims(ψ₀)])
    evals  = calc_expvals(e_ops,states)
    #probs  = levelprobs(states)
    return Result(sol.t,states,evals,sol.alg)
end

function lsolve_steady(L::Liouvillian,ψ₀::Ket,e_ops,alg;kwargs...)
    #alg = DynamicSS(odealg;dense=false,abstol=1E-10,reltol=1E-8,kwargs...))
    dimsmatch(L,ψ₀)
    f = ODEFunction(L,jac=(J,ψ,p,t)->L(Val{:jac},J,ψ,p,t))
    prob = SteadyStateProblem(f,complex(full(ψ₀)))
    sol  = solve(prob,alg;dense=false,abstol=1E-10,reltol=1E-8,kwargs...)
    states = [Ket(sol.u,dims(ψ₀))]
    evals  = calc_expvals(e_ops,states)
    #probs  = levelprobs(states)
    return Result(Float64[],states,evals,sol.alg)
end

function lsolve(L::Liouvillian,ρ₀::Operator,tspan,e_ops,alg;kwargs...)
    dimsmatch(L,ρ₀)
    prob = ODEProblem(L,vec(complex(full(ρ₀))),tspan)
    sol  = solve(prob,alg;dense=false,abstol=1E-10,reltol=1E-8,kwargs...)
    states = Operator.(unvec.(sol.u),[dims(ρ₀)])
    evals  = calc_expvals(e_ops,states)
    return Result(sol.t,states,evals,sol.alg)
end

function sesolve(H,ψ₀::Ket,tspan,e_ops=(),alg=Vern8();kwargs...)
    L = SchrodingerEvo(H)
    return lsolve(L,ψ₀,tspan,e_ops,alg;kwargs...)
end

function mesolve(H,C,ρ₀::Operator,tspan,e_ops=(),alg=Tsit5();kwargs...)
    L = LindbladEvo(H,C)
    return lsolve(L,ρ₀,tspan,e_ops,alg;kwargs...)
end

function mesolve(H,C,ψ₀::Ket,tspan,e_ops=(),alg=Tsit5();kwargs...)
    return mesolve(H,C,Operator(ψ₀),tspan,e_ops,alg;kwargs...)
end

function psolve(U::Propagator,ψ₀::Ket,steps,e_ops)
    dimsmatch(U,ψ₀)
    states = Vector{Ket{Vector{Complex{Float64}},length(dims(ψ₀))}}(undef,steps+1)
    states[1] = complex(dense(ψ₀))
    for s = 2:steps+1
        states[s] = U(states[s-1])
    end
    evals  = calc_expvals(e_ops,states)
    #probs  = levelprobs(states)
    t = collect(LinRange(0,steps*U.Δt,steps+1))
    return Result(t,states,evals,:SchrodingerPropSolver)
end

function psolve(U::Propagator,ρ₀::Operator,steps,e_ops)
    dimsmatch(U,ρ₀)
    states = Vector{Operator{Matrix{Complex{Float64}},length(dims(ρ₀))}}(undef,steps+1)
    states[1] = complex(dense(ρ₀))
    for s = 2:steps+1
        states[s] = U(states[s-1])
    end
    evals  = calc_expvals(e_ops,states)
    #probs  = levelprobs(states)
    t = collect(LinRange(0,steps*U.Δt,steps+1))
    return Result(t,states,evals,:SchrodingerPropSolver)
end

function psteady(U::Propagator,ρ₀::Operator,steps,e_ops)
    dimsmatch(U,ρ₀)
    states = Vector{Operator{Matrix{Complex{Float64}},length(dims(ρ₀))}}(undef,2)
    states[1] = complex(dense(ρ₀))
    states[2] = U(states[1],steps)
    evals  = calc_expvals(e_ops,states)
    #probs  = levelprobs(states)
    return Result([0.0,steps*U.Δt],states,evals,:SchrodingerPropSteady)
end

calc_expvals(o::Operator,states) = calc_expvals((o,),states)
function calc_expvals(e_ops,states)
    isempty(e_ops) && return Matrix{ComplexF64}(undef,0,0)
    M = length(e_ops); N = length(states)
    expvals = Matrix{ComplexF64}(undef,N,M)
    for (j,σ) in enumerate(e_ops)
        eval = view(expvals,:,j)
        expect!(eval, σ, states)
    end
    return expvals
end


text:
__precompile__(true)
"""
ChrBase package

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones, and contributors to julia
Licensed under MIT License, see LICENSE.md
In part based on code for Char in Julia
"""
module ChrBase

using ModuleInterfaceTools

@api extend! CharSetEncodings

@api public Chr

@api public! is_malformed, is_overlong

@api develop get_utf8_2, get_utf8_3, get_utf8_4, get_utf8_16, get_utf8_32, get_utf16, get_utf16_32,
             is_valid_continuation, is_surrogate_lead, is_surrogate_trail, is_surrogate_codeunit,
             LatinChars, ByteChars, WideChars, AbsChar, bytoff, chroff, chrdiff, utf_trail,
             codepoint_cse, codepoint_rng, codepoint_adj, utf8proc_error,
             write_utf8, write_utf16, _write_utf8_2, _write_utf8_3, _write_utf8_4, _write_ucs2,
             _lowercase_l, _uppercase_l,
             _is_lower_a, _is_lower_l, _is_lower_al, _is_lower_ch,
             _is_upper_a, _is_upper_l, _is_upper_al, _is_upper_ch

@api develop! _isvalid_chr

include("core.jl")
@static V6_COMPAT && include("compat.jl")
include("CaseTables.jl"); using .CaseTables
include("casefold.jl")
include("io.jl")
include("traits.jl")
include("unicode.jl")
include("utf8proc.jl")

@api freeze

end # module ChrBase


text:
# This file was generated, do not modify it. # hide
train = 1:findlast(X.Year .< 2005);
test = last(train)+1:length(y)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


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

"""
    Arrow.ArrowVector

An abstract type that subtypes `AbstractVector`. Each specific arrow array type
subtypes `ArrowVector`. See [`BoolVector`](@ref), [`Primitive`](@ref), [`List`](@ref),
[`Map`](@ref), [`FixedSizeList`](@ref), [`Struct`](@ref), [`DenseUnion`](@ref),
[`SparseUnion`](@ref), and [`DictEncoded`](@ref) for more details.
"""
abstract type ArrowVector{T} <: AbstractVector{T} end

Base.IndexStyle(::Type{A}) where {A <: ArrowVector} = Base.IndexLinear()
Base.similar(::Type{A}, dims::Dims) where {T, A <: ArrowVector{T}} = Vector{T}(undef, dims)
validitybitmap(x::ArrowVector) = x.validity
nullcount(x::ArrowVector) = validitybitmap(x).nc
getmetadata(x::ArrowVector) = x.metadata

function toarrowvector(x, i=1, de=Dict{Int64, Any}(), ded=DictEncoding[], meta=getmetadata(x); compression::Union{Nothing, LZ4FrameCompressor, ZstdCompressor}=nothing, kw...)
    @debug 2 "converting top-level column to arrow format: col = $(typeof(x)), compression = $compression, kw = $(kw.data)"
    @debug 3 x
    A = arrowvector(x, i, 0, 0, de, ded, meta; compression=compression, kw...)
    if compression isa LZ4FrameCompressor
        A = compress(Meta.CompressionType.LZ4_FRAME, compression, A)
    elseif compression isa ZstdCompressor
        A = compress(Meta.CompressionType.ZSTD, compression, A)
    end
    @debug 2 "converted top-level column to arrow format: $(typeof(A))"
    @debug 3 A
    return A
end

function arrowvector(x, i, nl, fi, de, ded, meta; dictencoding::Bool=false, dictencode::Bool=false, kw...)
    if !(x isa DictEncode) && !dictencoding && (dictencode || (x isa AbstractArray && DataAPI.refarray(x) !== x))
        x = DictEncode(x, dictencodeid(i, nl, fi))
    end
    S = maybemissing(eltype(x))
    return arrowvector(S, x, i, nl, fi, de, ded, meta; dictencode=dictencode, kw...)
end

# conversions to arrow types
arrowvector(::Type{Dates.Date}, x, i, nl, fi, de, ded, meta; kw...) =
    arrowvector(converter(DATE, x), i, nl, fi, de, ded, meta; kw...)
arrowvector(::Type{Dates.Time}, x, i, nl, fi, de, ded, meta; kw...) =
    arrowvector(converter(TIME, x), i, nl, fi, de, ded, meta; kw...)
arrowvector(::Type{Dates.DateTime}, x, i, nl, fi, de, ded, meta; kw...) =
    arrowvector(converter(DATETIME, x), i, nl, fi, de, ded, meta; kw...)
arrowvector(::Type{ZonedDateTime}, x, i, nl, fi, de, ded, meta; kw...) =
    arrowvector(converter(Timestamp{Meta.TimeUnit.MILLISECOND, Symbol(x[1].timezone)}, x), i, nl, fi, de, ded, meta; kw...)
arrowvector(::Type{P}, x, i, nl, fi, de, ded, meta; kw...) where {P <: Dates.Period} =
    arrowvector(converter(Duration{arrowperiodtype(P)}, x), i, nl, fi, de, ded, meta; kw...)

# fallback that calls ArrowType
function arrowvector(::Type{S}, x, i, nl, fi, de, ded, meta; kw...) where {S}
    if ArrowTypes.istyperegistered(S)
        meta = meta === nothing ? Dict{String, String}() : meta
        arrowtype = ArrowTypes.getarrowtype!(meta, S)
        if arrowtype === S
            return arrowvector(ArrowType(S), x, i, nl, fi, de, ded, meta; kw...)
        else
            return arrowvector(converter(arrowtype, x), i, nl, fi, de, ded, meta; kw...)
        end
    end
    return arrowvector(ArrowType(S), x, i, nl, fi, de, ded, meta; kw...)
end

arrowvector(::NullType, x, i, nl, fi, de, ded, meta; kw...) = MissingVector(length(x))
compress(Z::Meta.CompressionType, comp, v::MissingVector) =
    Compressed{Z, MissingVector}(v, CompressedBuffer[], length(v), length(v), Compressed[])

function makenodesbuffers!(col::MissingVector, fieldnodes, fieldbuffers, bufferoffset, alignment)
    push!(fieldnodes, FieldNode(length(col), length(col)))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    return bufferoffset
end

function writebuffer(io, col::MissingVector, alignment)
    return
end

"""
    Arrow.ValidityBitmap

A bit-packed array type where each bit corresponds to an element in an
[`ArrowVector`](@ref), indicating whether that element is "valid" (bit == 1),
or not (bit == 0). Used to indicate element missingness (whether it's null).

If the null count of an array is zero, the `ValidityBitmap` will be "emtpy"
and all elements are treated as "valid"/non-null.
"""
struct ValidityBitmap <: ArrowVector{Bool}
    bytes::Vector{UInt8} # arrow memory blob
    pos::Int # starting byte of validity bitmap
    ℓ::Int # # of _elements_ (not bytes!) in bitmap (because bitpacking)
    nc::Int # null count
end

Base.size(p::ValidityBitmap) = (p.ℓ,)
nullcount(x::ValidityBitmap) = x.nc

function ValidityBitmap(x)
    T = eltype(x)
    if !(T >: Missing)
        return ValidityBitmap(UInt8[], 1, length(x), 0)
    end
    len = length(x)
    blen = cld(len, 8)
    bytes = Vector{UInt8}(undef, blen)
    st = iterate(x)
    nc = 0
    b = 0xff
    j = k = 1
    for y in x
        if y === missing
            nc += 1
            b = setbit(b, false, j)
        end
        j += 1
        if j == 9
            @inbounds bytes[k] = b
            b = 0xff
            j = 1
            k += 1
        end
    end
    if j > 1
        bytes[k] = b
    end
    return ValidityBitmap(nc == 0 ? UInt8[] : bytes, 1, nc == 0 ? 0 : len, nc)
end

@propagate_inbounds function Base.getindex(p::ValidityBitmap, i::Integer)
    # no boundscheck because parent array should do it
    # if a validity bitmap is empty, it either means:
    #   1) the parent array null_count is 0, so all elements are valid
    #   2) parent array is also empty, so "all" elements are valid
    p.nc == 0 && return true
    # translate element index to bitpacked byte index
    a, b = fldmod1(i, 8)
    @inbounds byte = p.bytes[p.pos + a - 1]
    # check individual bit of byte
    return getbit(byte, b)
end

@propagate_inbounds function Base.setindex!(p::ValidityBitmap, v, i::Integer)
    x = convert(Bool, v)
    p.ℓ == 0 && !x && throw(BoundsError(p, i))
    a, b = fldmod1(i, 8)
    @inbounds byte = p.bytes[p.pos + a - 1]
    @inbounds p.bytes[p.pos + a - 1] = setbit(byte, x, b)
    return v
end

function writebitmap(io, col::ArrowVector, alignment)
    v = col.validity
    @debug 1 "writing validity bitmap: nc = $(v.nc), n = $(cld(v.ℓ, 8))"
    v.nc == 0 && return 0
    n = Base.write(io, view(v.bytes, v.pos:(v.pos + cld(v.ℓ, 8) - 1)))
    return n + writezeros(io, paddinglength(n, alignment))
end

include("compressed.jl")
include("primitive.jl")
include("bool.jl")
include("list.jl")
include("fixedsizelist.jl")
include("map.jl")
include("struct.jl")
include("unions.jl")
include("dictencoding.jl")


text:
doc"""
    Categorical(p)
A *Categorical distribution* is parameterized by a probability vector `p` (of length `K`).
$P(X = k) = p[k]  \quad \text{for } k = 1, 2, \ldots, K.$
```julia
Categorical(p)   # Categorical distribution with probability vector p
params(d)        # Get the parameters, i.e. (p,)
probs(d)         # Get the probability vector, i.e. p
ncategories(d)   # Get the number of categories, i.e. K
```
Here, `p` must be a real vector, of which all components are nonnegative and sum to one.
**Note:** The input vector `p` is directly used as a field of the constructed distribution, without being copied.
External links:
* [Categorical distribution on Wikipedia](http://en.wikipedia.org/wiki/Categorical_distribution)
"""

immutable Categorical{T<:Real} <: DiscreteUnivariateDistribution
    K::Int
    p::Vector{T}

    Categorical(p::Vector{T}, ::NoArgCheck) = new(length(p), p)

    function Categorical(p::Vector{T})
        @check_args(Categorical, isprobvec(p))
        new(length(p), p)
    end

    function Categorical(k::Integer)
        @check_args(Categorical, k >= 1)
        new(k, fill(1/k, k))
    end
end

Categorical{T<:Real}(p::Vector{T}, ::NoArgCheck) = Categorical{T}(p, NoArgCheck())
Categorical{T<:Real}(p::Vector{T}) = Categorical{T}(p)
Categorical(k::Integer) = Categorical{Float64}(k)

@distr_support Categorical 1 d.K

### Conversions

convert{T<:Real, S<:Real}(::Type{Categorical{T}}, p::Vector{S}) = Categorical(Vector{T}(p))
convert{T<:Real, S<:Real}(::Type{Categorical{T}}, d::Categorical{S}) = Categorical(Vector{T}(d.p))

### Parameters

ncategories(d::Categorical) = d.K
probs(d::Categorical) = d.p
params(d::Categorical) = (d.p,)
@inline partype{T<:Real}(d::Categorical{T}) = T

### Statistics

function categorical_mean{T<:Real}(p::AbstractArray{T})
    k = length(p)
    s = zero(T)
    for i = 1:k
        @inbounds s += p[i] * i
    end
    s
end

mean(d::Categorical) = categorical_mean(d.p)

function median{T<:Real}(d::Categorical{T})
    k = ncategories(d)
    p = probs(d)
    cp = zero(T)
    i = 0
    while cp < 1/2 && i <= k
        i += 1
        @inbounds cp += p[i]
    end
    i
end

function var{T<:Real}(d::Categorical{T})
    k = ncategories(d)
    p = probs(d)
    m = categorical_mean(p)
    s = zero(T)
    for i = 1:k
        @inbounds s += abs2(i - m) * p[i]
    end
    s
end

function skewness{T<:Real}(d::Categorical{T})
    k = ncategories(d)
    p = probs(d)
    m = categorical_mean(p)
    s = zero(T)
    for i = 1:k
        @inbounds s += (i - m)^3 * p[i]
    end
    v = var(d)
    s / (v * sqrt(v))
end

function kurtosis{T<:Real}(d::Categorical{T})
    k = ncategories(d)
    p = probs(d)
    m = categorical_mean(p)
    s = zero(T)
    for i = 1:k
        @inbounds s += (i - m)^4 * p[i]
    end
    s / abs2(var(d)) - 3
end

entropy(d::Categorical) = entropy(d.p)

function mgf{T<:Real}(d::Categorical{T}, t::Real)
    k = ncategories(d)
    p = probs(d)
    s = zero(T)
    for i = 1:k
        @inbounds s += p[i] * exp(t)
    end
    s
end

function cf{T<:Real}(d::Categorical{T}, t::Real)
    k = ncategories(d)
    p = probs(d)
    s = zero(T) + zero(T)*im
    for i = 1:k
        @inbounds s += p[i] * cis(t)
    end
    s
end

mode(d::Categorical) = indmax(probs(d))

function modes(d::Categorical)
    K = ncategories(d)
    p = probs(d)
    maxp = maximum(p)
    r = Array(Int, 0)
    for k = 1:K
        @inbounds if p[k] == maxp
            push!(r, k)
        end
    end
    r
end

### Evaluation

function cdf{T<:Real}(d::Categorical{T}, x::Int)
    k = ncategories(d)
    p = probs(d)
    x < 1 && return zero(T)
    x >= k && return one(T)
    c = p[1]
    for i = 2:x
        @inbounds c += p[i]
    end
    return c
end

pdf{T<:Real}(d::Categorical{T}, x::Int) = insupport(d, x) ? d.p[x] : zero(T)

logpdf(d::Categorical, x::Int) = insupport(d, x) ? log(d.p[x]) : -Inf

pdf(d::Categorical) = copy(d.p)

function _pdf!{T<:Real}(r::AbstractArray, d::Categorical{T}, rgn::UnitRange)
    vfirst = round(Int, first(rgn))
    vlast = round(Int, last(rgn))
    vl = max(vfirst, 1)
    vr = min(vlast, d.K)
    p = probs(d)
    if vl > vfirst
        for i = 1:(vl - vfirst)
            r[i] = zero(T)
        end
    end
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = p[v]
    end
    if vr < vlast
        for i = (vr - vfirst + 2):length(rgn)
            r[i] = zero(T)
        end
    end
    return r
end

function quantile(d::Categorical, p::Float64)
    0 <= p <= 1 || throw(DomainError())
    k = ncategories(d)
    pv = probs(d)
    i = 1
    v = pv[1]
    while v < p && i < k
        i += 1
        @inbounds v += pv[i]
    end
    i
end

# sampling

sampler(d::Categorical) = AliasTable(d.p)

### sufficient statistics

immutable CategoricalStats <: SufficientStats
    h::Vector{Float64}
end

function add_categorical_counts!{T<:Integer}(h::Vector{Float64}, x::AbstractArray{T})
    for i = 1 : length(x)
        @inbounds xi = x[i]
        h[xi] += 1.   # cannot use @inbounds, as no guarantee that x[i] is in bound
    end
    h
end

function add_categorical_counts!{T<:Integer}(h::Vector{Float64}, x::AbstractArray{T}, w::AbstractArray{Float64})
    n = length(x)
    if n != length(w)
        throw(ArgumentError("Inconsistent array lengths."))
    end
    for i = 1 : n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        h[xi] += wi   # cannot use @inbounds, as no guarantee that x[i] is in bound
    end
    h
end

function suffstats{T<:Integer}(::Type{Categorical}, k::Int, x::AbstractArray{T})
    CategoricalStats(add_categorical_counts!(zeros(k), x))
end

function suffstats{T<:Integer}(::Type{Categorical}, k::Int, x::AbstractArray{T}, w::AbstractArray{Float64})
    CategoricalStats(add_categorical_counts!(zeros(k), x, w))
end

typealias CategoricalData Tuple{Int, AbstractArray}

suffstats(::Type{Categorical}, data::CategoricalData) = suffstats(Categorical, data...)
suffstats(::Type{Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = suffstats(Categorical, data..., w)

# Model fitting

function fit_mle(::Type{Categorical}, ss::CategoricalStats)
    Categorical(pnormalize!(ss.h))
end

function fit_mle{T<:Integer}(::Type{Categorical}, k::Integer, x::AbstractArray{T})
    Categorical(pnormalize!(add_categorical_counts!(zeros(k), x)), NoArgCheck())
end

function fit_mle{T<:Integer}(::Type{Categorical}, k::Integer, x::AbstractArray{T}, w::AbstractArray{Float64})
    Categorical(pnormalize!(add_categorical_counts!(zeros(k), x, w)), NoArgCheck())
end

fit_mle(::Type{Categorical}, data::CategoricalData) = fit_mle(Categorical, data...)
fit_mle(::Type{Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = fit_mle(Categorical, data..., w)

fit_mle{T<:Integer}(::Type{Categorical}, x::AbstractArray{T}) = fit_mle(Categorical, maximum(x), x)
fit_mle{T<:Integer}(::Type{Categorical}, x::AbstractArray{T}, w::AbstractArray{Float64}) = fit_mle(Categorical, maximum(x), x, w)

fit(::Type{Categorical}, data::CategoricalData) = fit_mle(Categorical, data)
fit(::Type{Categorical}, data::CategoricalData, w::AbstractArray{Float64}) = fit_mle(Categorical, data, w)


text:
using Test

@testset "SAC tests" begin
    include("construction.jl")
    include("io.jl")
    include("operations.jl")
    include("diff.jl")
    include("integrate.jl")
    include("stack.jl")
    include("great_circle.jl")
    include("rotate_to_gcp.jl")
    include("statistics.jl")
    include("util.jl")
end


text:
@testset "pop" begin
    x = Index{:a}(UnitMRange(1, 10), UnitMRange(1, 10))
    y = Index{:a}(UnitMRange(2, 10), UnitMRange(2, 10))
    z = Index{:a}(UnitMRange(1, 9), UnitMRange(1, 9))
    @test popfirst(x) == y
    @test pop(x) == z

    popfirst!(x)
    @test x == y

    x = Index{:a}(UnitMRange(1, 10), UnitMRange(1, 10))
    pop!(x)
    @test x == z

    x = IArray(collect(1:10))
    y = IArray(collect(2:10))
    z = IArray(collect(1:9))

    @test popfirst(x) == y
    @test pop(x) == z

    popfirst!(x)
    @test x == y

    x = IArray(collect(1:10))
    pop!(x)
    @test x == z
end


text:
using Symbolics, Test

@sym x y z m ω t μ

@testset "   Euler-Lagrange Solver" begin
    function Γ(w)
       t -> UpTuple([t, w(t), D(w)(t)])
    end

    function Lagrange_Equations(L)
        w -> D(∂(3)(L)∘Γ(w)) - ∂(2)(L)∘Γ(w)
    end

    function L_SHO(local_tuple::UpTuple)
        t, q, qdot = local_tuple.data
        (0.5m)*qdot^2 - (0.5m*ω^2)*q^2
    end

    function L_free(local_tuple::UpTuple)
        t, q, qdot = local_tuple.data
        (0.5m)*qdot^2
    end

    function L_pendulum(local_tuple::UpTuple)
        t, q, qdot = local_tuple.data
        m*qdot^2/2+cos(q)
    end

    @test Lagrange_Equations(L_SHO)(x)(t)  == (D(D(x)))(t) * m + (x)(t) * m * ω ^ 2
    @test Lagrange_Equations(L_free)(x)(t) == (D(D(x)))(t) * m
    @test Lagrange_Equations(L_pendulum)(x)(t) == (D(D(x)))(t) * m + sin((x)(t))

    # Lagrangian for the free particle in space (3D)
    function L_free_3d(local_tuple::UpTuple)
        t, q, qdot = local_tuple.data
        (qdot[1]^2+qdot[2]^2+qdot[3]^2)/2
    end

    q = UpTuple([x, y, z])

    lagrange_eqs_fp3d = Lagrange_Equations(L_free_3d)(q)
    lagrange_eqs_fp3d_t = Lagrange_Equations(L_free_3d)(q)(t)
    @test typeof(lagrange_eqs_fp3d_t) == UpTuple
    fp3d_eqs = [(D(D(x)))(t), (D(D(y)))(t), (D(D(z)))(t)] # eqs of motion for the free particle in space
    @test lagrange_eqs_fp3d_t.data == fp3d_eqs
    @test lagrange_eqs_fp3d_t == up(fp3d_eqs)
    @test lagrange_eqs_fp3d_t == UpTuple(fp3d_eqs)
    @test lagrange_eqs_fp3d(t).data == fp3d_eqs
    @test lagrange_eqs_fp3d(t) == up(fp3d_eqs)
    @test lagrange_eqs_fp3d(t) == UpTuple(fp3d_eqs)

    # Lagrangian for the Kepler problem
    function L_kepler(local_tuple::UpTuple)
        t, q, qdot = local_tuple.data
        r = (q[1]^2+q[2]^2)^(1/2)
        (qdot[1]^2+qdot[2]^2)*(m/2)+μ/r
    end

    q = UpTuple([x, y])

    # Equivalent to Eq. (1.48) of SICM
    lagrange_eqs_kepler = Lagrange_Equations(L_kepler)(q)
    lagrange_eqs_kepler_t = Lagrange_Equations(L_kepler)(q)(t)
    @test typeof(lagrange_eqs_kepler_t) == UpTuple
    kepler_eqs = [m * (D(D(x)))(t) + ((x)(t) ^ 2 + (y)(t) ^ 2) ^ -1.5 * (x)(t) * μ, m * (D(D(y)))(t) + ((x)(t) ^ 2 + (y)(t) ^ 2) ^ -1.5 * (y)(t) * μ]
    @test lagrange_eqs_kepler_t.data == kepler_eqs
    @test lagrange_eqs_kepler_t == up(kepler_eqs)
    @test lagrange_eqs_kepler_t == UpTuple(kepler_eqs)
    @test lagrange_eqs_kepler(t).data == kepler_eqs
    @test lagrange_eqs_kepler(t) == up(kepler_eqs)
    @test lagrange_eqs_kepler(t) == UpTuple(kepler_eqs)

end


text:
# File used to benchmark the performance of individual functions and different versions of those functions to improve solver performance
# Useful for developers readOFBoundaryFile

# Install these packages first
using Profile
using ProfileView #Doesn't want to install on my work ubuntu desktop for some reason
using BenchmarkTools

include("shockTube.jl")
include("FvCFD.jl")

# Create mesh
meshPath = "test/OFshockTube_400"
OFmesh = OpenFOAMMesh(meshPath)
nCells = size(OFmesh.cells, 1)
mesh, cellPrimitives = initializeShockTube3DFVM(nCells...)
boundaryConditions = [ zeroGradientBoundary, [], emptyBoundary, [] ]
nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(OFmesh)

# Create solution state
nVars = 5
nFluxes = 3*nVars
cellState = encodePrimitives3D(cellPrimitives)
cellFluxes = zeros(nCells, nFluxes)
fluxResiduals = zeros(nCells, nVars)
faceFluxes = zeros(nFaces, nFluxes)
sln = SolutionState(cellState, cellFluxes, cellPrimitives, fluxResiduals, faceFluxes)
decodeSolution_3D(sln)

# Apply zero face fluxes to empty faces
zeroFlux = zeros(nFluxes)
for b in 1:nBoundaries
    if boundaryConditions[2*b-1] == emptyBoundary
        for f in OFmesh.boundaryFaces[b]
            sln.faceFluxes[f,:] = zeroFlux
        end
    end
end

#### greenGaussGrad_matrix ####
# # Create test data
# fDeltas = faceDeltas(OFmesh, solutionState)
#
# # Test
# greenGaussGrad_matrix(OFmesh, fDeltas, false)
# function run1000()
#     for i in 1:1000
#         greenGaussGrad_matrix(OFmesh, fDeltas, false)
#     end
# end
# @profview run1000()
# @benchmark greenGaussGrad_matrix(OFmesh, fDeltas, false)

#### greenGaussGrad ####
P = sln.cellPrimitives[:,1]
P = reshape(P, nCells, :)
function run1000()
    for i in 1:1000
        greenGaussGrad(OFmesh, P)
    end
end
# @profview run1000()
output = zeros(400, 1, 3)
@benchmark greenGaussGrad(OFmesh, P)
# output = zeros(2001, 1)
# @benchmark linInterp_3D(OFmesh, P, output)
# @code_warntype greenGaussGrad(OFmesh, P, false)

#### LinInterp ####
# @code_warntype linInterp_3D(mesh, sln.cellFluxes, sln.faceFluxes)

#### Entire Code ####
# @profview solve(OFmesh, meshPath, cellPrimitives, boundaryConditions, initDt=0.00001, endTime=0.14267, outputInterval=0.14267, targetCFL=0.1, silent=true)
# @btime solve(OFmesh, meshPath, cellPrimitives, boundaryConditions, initDt=0.00001, endTime=0.005, outputInterval=0.14267, targetCFL=0.01, silent=true, createRestartFile=false, createVTKOutput=false)
# @code_warntype solve(OFmesh, meshPath, cellPrimitives, boundaryConditions, initDt=0.00001, endTime=0.00001, outputInterval=0.14267, targetCFL=0.01, silent=true, createRestartFile=false, createVTKOutput=false)

max_stars_repo_name: henrystoldt/fvCFD


text:
using Test
using Delta

@testset "checking gradient" begin
    include("./checkgrad/0-softmax.jl")
    include("./checkgrad/1-pool.jl")
    include("./checkgrad/2-linear.jl")
    include("./checkgrad/3-mlp.jl")
    include("./checkgrad/4-chain.jl")
    include("./checkgrad/5-conv1d.jl")
    include("./checkgrad/6-addingProblem.jl")
    include("./checkgrad/7-scaler.jl")
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Part 2: another PDE
if !@isdefined USE_GPU
    const USE_GPU = length(ARGS) > 0 && ARGS[1] == "gpu"
end
using JLD
using SparseArrays
using SuiteSparse
using LinearAlgebra
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
if !ParallelStencil.is_initialized()
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 2)
    else
        @init_parallel_stencil(Threads, Float64, 2)
    end
end

include("./part2_utils.jl")
include("./multigrid.jl")

""" Initial condition (for temperature T or vorticity W). """
@enum Init_t begin
    cosine       # Initalize field as cosine
    random       # Initialize field as random
    W_from_file  # Read initial condition from file (for W, only for debugging)
end

""" Input parameters for Navier-Stokes simulation. """
mutable struct SimIn_t
    k               :: Float64  # thermal diffusivity (always set k = 1.0 since we use the non-dimensionalized Navier-Stokes equations)
    Ra              :: Float64  # Rayleigh number
    Pr              :: Float64  # Prandtl number
    nx              :: Int      # number of grid-points in x-direction
    ny              :: Int      # number of grid-points in y-direction
    ttot            :: Float64  # total time for simulation
    beta            :: Float64  # explicit (beta = 0.0) / implicit (beta = 1.0) / semi-implicit (beta = 0.5)
    niters          :: Int      # max number of iterations for iterative solver
    tol             :: Float64  # tolerance for linear solver
    a_dif           :: Float64  # diffusive timestep parameter
    a_adv           :: Float64  # advective timestep parameter
    T_init_strategy :: Init_t   # initialization strategy for T
    W_init_strategy :: Init_t   # initialization strategy for W

    SimIn_t() = new(1.0, 1.0e6, 1.0e-3, 257, 65, 0.1, 0.0, 50, 1.0e-3, 0.15, 0.4, cosine, random)
end

""" Output fields of Navier-Stokes simulation. """
struct SimOut_t
    T :: Matrix{Float64}    # Temperature
    W :: Matrix{Float64}    # Vorticity
    S :: Matrix{Float64}    # Streamfunction
    t_elapsed :: Float64    # Computation time
    timed_iters :: Float64  # Number of timesteps timed (i.e. minus warmup)
end

""" Initializes array M (i.e. T or W) given initialization scheme, mesh size h, and width. """
@views function init_array!(M, scheme::Init_t, h, width)
    nx, ny = size(M)
    if scheme == cosine
        for i in 1:nx
            M[i, :] .= 0.5 * (1.0 + cos((3.0*π*(i-1)*h)/width))
        end
    elseif scheme == random
        M .= @rand(nx, ny)
    elseif scheme == W_from_file # only for testing against fortran
        M .= open(joinpath(@__DIR__, "..", "test", "reftest-files", "fortran", "Winit.bin")) do io
            return load(io)
        end
    else
        error()
    end
end

""" Computes dt for explicit or implicit timestepping scheme. """
@views function compute_dt(v, vx, vy, dt_dif, a_dif, a_adv, h, beta)
    v_max = maximum(v)
    if (v_max == 0)
        dt = dt_dif
    else
        # compute advective timestep
        dt_adv = a_adv * min(h / maximum(abs.(vx)), h / maximum(abs.(vy)))
        # compute timestep
        dt = (beta >= 0.5 ? dt_adv : min(dt_dif, dt_adv))
    end
    return dt
end

""" Computes velocity given streamfunction. """
@parallel_indices (ix, iy) function compute_velocity!(S, hx, hy, vx, vy)
    if (1 < ix < size(S, 1) && 1 < iy < size(S, 2))
        vx[ix, iy] =  ( S[ix, iy+1] - S[ix, iy-1] ) / ( 2 * hy )
        vy[ix, iy] = -( S[ix+1, iy] - S[ix-1, iy] ) / ( 2 * hx )
    end
    return nothing
end

""" Computes Ra * (dT / dx) term, used in vorticity update. """
@parallel_indices (ix, iy) function compute_Ra_dTdx!(Ra, hx, T, Ra_dTdx)
    if (1 < ix < size(T, 1) && 1 < iy < size(T, 2))
        Ra_dTdx[ix, iy] =  Ra * ( T[ix+1, iy] - T[ix-1, iy] ) / ( 2 * hx )
    end
    return nothing
end

""" Computes diffusion of 2D field. """
@parallel_indices (ix, iy) function compute_diffusion2d!(T, hx, hy, k, dT2)
    if (1 < ix < size(T, 1) && 1 < iy < size(T, 2))
        dT2[ix, iy] = k * (( T[ix+1, iy] - 2*T[ix, iy]+ T[ix-1, iy] ) / hx^2 +
                           ( T[ix, iy+1] - 2*T[ix, iy]+ T[ix, iy-1] ) / hy^2)
    end
    return nothing
end

""" Computes advection of 2D field in x-dimension (with upwinding). """
@parallel_indices (ix, iy) function compute_advection2d_x!(T, hx, vx, dTx)
    if (1 < ix < size(T, 1) && 1 < iy < size(T, 2))
        if vx[ix, iy] > 0
            dTx[ix, iy] = vx[ix, iy] * (T[ix, iy] - T[ix-1, iy]) / hx
        else
            dTx[ix, iy] = vx[ix, iy] * (T[ix+1, iy] - T[ix, iy]) / hx
        end
    end
    return nothing
end

""" Computes advection of 2D field in y-dimension (with upwinding). """
@parallel_indices (ix, iy) function compute_advection2d_y!(T, hy, vy, dTy)
    if (1 < ix < size(T, 1) && 1 < iy < size(T, 2))
        if vy[ix, iy] > 0
            dTy[ix, iy] = vy[ix, iy] * (T[ix, iy] - T[ix, iy-1]) / hy
        else
            dTy[ix, iy] = vy[ix, iy] * (T[ix, iy+1] - T[ix, iy]) / hy
        end
    end
    return nothing
end

""" Runs Navier-Stokes simulation in 2D. """
@views function navier_stokes_2D(; opt :: SimIn_t = SimIn_t(), verbose=true, do_vis=false, testmode=false)
    nx, ny = opt.nx, opt.ny

    # derived numerics
    h       = 1.0 / (ny - 1.0)
    width   = (nx - 1.0) / (ny - 1.0)
    dt_dif  = (opt.a_dif * min(h, h)^2) / max(opt.k, opt.Pr)

    # allocate field arrays
    S       = @zeros(nx, ny)
    vx      = @zeros(nx, ny)
    vy      = @zeros(nx, ny)
    v       = @zeros(nx, ny)
    T       = @zeros(nx, ny)
    dT2     = @zeros(nx, ny)
    dTx     = @zeros(nx, ny)
    dTy     = @zeros(nx, ny)
    T_rhs   = @zeros(nx, ny)
    W       = @zeros(nx, ny)
    dW2     = @zeros(nx, ny)
    dWx     = @zeros(nx, ny)
    dWy     = @zeros(nx, ny)
    W_rhs   = @zeros(nx, ny)
    Ra_dTdx = @zeros(nx, ny)

    init_array!(T, opt.T_init_strategy, h, width)
    init_array!(W, opt.W_init_strategy, h, width)

    # preallocate all the buffers necessary for multigrid
    prealloc_dict = preallocate_buffers(nx, ny)

    hx, hy = h, h

    # plotting
    T_storage = Vector{Matrix{Float64}}()
    W_storage = Vector{Matrix{Float64}}()
    S_storage = Vector{Matrix{Float64}}()

    tic = 0.0
    sim_time = 0.0
    step = 0
    while (sim_time < opt.ttot)
        if step == 3
            tic = time()
        end

        # solve for streamfunction S: D S = W (Dirichlet BCs = 0)
        r_rms = MGsolve_2DPoisson!(S, W, h, 0.0, opt.tol, opt.niters, false, prealloc_dict=prealloc_dict)

        # compute velocity field (vx, vy) from streamfunction S
        @parallel compute_velocity!(S, hx, hy, vx, vy)

        # compute velocity magnitude v
        v = sqrt.(vx.^2 + vy.^2)

        # compute timestep dt
        dt = compute_dt(v, vx, vy, dt_dif, opt.a_dif, opt.a_adv, h, opt.beta)

        # apply boundary conditions
        apply_boundary_conditions!(T)

        # compute Ra * dT / dx
        @parallel compute_Ra_dTdx!(opt.Ra, hx, T, Ra_dTdx)

        # diffusion terms for temperature T and vorticity W
        if (opt.beta ≉ 1.0)
            @parallel compute_diffusion2d!(T, hx, hy, opt.k, dT2)
            @parallel compute_diffusion2d!(W, hx, hy, opt.Pr, dW2)
        end

        # advection terms for temperature T and vorticity W
        @parallel compute_advection2d_x!(T, hx, vx, dTx)
        @parallel compute_advection2d_y!(T, hy, vy, dTy)
        @parallel compute_advection2d_x!(W, hx, vx, dWx)
        @parallel compute_advection2d_y!(W, hy, vy, dWy)

        # Euler step for temperature T and vorticity W
        if (opt.beta > 0.0)
            # semi-implicit step for temperature T
            c = 1.0 / (opt.beta * dt)
            T_rhs .= -c * ( T + dt * ( (1.0 - opt.beta) * dT2 - dTx - dTy ) )
            r_rms = MGsolve_2DPoisson!(T, T_rhs, h, c, opt.tol, opt.niters, true, prealloc_dict=prealloc_dict)

            # semi-implicit step for vorticity W
            c = c / opt.Pr
            W_rhs .= -c * ( W + dt * ( (1.0 - opt.beta) * dW2 - dWx - dWy - opt.Pr * Ra_dTdx ) )
            r_rms = MGsolve_2DPoisson!(W, W_rhs, h, c, opt.tol, opt.niters, false, prealloc_dict=prealloc_dict)
        else
            # explicit step for temperature T and vorticity W
            T .= T + dt * ( dT2 - dTx - dTy )
            W .= W + dt * ( dW2 - dWx - dWy - opt.Pr * Ra_dTdx )
        end

        sim_time += dt
        step += 1

        if ((step-1) % 20 == 0)
            if verbose
                println("time, step: $(sim_time) $(step)")
            end
            if do_vis
                push!(T_storage, copy(Array(T)))
                push!(W_storage, copy(Array(W)))
                push!(S_storage, copy(Array(S)))
            end
        end

        if testmode
            break
        end
    end
    toc = time()
    t_elapsed = toc - tic
    timed_iters = step - 3

    if do_vis
        save("/tmp/sim_results.jld", "T_storage", T_storage, "W_storage", W_storage, "S_storage", S_storage)
    end

    println("time, step: $(sim_time) $(step)")

    return SimOut_t(T, W, S, t_elapsed, timed_iters)
end

if abspath(PROGRAM_FILE) == @__FILE__
    opt = SimIn_t()
    opt.beta = 0.5
    opt.Pr = 1.0e-1
    opt.tol = 1.0e-7
    opt.nx = 1025
    opt.ny = 257
    navier_stokes_2D(opt=opt, do_vis=false)
end


text:
OPTIMIZER = nothing

"""
Set optimizer for all quadratic and linear programming subproblems.
"""
function qap_set_optimizer!(o)
	global OPTIMIZER = o
end

function get_optimizer()
	if isnothing(OPTIMIZER)
		error("Need to call set_optimizer! before using this solver.")
	else
		return OPTIMIZER
	end
end


text:
##
using Revise
using GLMakie
import Pkg
Pkg.activate(".")
##
using Images
using FileIO
using StaticArrays
import Statistics
using NearestNeighbors
using ProgressMeter
##
function preprocess_image(img_path)
    img = load(img_path)
    return @. Float32(Images.Gray(img))
end
##
desmin, actin = map(preprocess_image, ["./Data/SR2_desmin_image.jpg", "./Data/SR3_actin_image.jpg"])

function pointcloud_from_image(image::AbstractArray{T, N}) where {T, N}
    coords = map(findall(x -> x > zero(T), image)) do ind
        ind |> Tuple |> SVector{N, T}
    end
    return coords
end
##
import ImageBinarization: binarize, Otsu
algo = Otsu()
threshold(img) = binarize(img, algo)
##
let
    fig = Figure(resolution = (1200, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Desmin")
    ax2 = GLMakie.Axis(fig[1, 2], title = "Thresholded (Otsu)")
    image!(ax1, desmin, colormap = :Reds_9)
    image!(ax2, threshold(desmin))
    fig
end
##
let
    fig = Figure(resolution = (1200, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Actin")
    ax2 = GLMakie.Axis(fig[1, 2], title = "Thresholded (Otsu)")
    image!(ax1, actin, colormap = :Greens_9)
    image!(ax2, threshold(actin))
    fig
end
##
let
    xrange, yrange = 1900:2900, 1000:3000
    fig = Figure(resolution = (1500, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Desmin (Thresholded)")
    ax2 = GLMakie.Axis(fig[1, 2], title = "Actin (Thresholded)")
    image!(ax1, threshold(desmin)[xrange, yrange])
    image!(ax2, threshold(actin)[xrange, yrange])
    fig
end

##
xrange, yrange = 1900:2900, 1000:3000
desmin_coord = threshold(desmin)[xrange, yrange] |> pointcloud_from_image;
actin_coord = threshold(actin)[xrange, yrange] |> pointcloud_from_image;
##
size(actin_coord)
##
import ThreadTools: tmap1
function RipleyK(r::T, tree::NNTree, area) where {T}
    counts = Statistics.mean([length(inrange(tree, p, r)) for p in tree.data])
    N = length(tree.data)
    return counts*area/N
end

function RipleyK(r, T, area, p::ProgressMeter.AbstractProgress)
    K = RipleyK(r, T, area)
    next!(p)
    return K
end

function RipleyK(r::AbstractVector, tree::NNTree, area)
    p = Progress(length(r))
    return tmap1(x -> RipleyK(x, tree, area, p), r)
end

RipleyK(r, points, area) = RipleyK(r, KDTree(points), area)

##
r = 0:5:100
trees = map(KDTree, [desmin_coord, actin_coord])
##
K_desmin = RipleyK(r, trees[1], 1000*2000)
##
K_actin = RipleyK(r, trees[2], 1000*2000)
##
H(r, K) = sqrt.(K/π) - r
##
H_desmin, H_actin = H(r, K_desmin), H(r, K_actin)
##
let
    fig = Figure(resolution = (1200, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Ripley K function")
    l1 = lines!(r, K_desmin)
    l2 = lines!(r, K_actin)
    Legend(fig[1, 2], [l1, l2], ["Desmin", "Actin"])

    ax1 = GLMakie.Axis(fig[1, 3], title = "Ripley H function")
    l3 = lines!(r, H_desmin)
    l4 = lines!(r, H_actin)
    Legend(fig[1, 4], [l3, l4], ["Desmin", "Actin"])

    fig
end
##
lines(r, H)
##

relative_angle(coord1::SVector{2, T}, coord2::SVector{2, T}) where {T} = atan(coord2[2]-coord1[2], coord2[1]-coord1[1])

function get_inrange_angles(r::T, tree, coord::SVector{2, T}) where {T}
    idx = inrange(tree, coord, r)
    angles = [relative_angle(coord, tree.data[i]) for i in idx]
    return angles
end

import StatsBase: Histogram, fit
function get_angle_bin_counts(angles, angles_bins)
    h = fit(Histogram, angles, angle_bins, closed = :left)
    return h.weights
end

r = 20.0f0
Δα = 5π/180
A = 1000*2000
angle_bins = float32.(collect(-π:Δα:π))
##
angle_counts_desmin = (2π/Δα)*(A/length(desmin_coord)^2)*sum(tmap1(c -> get_angle_bin_counts(get_inrange_angles(r, trees[1], c), angle_bins), desmin_coord))
angle_counts_actin = (2π/Δα)*(A/length(actin_coord)^2)*sum(tmap1(c -> get_angle_bin_counts(get_inrange_angles(r, trees[2], c), angle_bins), actin_coord))

##
extrema(angle_counts_desmin)
##

##
angle_bins
##
let
    fig = Figure()
    xticks = ([-π, -π/2, 0,  π/2, π], ["-π", "-π/2", "0",  "π/2", "π"])
    ax = GLMakie.Axis(fig[1, 1], xticks = xticks, xlabel = "Angle (rad)", title = "Angular Ripley's function K(α; r = $(r))")
    A = 1000*2000
    l1 = lines!(ax, angle_bins[1:end-1], sqrt.(angle_counts_desmin/π).-r)
    l2 = lines!(ax, angle_bins[1:end-1], sqrt.(angle_counts_actin/π).-r)

    Legend(fig[1, 2], [l1, l2], ["Desmin", "Actin"])
    fig
end
##


text:
include(joinpath("dataAnalysis", "Occurrence.jl"))
include(joinpath("dataAnalysis", "Duration.jl"))
include(joinpath("dataAnalysis", "testFunctions.jl"))

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
# src/fiat.jl
# Fiat Currencies

export f

@reexport module f

using Unitful
using .._unitfulcurrency

export USD, EUR, GBP, CAD, DKK, JPY, CNY, KRW

@currencyunit AED UAEDirham
@currencyunit AFN Afghani
@currencyunit ALL Lek
@currencyunit AMD ArmenianDram
@currencyunit ANG NetherlandsAntilleanGuilder
@currencyunit AOA Kwanza
@currencyunit ARS ArgentinePeso
@currencyunit AUD AustralianDollar
@currencyunit AWG ArubanFlorin
@currencyunit AZN AzerbaijanManat
@currencyunit BAM ConvertibleMark
@currencyunit BBD BarbadosDollar
@currencyunit BDT Taka
@currencyunit BGN BulgarianLev
@currencyunit BHD BahrainiDinar
@currencyunit BIF BurundiFranc
@currencyunit BMD BermudianDollar
@currencyunit BND BruneiDollar
@currencyunit BOB Boliviano
@currencyunit BRL BrazilianReal
@currencyunit BSD BahamianDollar
@currencyunit BTN Ngultrum
@currencyunit BWP Pula
@currencyunit BYN BelarusianRuble
@currencyunit BZD BelizeDollar
@currencyunit CAD CanadianDollar
@currencyunit CDF CongoleseFranc
@currencyunit CHF SwissFranc
@currencyunit CLP ChileanPeso
@currencyunit CNY YuanRenminbi
@currencyunit COP ColombianPeso
@currencyunit CRC CostaRicanColon
@currencyunit CUC PesoConvertible
@currencyunit CUP CubanPeso
@currencyunit CVE CaboVerdeEscudo
@currencyunit CZK CzechKoruna
@currencyunit DJF DjiboutiFranc
@currencyunit DKK DanishKrone
@currencyunit DOP DominicanPeso
@currencyunit DZD AlgerianDinar
@currencyunit EGP EgyptianPound
@currencyunit ERN Nakfa
@currencyunit ETB EthiopianBirr
@currencyunit EUR Euro
@currencyunit FJD FijiDollar
@currencyunit GBP PoundSterling
@currencyunit GEL Lari
@currencyunit GHS GhanaCedi
@currencyunit GIP GibraltarPound
@currencyunit GMD Dalasi
@currencyunit GNF GuineanFranc
@currencyunit GTQ Quetzal
@currencyunit GYD GuyanaDollar
@currencyunit HKD HongKongDollar
@currencyunit HNL Lempira
@currencyunit HRK Kuna
@currencyunit HTGUSD GourdeUSDollar
@currencyunit HUF Forint
@currencyunit IDR Rupiah
@currencyunit ILS NewIsraeliSheqel
@currencyunit INR IndianRupee
@currencyunit IQD IraqiDinar
@currencyunit IRR IranianRial
@currencyunit ISK IcelandKrona
@currencyunit JMD JamaicanDollar
@currencyunit JOD JordanianDinar
@currencyunit JPY Yen
@currencyunit KES KenyanShilling
@currencyunit KGS Som
@currencyunit KHR Riel
@currencyunit KMF ComorianFranc
@currencyunit KPW NorthKoreanWon
@currencyunit KRW Won
@currencyunit KWD KuwaitiDinar
@currencyunit KZT Tenge
@currencyunit LAK LaoKip
@currencyunit LBP LebanesePound
@currencyunit LKR SriLankaRupee
@currencyunit LRD LiberianDollar
@currencyunit LSLZAR LotiRand
@currencyunit LYD LibyanDinar
@currencyunit MAD MoroccanDirham
@currencyunit MDL MoldovanLeu
@currencyunit MGA MalagasyAriary
@currencyunit MKD Denar
@currencyunit MMK Kyat
@currencyunit MNT Tugrik
@currencyunit MOP Pataca
@currencyunit MRO Ouguiya
@currencyunit MUR MauritiusRupee
@currencyunit MVR Rufiyaa
@currencyunit MWK MalawiKwacha
@currencyunit MXN MexicanPeso
@currencyunit MYR MalaysianRinggit
@currencyunit MZN MozambiqueMetical
@currencyunit NADZAR NamibiaDollarRand
@currencyunit NGN Naira
@currencyunit NIO CordobaOro
@currencyunit NOK NorwegianKrone
@currencyunit NPR NepaleseRupee
@currencyunit NZD NewZealandDollar
@currencyunit OMR RialOmani
@currencyunit PABUSD BalboaUSDollar
@currencyunit PEN Sol
@currencyunit PGK Kina
@currencyunit PHP PhilippinePiso
@currencyunit PKR PakistanRupee
@currencyunit PLN Zloty
@currencyunit PYG Guarani
@currencyunit QAR QatariRial
@currencyunit RON RomanianLeu
@currencyunit RSD SerbianDinar
@currencyunit RUB RussianRuble
@currencyunit RWF RwandaFranc
@currencyunit SAR SaudiRiyal
@currencyunit SBD SolomonIslandsDollar
@currencyunit SCR SeychellesRupee
@currencyunit SDG SudanesePound
@currencyunit SEK SwedishKrona
@currencyunit SGD SingaporeDollar
@currencyunit SHP SaintHelenaPound
@currencyunit SLL Leone
@currencyunit SOS SomaliShilling
@currencyunit SRD SurinamDollar
@currencyunit SSP SouthSudanesePound
@currencyunit STD Dobra
@currencyunit SVC ElSalvadorColon
@currencyunit SYP SyrianPound
@currencyunit SZL Lilangeni
@currencyunit THB Baht
@currencyunit TJS Somoni
@currencyunit TMT TurkmenistanNewManat
@currencyunit TND TunisianDinar
@currencyunit TOP Paanga  # Pa\u2019anga
@currencyunit TRY TurkishLira
@currencyunit TTD TrinidadandTobagoDollar
@currencyunit TZS TanzanianShilling
@currencyunit UAH Hryvnia
@currencyunit UGX UgandaShilling
@currencyunit USD USDollar
@currencyunit UYU PesoUruguayo
@currencyunit UZS UzbekistanSum
@currencyunit VEF Bolivar  # Bol\u00edvar
@currencyunit VND Dong
@currencyunit VUV Vatu
@currencyunit WST Tala
@currencyunit XAF CFAFrancBEAC
@currencyunit XCD EastCaribbeanDollar
@currencyunit XOF CFAFrancBCEAO
@currencyunit XPF CFPFranc
@currencyunit XYD CaymanIslandsDollar
@currencyunit YER YemeniRial
@currencyunit ZAR Rand
@currencyunit ZMW ZambianKwacha
@currencyunit ZWL ZimbabweDollar

end  # module f


text:
"""
    ICVI-Exmamples.jl

Description:
    Placeholder for the package file for the sake of testing because this is a project and not a package.

Author:
    Sasha Petrenko <sap625@mst.edu>
"""


text:
module _vscodeserver

using REPL, Sockets, Base64, Pkg, UUIDs
import Base: display, redisplay
import Dates

include("../languageserver/packages/JSON/src/JSON.jl")

include("gridviewer.jl")

module JSONRPC
    import ..JSON
    import ..UUIDs

    include("../packages/JSONRPC/src/core.jl")
end

include("repl.jl")
include("../debugger/debugger.jl")

struct InlineDisplay <: AbstractDisplay end

repl_pipename = Base.ARGS[1]

!(Sys.isunix() || Sys.iswindows()) && error("Unknown operating system.")

function ends_with_semicolon(x)
    return REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])
end

repl_conn = connect(repl_pipename)

conn_endpoint = JSONRPC.JSONRPCEndpoint(repl_conn, repl_conn)

function sendDisplayMsg(kind, data)
    JSONRPC.send_notification(conn_endpoint, "display", Dict{String,String}("kind"=>kind, "data"=>data))
end

run(conn_endpoint)

@async begin

    while true
        msg = JSONRPC.get_next_message(conn_endpoint)

        if msg["method"] == "repl/runcode"
            params = msg["params"]

            source_filename = params["filename"]
            code_line = params["line"]
            code_column = params["column"]
            source_code = params["code"]

            hideprompt() do
                if isdefined(Main, :Revise) && isdefined(Main.Revise, :revise) && Main.Revise.revise isa Function
                    let mode = get(ENV, "JULIA_REVISE", "auto")
                        mode == "auto" && Main.Revise.revise()
                    end
                end
                for (i,line) in enumerate(eachline(IOBuffer(source_code)))
                    if i==1
                        printstyled("julia> ", color=:green)
                        print(' '^code_column)
                    else
                        # Indent by 7 so that it aligns with the julia> prompt
                        print(' '^7)
                    end

                    println(line)
                end

                try
                    withpath(source_filename) do
                        res = Base.invokelatest(include_string, Main, '\n'^code_line * ' '^code_column *  source_code, source_filename)

                        if res !== nothing && !ends_with_semicolon(source_code)
                            Base.invokelatest(display, res)
                        end
                    end
                catch err
                    Base.display_error(stderr, err, catch_backtrace())
                end
            end
        elseif msg["method"] == "repl/startdebugger"
            hideprompt() do
                debug_pipename = msg["params"]
                try
                    VSCodeDebugger.startdebug(debug_pipename)
                catch err
                    VSCodeDebugger.global_err_handler(err, catch_backtrace(), ARGS[4])
                end
            end
        end
    end
end

function display(d::InlineDisplay, ::MIME{Symbol("image/png")}, x)
    payload = stringmime(MIME("image/png"), x)
    sendDisplayMsg("image/png", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/png")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)
    sendDisplayMsg("image/svg+xml", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("text/html")}, x)
    payload = stringmime(MIME("text/html"), x)
    sendDisplayMsg("text/html", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("text/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}, x)
    payload = stringmime(MIME("juliavscode/html"), x)
    sendDisplayMsg("juliavscode/html", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("juliavscode/html")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v2+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v3+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v4+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v3+json"), x)
    sendDisplayMsg("application/vnd.vega.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v4+json"), x)
    sendDisplayMsg("application/vnd.vega.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v5+json"), x)
    sendDisplayMsg("application/vnd.vega.v5+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}, x)
    payload = stringmime(MIME("application/vnd.plotly.v1+json"), x)
    sendDisplayMsg("application/vnd.plotly.v1+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}, x)
    payload = stringmime(MIME("application/vnd.dataresource+json"), x)
    sendDisplayMsg("application/vnd.dataresource+json", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("application/vnd.dataresource+json")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}) = true

function display(d::InlineDisplay, x)
    if showable("application/vnd.vegalite.v4+json", x)
        display(d,"application/vnd.vegalite.v4+json", x)
    elseif showable("application/vnd.vegalite.v3+json", x)
        display(d,"application/vnd.vegalite.v3+json", x)
    elseif showable("application/vnd.vegalite.v2+json", x)
        display(d,"application/vnd.vegalite.v2+json", x)
    elseif showable("application/vnd.vega.v5+json", x)
        display(d,"application/vnd.vega.v5+json", x)
    elseif showable("application/vnd.vega.v4+json", x)
        display(d,"application/vnd.vega.v4+json", x)
    elseif showable("application/vnd.vega.v3+json", x)
        display(d,"application/vnd.vega.v3+json", x)
    elseif showable("application/vnd.plotly.v1+json", x)
        display(d,"application/vnd.plotly.v1+json", x)
    elseif showable("juliavscode/html", x)
        display(d,"juliavscode/html", x)
    # elseif showable("text/html", x)
    #     display(d,"text/html", x)
    elseif showable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif showable("image/png", x)
        display(d,"image/png", x)
    else
        throw(MethodError(display,(d,x)))
    end
end

function _display(d::InlineDisplay, x)
    if showable("application/vnd.dataresource+json", x)
        display(d, "application/vnd.dataresource+json", x)
    else
        try
            display(d, x)
        catch err
            if err isa MethodError
                @warn "Cannot display values of type $(typeof(x)) in VS Code."
            else
                rethrow(err)
            end
        end
    end
end

if length(Base.ARGS) >= 3 && Base.ARGS[3] == "true"
    atreplinit(i->Base.Multimedia.pushdisplay(InlineDisplay()))
end

# Load revise?
load_revise = Base.ARGS[2] == "true"

const tabletraits_uuid = UUIDs.UUID("3783bdb8-4a98-5b6b-af9a-565f29a5fe9c")
const datavalues_uuid = UUIDs.UUID("e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5")

global _isiterabletable = i -> false
global _getiterator = i -> i

function pkgload(pkg)
    if pkg.uuid==tabletraits_uuid
        x = Base.require(pkg)

        global _isiterabletable = x.isiterabletable
        global _getiterator = x.getiterator
    elseif pkg.uuid==datavalues_uuid
        x = Base.require(pkg)

        eval(quote
            function JSON_print_escaped(io, val::$(x.DataValue))
                $(x.isna)(val) ? print(io, "null") : JSON_print_escaped(io, val[])
            end

            julia_type_to_schema_type(::Type{T}) where {S, T<:$(x.DataValue){S}} = julia_type_to_schema_type(S)
        end)
    end
end

push!(Base.package_callbacks, pkgload)

function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

end

function vscodedisplay(x)
    if showable("application/vnd.dataresource+json", x)
        _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
    elseif _vscodeserver._isiterabletable(x)===true
        buffer = IOBuffer()
        io = IOContext(buffer, :compact=>true)
        _vscodeserver.printdataresource(io, _vscodeserver._getiterator(x))
        buffer_asstring = _vscodeserver.CachedDataResourceString(String(take!(buffer)))
        _vscodeserver._display(_vscodeserver.InlineDisplay(), buffer_asstring)
    elseif _vscodeserver._isiterabletable(x)===missing
        try
            buffer = IOBuffer()
            io = IOContext(buffer, :compact=>true)
            _vscodeserver.printdataresource(io, _vscodeserver._getiterator(x))
            buffer_asstring = _vscodeserver.CachedDataResourceString(String(take!(buffer)))
            _vscodeserver._display(_vscodeserver.InlineDisplay(), buffer_asstring)
        catch err
            _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
        end
    elseif x isa AbstractVector || x isa AbstractMatrix
        buffer = IOBuffer()
        io = IOContext(buffer, :compact=>true)
        _vscodeserver.print_array_as_dataresource(io, _vscodeserver._getiterator(x))
        buffer_asstring = _vscodeserver.CachedDataResourceString(String(take!(buffer)))
        _vscodeserver._display(_vscodeserver.InlineDisplay(), buffer_asstring)
    else
        _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
    end
end

vscodedisplay() = i -> vscodedisplay(i)

if _vscodeserver.load_revise
    try
        @eval using Revise
        Revise.async_steal_repl_backend()
    catch err
    end
end

macro enter(command)
    _vscodeserver.remove_lln!(command)
    :(_vscodeserver.JSONRPC.send_notification(_vscodeserver.conn_endpoint, "debugger/enter", $(string(command))))
end

macro run(command)
    _vscodeserver.remove_lln!(command)
    :(_vscodeserver.JSONRPC.send_notification(_vscodeserver.conn_endpoint, "debugger/run", $(string(command))))
end

max_stars_repo_name: cloudgiant/julia-vscode

max_issues_repo_name: cloudgiant/julia-vscode

max_forks_repo_name: cloudgiant/julia-vscode


text:
using ExactWrightFisher, Distributions

@testset "test full inference CIR" begin
    Random.seed!(0)
    times_sim = range(0, stop = 20, length = 20)
    X = DualOptimalFiltering.generate_CIR_trajectory(times_sim, 3, 3., 0.5, 1);
    Y = map(λ -> rand(Poisson(λ),1), X);
    data = Dict(zip(times_sim, Y))
    δ, γ, σ, λ = 3., 0.5, 1., 1.
    Λ_of_t, logwms_of_t, θ_of_t, = DualOptimalFiltering.filter_CIR_logweights(δ, γ, σ, λ, data);

    times = Λ_of_t |> keys |> collect |> sort

    @test_nowarn DualOptimalFiltering.hidden_signal_loglikelihood_CIR(X, times, δ, γ, σ)

    @test_nowarn DualOptimalFiltering.emitted_data_conditional_loglikelihood_CIR(data, X, times, λ)

    @test_nowarn DualOptimalFiltering.joint_loglikelihood_CIR(data, X, times, δ, γ, σ, λ)

    prior_δ = truncated(Normal(5, 4), 0, Inf)
    prior_γ = truncated(Normal(5, 4), 0, Inf)
    prior_σ = truncated(Normal(5, 4), 0, Inf)
    prior_logpdf(δi, γi, σi) = logpdf(prior_δ, δi) + logpdf(prior_γ, γi) + logpdf(prior_σ, σi)

    @test_nowarn DualOptimalFiltering.joint_sampler_CIR_pruning_precompute(data, λ, prior_logpdf, [1.,1.,1.], 10, (x, y) -> (x, y), final_chain_length = 10)

    @test_nowarn DualOptimalFiltering.joint_sampler_CIR_keep_fixed_number_precompute(data, λ, prior_logpdf, [1.,1.,1.], 10, 2; final_chain_length = 10)

    @test_nowarn DualOptimalFiltering.joint_sampler_CIR_reparam_keep_fixed_number_precompute(data, λ, prior_logpdf, 10, 2; final_chain_length = 10)

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# functions for convenience of use

## diagonal aview

function diagview(a::DenseArray{T,2}) where T
    m, n = size(a)
    s1, s2 = astrides(a)
    len = min(m, n)
    StridedView(parent(a), offset(a), (len,), ContRank{0}, (s1 + s2,))
end

## row vector aview

function rowvec_view(a::DenseArray{T,2}, i::Integer) where T
    m, n = size(a)
    s1, s2 = astrides(a)
    StridedView(parent(a), offset(a) + (i-1) * s1, (n,), ContRank{0}, (s2,))
end

## flatten_view

flatten_view(a::ContiguousArray) =
    ContiguousView(parent(a), offset(a), (length(a),))

## reshape_view

function reshape_view(a::ContiguousArray, shp::NTuple{N,Int}) where N
    prod(shp) == length(a) || throw(DimensionMismatch("Inconsistent array size."))
    ContiguousView(parent(a), offset(a), shp)
end

## ellipview

ellipview(a::DenseArray{T,2}, i::Union{Integer, UnitRange}) where {T} = aview(a, :, i)
ellipview(a::DenseArray{T,3}, i::Union{Integer, UnitRange}) where {T} = aview(a, :, :, i)
ellipview(a::DenseArray{T,4}, i::Union{Integer, UnitRange}) where {T} = aview(a, :, :, :, i)
ellipview(a::DenseArray{T,5}, i::Union{Integer, UnitRange}) where {T} = aview(a, :, :, :, :, i)
ellipview(a::DenseArray{T,N}, i::Union{Integer, UnitRange}) where {T,N} = aview(a, ntuple(i->Colon(),N-1)..., i)