text:
# Tree implementations, including TokenTree

"""
A tree having nodes of type NodeType.

Any subtype must implement the tree operations defined in this source file.
"""
abstract type Tree{NodeType}
end

"""
Tree implementation using indices into a node vector.

Instance vields are considered private, they must not be changed
outside the tree API, or tree corruption could happen.

idx stores the tree structure with the following conventions:

idx[n] == 0: marks the end of a subtree, no valid node index

idx[n]>0: base[idx[n]] is a leaf node in the tree

idx[n]<0: base[-idx[n]] is a node of a subtree having 0 or more children
idx[n+1]..idx[n+k] are the subtrees of all children (k=0 if no children present)
idx[n+k+1] is 0, marking the end if the list of children

"""
struct ITree{I<:Integer, NodeType} <: Tree{NodeType}
    idx ::Vector{I}
    base :: Vector{NodeType}
end

"""
A tree variant of TokenVector.

A TokenTree is a TokenVector organized as a tree.
In contrast to  but it offers an additional tree API.

A tree node is identified by an index into the token array.
Depending on its category, it is a leaf node (category <= T_END),
 or it has a (probably empty) subtree (category > T_END): its children are
the following nodes in the token array, until termination by a T_END token.
T_END tokens are no tree nodes not part of the tree structure, itThe index of a T_END token

"""
mutable struct TokenTree <: AbstractTokenVector{HybridFly}
    vec :: Vector{HybridFly}
    buffer :: String # private storage shared by all elements
    used :: UInt32 # no of bytes currently used
    shared :: UInt32 # no of bytes in buffer shared with other objects
    flags :: Int
end

# alternative: ein tree ist ein int32 array. Wert >0: leaf node.
# wert <0: subtree
# wert =0: end of subtree
#
# Vorteil: kann über JEDEM Vector gebaut werden (!!)
# vorteil: halbe größe der END marker.
# nachteil: zusätzliches array + 1 indirektion mehr


text:
"""
Liu S., Yang Y., Forrest J. (2017) Grey Numbers and Their Operations. In: Grey Data Analysis.
Computational Risk Management. Springer, Singapore. https://doi.org/10.1007/978-981-10-1841-1_3
"""
module GreyNumbers

struct GreyNumber{T<:Real}
    a::T
    b::T
    GreyNumber{T}(x::T, y::T) where {T <: Real} = new{T}(min(x, y), max(x, y))
end

GreyNumber(x::T, y::T) where {T<:Real} = GreyNumber{T}(x,y)

GreyNumber()::GreyNumber{Float64} = GreyNumber(0.0, 0.0)

GreyNumber(x::T) where {T <: Real} = GreyNumber{T}(x, x)

GreyNumber{T}(x::GreyNumber{T}) where {T <: Real} = GreyNumber{T}(x.a, x.b)

GreyNumber{T}(x::Bool)  where {T <: Real}= GreyNumber(Float64(x), Float64(x))

Base.transpose(x::GreyNumber{T}) where {T <: Real} = x

Base.adjoint(x::GreyNumber{Float64}) where {T <: Real} = x

Base.copy(x::GreyNumber{T}) where {T <: Real} = GreyNumber{T}(x.a, x.b)

Base.length(g::GreyNumber) = 1

function whitenizate(g::GreyNumber; t::Float64 = 0.5)::Float64
    @assert 0.0 <= t <= 1.0
    return g.a * t + g.b * (1.0 - t)
end

function Base.:+(g1::GreyNumber, g2::GreyNumber)::GreyNumber
    return GreyNumber(g1.a + g2.a, g1.b + g2.b)
end

function Base.:+(g::GreyNumber{T}, k::T)::GreyNumber{T} where T <: Real
    return GreyNumber(g.a + k, g.b + k)
end

function Base.:+(k::T, g::GreyNumber{T})::GreyNumber{T} where T <: Real
    return g + k
end

function Base.:-(g::GreyNumber)::GreyNumber
    minusone = -1 * one(g.a)
    return GreyNumber(minusone * g.b, minusone * g.a)
end

function Base.:-(g1::GreyNumber, g2::GreyNumber)::GreyNumber
    return GreyNumber(g1.a - g2.b, g1.b - g2.a)
end

function Base.:-(g::GreyNumber{T}, k::T)::GreyNumber{T} where T <: Real
    return g - GreyNumber(k, k)
end

function Base.:-(k::T, g::GreyNumber{T})::GreyNumber{T} where T <: Real
    return GreyNumber(k, k) - g
end

function Base.:*(g1::GreyNumber, g2::GreyNumber)::GreyNumber
    v = [g1.a * g2.a, g1.a * g2.b, g1.b * g2.a, g1.b * g2.b]
    return GreyNumber(
        minimum(v),
        maximum(v),
    )
end

function Base.:*(k::T where {T <: Number}, g1::GreyNumber)::GreyNumber
    a = g1.a * k
    b = g1.b * k
    return GreyNumber(min(a, b), max(a,b))
end

function Base.:*(g1::GreyNumber, k::T where {T <: Number})::GreyNumber
    return k * g1
end

function Base.inv(g::GreyNumber)::GreyNumber
    newa = min(inv(g.a), inv(g.b))
    newb = max(inv(g.a), inv(g.b))
    return GreyNumber(newa, newb)
end

function Base.:/(g1::GreyNumber, g2::GreyNumber)::GreyNumber
    #@assert g2.a != 0
    #@assert g2.b != 0
    #@assert g2.a * g2.b > 0
    v = [g1.a / g2.a, g1.a / g2.b, g1.b / g2.a, g1.b / g2.b]
    return GreyNumber(
        minimum(v),
        maximum(v)
    )
end

function Base.:/(k::T where {T <: Real}, g1::GreyNumber)::GreyNumber
    return g1 / GreyNumber(k, k)
end

function Base.:/(g1::GreyNumber, k::T where {T <: Real})::GreyNumber
    return g1 / GreyNumber(k, k)
end

function Base.:^(g::GreyNumber, k::T where {T <: Real})::GreyNumber
    @assert k > 0
    @assert g.a <= g.b
    anew = g.a^k
    bnew = g.b^k
    return GreyNumber(min(anew, bnew), max(anew, bnew))
end

function Base.isless(g1::GreyNumber, g2::GreyNumber)::Bool
    if g1.a < g2.a
        return true
    elseif g1.a == g2.a
        return g1.b < g2.b
    else
        return false
    end
end

function Base.:≤(g1::GreyNumber, g2::GreyNumber)::Bool
    if g1.a <= g2.a
        return true
    elseif g1.a == g2.a
        return g1.b <= g2.b
    else
        return false
    end
end

function Base.isgreater(g1::GreyNumber, g2::GreyNumber)::Bool
    if g1.a > g2.a
        return true
    elseif g1.a == g2.a
        return g1.b > g2.b
    else
        return false
    end
end

function Base.:≥(g1::GreyNumber, g2::GreyNumber)::Bool
    if g1.a >= g2.a
        return true
    elseif g1.a == g2.a
        return g1.b >= g2.b
    else
        return false
    end
end

function Base.isequal(g1::GreyNumber, g2::GreyNumber)::Bool
    return (g1.a == g2.a) && (g1.b == g2.b)
end

function Base.isreal(g::GreyNumber)::Bool
    return isreal(g.a)
end

function Base.isinteger(g::GreyNumber)::Bool
    return isinteger(g.a)
end

function Base.isinf(g::GreyNumber)::Bool
    return isinf(g.a) || isinf(g.b)
end

function Base.iszero(g::GreyNumber)::Bool
    return iszero(g.a) && iszero(g.b)
end

function Base.isone(g::GreyNumber)::Bool
    return isone(g.a) && isone(g.b)
end

function Base.iseven(g::GreyNumber)::Bool
    @assert eltype(g.a) <: Int
    return iseven(g.a) && iseven(g.b)
end

function Base.isodd(g::GreyNumber)::Bool
    @assert eltype(g.a) <: Int
    return isodd(g.a) && isodd(g.b)
end

function Base.sum(g::GreyNumber{T})::T where {T <: Number}
    return g.a + g.b
end

function Base.eltype(g::GreyNumber)::Type
    return eltype(g.a)
end

function Base.iterate(g::GreyNumber, state=1)
    if state == 1
        (g.a, state + 1)
    elseif  state == 2
        (g.b, state + 1)
    else
        nothing
    end
end

function Base.sqrt(g::GreyNumber)::GreyNumber
    return g^0.5
end

function Base.cbrt(g::GreyNumber)::GreyNumber
    return GreyNumber(cbrt(g.a), cbrt(g.b))
end

function Base.abs(g::GreyNumber)::GreyNumber
    return GreyNumber(abs(g.a), abs(g.b))
end

function Base.getindex(g::GreyNumber{T}, index::Int64)::T where T
    @assert index == 1 || index == 2
    if index == 1
        g.a
    else
        g.b
    end
end

function Base.last(g::GreyNumber{T})::T where T
    return g.b
end

function kernel(g::GreyNumber)::Float64
    return (g.a + g.b) / 2.0
end

end # end of module GreyNumbers

ext: jl


text:
using LinearAlgebra
include("Pauli.jl")
include("BlochRotation.jl")

Hadammard = [ 1  1 ; 1 -1 ];
CNOT12= [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0];


text:
###############################################
# FUNCTIONS TO LOAD MODEL IMPLEMENTATION CODE #
###############################################

## HELPERS

function _append!(program, ex, doprint::Bool, tick_early::Bool)
    str = string(ex)
    doprint && push!(program.args, :(print($str)))
    doprint && !tick_early && push!(program.args, :(println(" \u2714")))
    push!(program.args, ex)
    tick_early && doprint && push!(program.args, :(println(" \u2714")))
    return program
end

function _import(modl, api_pkg, pkg, doprint)
    # can be removed once MLJModel #331 is resolved:
    if pkg == :NearestNeighbors
        doprint && print("import NearestNeighbors")
        try
            modl.eval(:(import MLJModels))
        catch
            try
                modl.eval(:(import MLJ.MLJModels))
            catch
                error("Problem putting MLJModels into scope. ")
            end
        end
        modl.eval(:(import NearestNeighbors))
        doprint && println(" \u2714")
    else
        doprint && print("import $api_pkg")
        modl.eval(:(import $api_pkg))
        doprint && println(" \u2714")
    end
end

function _eval(modl, path::Union{Expr,Symbol})
    modl.eval(path)
end

## OVERLOADING load_path

"""
    load_path(model::String, pkg=nothing)

Return the load path for model type with name `model`, specifying the
package name `pkg` to resolve name conflicts if necessary.

    load_path(model)

Return the load path of a `model` instance or type. Usually requires
necessary model code to have been separately loaded. Supply a string
as above if code is not loaded.

"""
function MLJModelInterface.load_path(proxy::ModelProxy)
    handle = (name=proxy.name, pkg=proxy.package_name)
    return INFO_GIVEN_HANDLE[handle][:load_path]
end
function MLJModelInterface.load_path(name::String; pkg=nothing)
    proxy = info(name; pkg=pkg)
    return load_path(proxy)
end

# to also get pkg, which could be different from glue code pkg
# appearing in load_path:
function load_path_and_pkg(name::String; pkg=nothing, interactive=false)
    proxy = info(name; pkg=pkg, interactive=interactive)
    handle = (name=proxy.name, pkg=proxy.package_name)
    _info = INFO_GIVEN_HANDLE[handle]

    return _info[:load_path], _info[:package_name]
end

## THE CODE LOADING MACROS

"""
    @load ModelName pkg=nothing verbosity=0 add=false

Import the model type the model named in the first argument into the
calling module, specfying `pkg` in the case of an ambiguous name (to
packages providing a model type with the same name). Returns the model
type.

**Warning** In older versions of MLJ/MLJModels, `@load` returned an
*instance* instead.

To automatically add required interface packages to the current
environment, specify `add=true`. For interactive loading, use
`@iload` instead.

### Examples

    Tree = @load DecisionTreeRegressor
    tree = Tree()
    tree2 = Tree(min_samples_split=6)

    SVM = @load SVC pkg=LIBSVM
    svm = SVM()

See also [`@iload`](@ref)

"""
macro load(name_ex, kw_exs...)
    _load(__module__, name_ex, kw_exs...)
end

"""
    @iload ModelName

Interactive alternative to @load. Provides user with an optioin to
install (add) the required interface package to the current
environment, and to choose the relevant model-providing package in
ambiguous cases.  See [`@load`](@ref)

"""
macro iload(name_ex, kw_exs...)
    _load(__module__, name_ex, kw_exs...; interactive=true)
end

# builds the program to be evaluated by the @load/@iload macros:
function _load(modl, name_ex, kw_exs...; interactive=false)

    # initialize:
    program = quote end

    # fallbacks:
    pkg_str = nothing
    verbosity = 1
    scope = :global
    install_pkgs = false

    # parse name_ex:
    name_ = string(name_ex)
    # get rid of parentheses in `name_`, as in
    # "(MLJModels.Clustering).KMedoids":
    name_str = filter(name_) do c !(c in ['(',')']) end
    name = Symbol(name_str)

    # parse kwargs:
    warning = "Invalid @load or @iload  syntax.\n "*
        "Sample usage: @load PCA pkg=MultivariateStats verbosity=0 install=true"
    for ex in kw_exs
        ex.head == :(=) || throw(ArgumentError(warning))
        variable_ex = ex.args[1]
        value_ex = ex.args[2]
        if variable_ex == :pkg
            pkg_str = string(value_ex)
        elseif variable_ex == :verbosity
            verbosity = value_ex
        elseif variable_ex in [:install_pkgs, :install, :add]
            install_pkgs = value_ex
        else
            throw(ArgumentError(warning))
        end
    end

    # are we printing stuff to stdout?
    doprint = verbosity > 0

    # next expression defines "program" variables `path`, `path_ex`,
    # `pkg_str`, `pkg`, `type_already_loaded`:
    ex  = quote
        $doprint && @info "For silent loading, specify `verbosity=0`. "
        proxy = MLJModels.info($name_str;
                               pkg=$pkg_str,
                               interactive=$interactive)
        handle = MLJModels.Handle(proxy.name, proxy.package_name)
        dic = MLJModels.INFO_GIVEN_HANDLE[handle]
        path_str, pkg_str = dic[:load_path], dic[:package_name]
        path = path_str |> Meta.parse
        pkg = Symbol(pkg_str)
    end
    push!(program.args, ex)

    ex = quote
        path_components = Symbol.(split(path_str, '.') )

        # get pkg containing implementation of model API implementation:
        api_pkg = path_components[1]
        api_pkg_str = string(api_pkg)

        if $install_pkgs || $interactive
            try
                MLJModels._import($modl, api_pkg, pkg, false)
            catch
                if $interactive
                    MLJModels.request(
                        "The package providing an interface "*
                        "to `$($name_str)` "*
                        "is not in your "*
                        "current environment.\n"*
                        "What do you want to do? ",
                        # choices:
                        "Install $api_pkg_str in current environment.",
                        "Abort.") == 1 || throw(InterruptException)
                end
                MLJModels._import($modl, :Pkg, :Pkg, false)
                Pkg.add(api_pkg_str)
            end
        end
        MLJModels._import($modl, api_pkg, pkg, $doprint)
        MLJModels._eval($modl, path)
    end
    push!(program.args, ex)

    return program

end

## NO LONGER SUPPORTED

_deperror() = error(
    "The `load` function is no longer supported. "*
    "Use the `@load` macro instead, as in "*
    "`@load RandomForestRegressor pkg = DecisionTree`.\n"*
    "For explicit importing, you can discover a model's "*
    "full load path with the `load_path` function, as in "*
    "`load_path(\"RandomForestRegressor\", pkg=\"DecisionTree\")`. )")

load(proxy::ModelProxy; kwargs...) = _deperror()
load(name::String; kwargs...) = _deperror()


text:
import DielectricUpdateTechique as DUTCH
using Test

f = 300f6; #300MHz
m = DUTCH.get_constants(f);

pif32 = π |> Float32;

@test m[:c]  == 299792458
@test m[:μ₀] == 4f-7*pif32
@test m[:ϵ₀] == 1/(m[:c]^2*m[:μ₀])
@test m[:f₀] == f
@test m[:ω]  == 2*pif32*f
@test m[:kb] == 2*pif32*f/m[:c]
@test m[:λ]  == m[:c]/f

f = 1337f6; #1337MHz
m = DUTCH.get_constants(f);
@test m[:f₀] == f
@test m[:ω]  == 2*pif32*f
@test m[:kb] == 2*pif32*f/m[:c]
@test m[:λ]  == m[:c]/f


text:
function print_indent(io, n)
    for d in 1:n
        print(io, "  ")
    end
end

function show(io::IO, cursor::RecordCursor)
    par = cursor.par
    rows = cursor.colcursors[1].rows
    println(io, "Record Cursor on $(par.path)")
    println(io, "    rows: $rows")

    colpaths = [join(colname, '.') for colname in cursor.colnames]
    println(io, "    cols: $(join(colpaths, ", "))")
end

function show(io::IO, cursor::BatchedColumnsCursor)
    par = cursor.par
    rows = cursor.colcursors[1].rows
    println(io, "Batched Columns Cursor on $(par.path)")
    println(io, "    rows: $rows")
    println(io, "    batches: $(length(cursor))")

    colpaths = [join(colname, '.') for colname in cursor.colnames]
    println(io, "    cols: $(join(colpaths, ", "))")
end

function show(io::IO, schema::SchemaElement, indent::AbstractString="", nchildren::Vector{Int}=Int[])
    print(io, indent)
    lchildren = length(nchildren)
    print_indent(io, lchildren)
    if hasproperty(schema, :repetition_type)
        r = schema.repetition_type
        print(io, (r == FieldRepetitionType.REQUIRED) ? "required" : (r == FieldRepetitionType.OPTIONAL) ? "optional" : "repeated", " ");
    end
    hasproperty(schema, :_type) && print(io, Thrift.enumstr(_Type, schema._type), " ")

    print(io, schema.name)
    hasproperty(schema, :field_id) && print(io, " (", schema.field_id, ")")

    if hasproperty(schema, :converted_type)
        print(io, " # (from ", Thrift.enumstr(ConvertedType, schema.converted_type))
        if schema.converted_type == ConvertedType.DECIMAL
            print(io, "(", schema.precision, ",", schema.scale, ")")
        end
        print(io, ") ")
    end

    if hasproperty(schema, :num_children) && (getproperty(schema, :num_children) > 0)
        push!(nchildren, schema.num_children)
        print(io, " {")
    elseif lchildren > 0
        nchildren[lchildren] -= 1
        if nchildren[lchildren] == 0
            pop!(nchildren)
            println(io, "")
            print_indent(io, length(nchildren))
            print(io, indent, "}")
        end
    end

    println(io, "")
end

function show(io::IO, schema::Vector{SchemaElement}, indent::AbstractString="")
    println(io, indent, "Schema:")
    nchildren=Int[]
    for schemaelem in schema
        show(io, schemaelem, indent * "    ", nchildren)
    end
end

show(io::IO, schema::Schema, indent::AbstractString="") = show(io, schema.schema, indent)

function show(io::IO, kvmeta::KeyValue, indent::AbstractString="")
    println(io, indent, kvmeta.key, " => ", kvmeta.value)
end

function show(io::IO, kvmetas::Vector{KeyValue}, indent::AbstractString="")
    isempty(kvmetas) && return
    println(io, indent, "Metadata:")
    for kvmeta in kvmetas
        show(io, kvmeta, indent * "  ")
    end
end

function show_encodings(io::IO, encodings::Vector{Int32}, indent::AbstractString="")
    isempty(encodings) && return
    print(io, indent, "Encodings: ")
    pfx = ""
    for encoding in encodings
        print(io, pfx, Thrift.enumstr(Encoding, encoding))
        pfx = ", "
    end
    println(io, "")
end

show(io::IO, hdr::IndexPageHeader, indent::AbstractString="") = nothing
function show(io::IO, page::DictionaryPageHeader, indent::AbstractString="")
    println(io, indent, page.num_values, " values")
end

function show(io::IO, hdr::DataPageHeader, indent::AbstractString="")
    println(io, indent, hdr.num_values, " values")
    println(io, indent, "encodings: values as ", Thrift.enumstr(Encoding, hdr.encoding), ", definitions as ", Thrift.enumstr(Encoding, hdr.definition_level_encoding), ", repetitions as ", Thrift.enumstr(Encoding, hdr.repetition_level_encoding))
    hasproperty(hdr, :statistics) && show(io, hdr.statistics, indent)
end

function show(io::IO, hdr::DataPageHeaderV2, indent::AbstractString="")
    compressed = hasproperty(hdr, :is_compressed) ? hdr.is_compressed : true
    println(io, indent, hdr.num_values, " values, ", hdr.num_nulls, " nulls, ", hdr.num_rows, " rows, compressed:", compressed)
    println(io, indent, "encoding:", Thrift.enumstr(Encoding, hdr.encoding), ", definition:", Thrift.enumstr(Encoding, hdr.definition_level_encoding), ", repetition:", Thrift.enumstr(Encoding, hdr.repetition_level_encoding))
    hasproperty(hdr, :statistics) && show(io, hdr.statistics, indent)
end

function show(io::IO, page::PageHeader, indent::AbstractString="")
    println(io, indent, Thrift.enumstr(PageType, page._type), " compressed bytes:", page.compressed_page_size, " (", page.uncompressed_page_size, " uncompressed)")
    hasproperty(page, :data_page_header) && show(io, page.data_page_header, indent * "  ")
    hasproperty(page, :data_page_header_v2) && show(io, page.data_page_header_v2, indent * "  ")
    hasproperty(page, :index_page_header) && show(io, page.index_page_header, indent * "  ")
    hasproperty(page, :dictionary_page_header) && show(io, page.dictionary_page_header, indent * "  ")
end

function show(io::IO, pages::Vector{PageHeader}, indent::AbstractString="")
    println(io, indent, "Pages:")
    for page in pages
        show(io, page, indent * "  ")
    end
end

show(io::IO, page::Page, indent::AbstractString="") = show(io, page.hdr, indent)
show(io::IO, pages::Vector{Page}, indent::AbstractString="") = show(io, [page.hdr for page in pages], indent)

function show(io::IO, stat::Statistics, indent::AbstractString="")
    println(io, indent, "Statistics:")
    if hasproperty(stat, :min) && hasproperty(stat, :max)
        println(io, indent, "  range:", stat.min, ":", stat.max)
    elseif hasproperty(stat, :min)
        println(io, indent, "  min:", stat.min)
    elseif hasproperty(stat, :max)
        println(io, indent, "  max:", stat.max)
    end
    hasproperty(stat, :null_count) && println(io, indent, "  null count:", stat.null_count)
    hasproperty(stat, :distinct_count) && println(io, indent, "  distinct count:", stat.distinct_count)
end

function show(io::IO, page_enc::PageEncodingStats, indent::AbstractString="")
    println(io, indent, page_enc.count, " ", Thrift.enumstr(Encoding, page_enc.encoding), " encoded ", Thrift.enumstr(PageType, page_enc.page_type), " pages")
end

function show(io::IO, page_encs::Vector{PageEncodingStats}, indent::AbstractString="")
    isempty(page_encs) && return
    println(io, indent, "Page encoding statistics:")
    for page_enc in page_encs
        show(io, page_enc, indent * "  ")
    end
end

function show(io::IO, colmeta::ColumnMetaData, indent::AbstractString="")
    println(io, indent, Thrift.enumstr(_Type, coltype(colmeta)), " ", join(colname(colmeta), '.'), ", num values:", colmeta.num_values)
    show_encodings(io, colmeta.encodings, indent)
    if colmeta.codec != CompressionCodec.UNCOMPRESSED
        println(io, indent, Thrift.enumstr(CompressionCodec, colmeta.codec), " compressed bytes:", colmeta.total_compressed_size, " (", colmeta.total_uncompressed_size, " uncompressed)")
    else
        println(io, indent, Thrift.enumstr(CompressionCodec, colmeta.codec), " bytes:", colmeta.total_compressed_size)
    end

    print(io, indent, "offsets: data:", colmeta.data_page_offset)
    hasproperty(colmeta, :index_page_offset) && print(io, ", index:", colmeta.index_page_offset)
    hasproperty(colmeta, :dictionary_page_offset) && print(io, ", dictionary:", colmeta.dictionary_page_offset)
    println(io, "")
    hasproperty(colmeta, :statistics) && show(io, colmeta.statistics, indent)
    hasproperty(colmeta, :encoding_stats) && show(io, colmeta.encoding_stats, indent)
    hasproperty(colmeta, :key_value_metadata) && show(io, colmeta.key_value_metadata, indent)
end

function show(io::IO, columns::Vector{ColumnChunk}, indent::AbstractString="")
    for col in columns
        path = hasproperty(col, :file_path) ? col.file_path : ""
        println(io, indent, "Column at offset: ", path, "#", col.file_offset)
        show(io, col.meta_data, indent * "  ")
    end
end

function show(io::IO, grp::RowGroup, indent::AbstractString="")
    println(io, indent, "Row Group: ", grp.num_rows, " rows in ", grp.total_byte_size, " bytes")
    show(io, grp.columns, indent * "  ")
end

function show(io::IO, row_groups::Vector{RowGroup}, indent::AbstractString="")
    println(io, indent, "Row Groups:")
    for grp in row_groups
        show(io, grp, indent * "  ")
    end
end

function show(io::IO, meta::FileMetaData, indent::AbstractString="")
    println(io, indent, "version: ", meta.version)
    println(io, indent, "nrows: ", meta.num_rows)
    println(io, indent, "created by: ", meta.created_by)

    show(io, meta.schema, indent)
    show(io, meta.row_groups, indent)
    hasproperty(meta, :key_value_metadata) && show(io, meta.key_value_metadata, indent)
end

function show(io::IO, par::Parquet.File)
    println(io, "Parquet file: $(par.path)")
    meta = par.meta
    println(io, "    version: $(meta.version)")
    println(io, "    nrows: $(meta.num_rows)")
    println(io, "    created by: $(meta.created_by)")
    println(io, "    cached: $(length(par.page_cache.refs)) column chunks")
end

function show(io::IO, table::Parquet.Table)
    print(io, "Parquet.Table(\"$(getfield(table, :path))\")")
end

function show(io::IO, dataset::Parquet.Dataset)
    print(io, "Parquet.Dataset(\"$(getfield(dataset, :path))\")")
end


text:
abstract type Engine end
abstract type OSAbstraction end

findfileswithextension(::OSAbstraction, ::String, ::String) = error("override this method")
readfile(::OSAbstraction, ::String) = error("override this method")
fileexists(::OSAbstraction, ::String) = error("override this method")

struct ExecutorEngine <: Engine
    accumulator::ResultAccumulator
    executor::Executor
    matcher::StepDefinitionMatcher

    function ExecutorEngine(realtimepresenter::RealTimePresenter;
                            executionenv=NoExecutionEnvironment())
        matcher = CompositeStepDefinitionMatcher()
        executor = Executor(matcher, realtimepresenter; executionenv=executionenv)
        new(ResultAccumulator(), executor, matcher)
    end
end

addmatcher!(engine::ExecutorEngine, matcher::StepDefinitionMatcher) = addmatcher!(engine.matcher, matcher)
function runfeature!(engine::ExecutorEngine, feature::Feature)
    result = executefeature(engine.executor, feature)
    accumulateresult!(engine.accumulator, result)
end

finish(engine::ExecutorEngine) = engine.accumulator

struct Driver
    os::OSAbstraction
    engine::Engine

    Driver(os::OSAbstraction, engine::Engine) = new(os, engine)
end

function readstepdefinitions!(driver::Driver, path::String)
    stepdefinitionfiles = findfileswithextension(driver.os, path, ".jl")
    for f in stepdefinitionfiles
        addmatcher!(driver.engine, FromMacroStepDefinitionMatcher(readfile(driver.os, f)))
    end
end

function runfeatures!(driver::Driver, path::String)
    featurefiles = findfileswithextension(driver.os, path, ".feature")
    for featurefile in featurefiles
        featureparseresult = parsefeature(readfile(driver.os, featurefile))
        feature = featureparseresult.value
        runfeature!(driver.engine, feature)
    end
    finish(driver.engine)
end


text:
#=

	This file calculates various aircraft parameters
	horizontal stabilizer coefficient				August 20, 2021
	vertical stabilizer coefficient					September 11, 2021
	initial unit conversion, AR, Sw and MAC functions		October 1, 2021
	added AC calculations						October 7, 2021

	This code is licensed under MIT license (see LICENSE.txt for details)

=#

using DataFrames, ArgParse, Printf					# add functions
import XLSX								# add the Excel access routines

# data structure for keeping track of variables
mutable struct req_vars
		name  ::String
		value ::Float32
		units ::String
		stored::Int8
	end

include("aircraft_dynamics.jl")						# basic routines for dynamics calculations
include("aircraft_drag.jl")						# basic drag calculations

#=  ###########################################################################################
	Parse the command line
	Useage: aircraft_calculations arguments
	Arguments:
			input_file_name is required
			-h or --help for help
			-o or --output_file filename (default no output)
			-s or --sheet_name is the name of the sheet to retrieve the data from
=#

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--output_file", "-o"
            default = "null"
        "--sheet_name", "-s"
            default = "critical_data"
        "input_file"
            help = "Excel spreadsheet name to read data from."
            required = true
    end

    return parse_args(s)
end

#   ###########################################################################################
#	Main function

function main()

    parsed_args = parse_commandline() # parse the input arguments

	# assign the parsed arguments
	AC_data = get(parsed_args, "input_file", "null")						# get the name of the input Excel filename
	output_file = get(parsed_args, "output_file", "null")					# get the name of the output file
	sheet_name = get(parsed_args, "sheet_name", "critical_data") 			# get the name of the sheet within the spreadsheet

	# read the Excel data into a matrix
	xf = XLSX.readxlsx(AC_data)												# open the aircraft data file
	df = (DataFrame(XLSX.readtable(AC_data, sheet_name)...))				# read the spreadsheet data into a data frame
	df_size=size(df)														# read the dataframe size

	if df_size[2] < 3														# check for at least 3 columns:  symbol, value, units
		print("\n\tError: Inadequate number of columns in spreadsheet!\n")	# print the error message if not enough columns
		exit(86)
	end

	## Call the routines that do the work on the input array and append result to data frame
	## Note that the order of these subroutines matters since later ones rely on results of
	## earlier calculations

	push!(df, ["Cg",  		Cg(df),			"m",   		"Center of gravity"])

	push!(df, ["cMAC",  	MAC(df),		"m",   		"Mean wing chord"])
	push!(df, ["Sw",    	Sw(df), 		"m^2", 		"Wing area"])
	push!(df, ["AR",    	AR(df),  		" ",   		"Wing aspect ratio"])
	push!(df, ["AC",    	AC(df),  		"m",   		"Wing AC from LE"])
	push!(df, ["d",			d(df),   		"m",   		"Wing MAC Distance to root"])

	push!(df, ["cMACh", 	MACh(df), 		"m",   		"Horizontal stab mean chord"])
	push!(df, ["Sht",   	Sht(df), 		"m^2", 		"Horizontal stab area"])
	push!(df, ["ARh",   	ARh(df), 		" ",   		"Horizontal stab aspect ratio"])
	push!(df, ["ACh",   	ACh(df),  		"m",   		"Horizontal stab AC from LE"])

	push!(df, ["cMACv", 	MACh(df), 		"m",   		"Vertical stab mean chord"])
	push!(df, ["Svt",   	Svt(df), 		"m^2", 		"Vertical stab area"])
	push!(df, ["ARv",   	ARv(df), 		" ",   		"Vertical stab aspect ratio"])
	push!(df, ["ACv",   	ACv(df),  		"m",   		"Vertical stab AC from LE"])

	push!(df, ["Lht",   	Lht(df), 		"m",   		"Length ACs of wing to h stab"])
	push!(df, ["Lvt",   	Lvt(df), 		"m",   		"Length ACs of wing to v stab"])
	push!(df, ["Cht",   	Cht(df), 		" ",   		"Coefficient of horizontal tail"])
	push!(df, ["Cvt",   	Cvt(df), 		" ",   		"Coefficient of vertical tail"])

	tup_1 = WL(df)
	push!(df, ["W",			tup_1[1],  		"kg",  		"Wing load (level flight)"])
	push!(df, ["Wh",		tup_1[2],  		"kg",  		"Horizontal stab load"])
	push!(df, ["Wload",    	Wload(df),  	"kg/m^2",   "Wing loading"])

	push!(df, ["Vstall", 	Vstall(df), 	"m/s", 		"Wing stall speed"])
	push!(df, ["Vhstall", 	Vhstall(df),	"m/s", 		"Horizontal stab stall speed"])

#	push!(df, ["c", 		c(df),			"m/s", 		"Speed of Sound"])

	drag_sim(df)

	# Print the dataframe
	print("\n\n", df)

end

#   ###########################################################################################
#	Main function (called internally)

main()

print("\n\n\t>> End of Script <<\n")

#   ###########################################################################################
#	End


text:
using Documenter, StructuredOptimization,
LinearAlgebra, DSP, FFTW, AbstractOperators, ProximalAlgorithms

makedocs(
  modules = [StructuredOptimization,ProximalAlgorithms],
  format = Documenter.HTML(),
  sitename = "StructuredOptimization",
  authors = "Niccolò Antonello and Lorenzo Stella",
  pages = [
  "Home"                  => "index.md",
  "Quick Tutorial Guide"  => "tutorial.md",
  "Expressions"           => "expressions.md",
  "Functions"             => "functions.md",
  "Solvers"               => "solvers.md",
  "Demos"                 => "demos.md",
  ],
)

deploydocs(
  repo   = "github.com/kul-forbes/StructuredOptimization.jl.git",
  target = "build",
)


text:
module Scatt

##########################################################################

module Auxiliaries
export read_data_from_file, save_data_to_file

function read_data_from_file(FileName::String)

	dataString = read(FileName, String)

	dataString = replace(dataString, "\r\n" => "\n")
	dataString = replace(dataString, "\r" => "\n")
	dataString = replace(dataString, "\t" => "")
	dataString = replace(dataString, " " => "")
	dataString = String.(split(dataString, "\n"; keepempty = false))
	dataString = dataString[2:end]

	qValues = zeros(length(dataString))
	intensities = zeros(length(qValues))

	for i in 1:length(dataString)
		s = dataString[i]
		qValues[i] = parse(Float64, split(s, ",")[1])
		intensities[i] = parse(Float64, split(s, ",")[2])
	end

	return qValues::Array{Float64, 1}, intensities::Array{Float64, 1}

end

function save_data_to_file(data::Tuple{Array{Float64, 1}, Array{Float64, 1}}, outputFileName::String)

	io = open(outputFileName, "a")
	if isempty(read(outputFileName, String))
		println(io, "q-Value [nm^-1], Intensity [nm^-1]")
	end

	for i in 1:length(data[1])

		local qValue = data[1][i]
		local intensity = data[2][i]

		println(io, qValue, ", ", intensity)

	end

	close(io)

	return nothing

end

end # Auxiliaries

##########################################################################

module Finitizer
export finitize_data, intensity_at_zero
using ..Auxiliaries, FFTW, Interpolations

function expand_data(data::Tuple{Array{Float64, 1}, Array{Float64, 1}}, qValues::Array{Float64, 3})

	intensity = zeros(size(qValues))
	intensity[:] = LinearInterpolation(data[1], data[2]; extrapolation_bc = ((Line(), Flat()),))[qValues[:]]

	return intensity::Array{Float64, 3}

end

function finitize_data(filePath::String, qRange::Array{Float64, 1}, sampleRadiusNm::Float64)

	data = read_data_from_file(filePath)
	(qs, ints) = finitize_data(data, qRange, sampleRadiusNm)
	rm("../output/finitized.data"; force = true)
	save_data_to_file((qs, ints), "../output/finitized.data")

	return nothing

end

function finitize_data(data::Tuple{Array{Float64, 1}, Array{Float64, 1}}, qs::Array{Float64, 1}, sampleRadiusNm::Float64)

	println("\nFinitizing data...")

	function reduce_data(q::Array{Float64, 3}, data::Array{Float64, 3}, px_per_inv_nm::Int, qRange::Array{Float64, 1})

		sideLength = size(data, 1)
		midpoint = Int((sideLength-1)/2 + 1)

		lowQ = [q[midpoint:midpoint+10, midpoint:midpoint+10, midpoint:midpoint+10][:] data[midpoint:midpoint+10, midpoint:midpoint+10, midpoint:midpoint+10][:]]
		lowQ = sortslices(lowQ; dims = 1)
		lowQ = lowQ[indexin(unique(lowQ[:, 1]), lowQ[:, 1]), :]
		lowQ = lowQ[lowQ[:, 1] .< 11/px_per_inv_nm, :]
		rs = [lowQ[:, 1]; ((11:Int((sideLength-1)/2)) ./ px_per_inv_nm)]
		intensities = [lowQ[:, 2]; data[midpoint+11:end, midpoint, midpoint]]
		intensities = LinearInterpolation(rs, intensities)[qRange]

		return intensities::Array{Float64, 1}

	end

	scale = 1
	r = sampleRadiusNm
	background = 0

	px_per_inv_nm = round(Int, 2*r)

	sideLength = 401
	x1 = repeat(1:sideLength, 1, sideLength, sideLength)
	x2 = permutedims(x1, [2 1 3])
	x3 = permutedims(x1, [3 2 1])
	midpoint = Int((sideLength-1)/2 + 1)
	q = sqrt.( (x1.-midpoint).^2 + (x2.-midpoint).^2 + (x3.-midpoint).^2 ) ./ px_per_inv_nm # nm^-1

	qr = q .* r
	mask = scale .* ((3*1).*(sin.(qr) - qr .* cos.(qr)) ./ qr.^3).^2 .+ background # nm^2
	mask[midpoint, midpoint, midpoint] = 1

	mask[abs.(q .* r) .> 98.950063] .= 0

	mask = mask ./ sum(mask)

	intensity = expand_data(data, q)

	convolution = abs.(ifftshift(ifft(fft(mask) .* fft(intensity))))

	ints = reduce_data(q, convolution, px_per_inv_nm, qs)

	println("...done.")

	return qs::Array{Float64, 1}, ints::Array{Float64, 1}

end

function intensity_at_zero(filePath::String)

	data = read_data_from_file(filePath)

	intensity = 0
	for i in 2:min(10, length(data[1]))
		intensity += data[2][1] + (data[2][i] - data[2][1])/(data[1][i] - data[1][1]) * (0 - data[1][1])
	end
	intensity /= min(10, length(data[1])) - 1

	return intensity::Float64

end

end # Finitizer

##########################################################################

module Scatterer
export get_sampsum_for_offset, measure
using Serialization

function load_meta(samplePath::String)

	let dataString = read(samplePath*"scale.meta", String)

		dataString = replace(dataString, "\r\n" => "\n")
		dataString = replace(dataString, "\r" => "\n")
		dataString = replace(dataString, "\t" => "")
		dataString = replace(dataString, " " => "")
		dataString = String.(split(dataString, "\n"; keepempty = false))

		global PX_PER_NM = parse(Float64, split(dataString[1], ":")[2])
		global SAMPLE_RADIUS_NM = parse(Float64, split(dataString[2], ":")[2])

		global SAMPLE_RADIUS_PX = PX_PER_NM * SAMPLE_RADIUS_NM

	end

	let dataString = read(samplePath*"offsets.meta", String)

		dataString = replace(dataString, "\r\n" => "\n")
		dataString = replace(dataString, "\r" => "\n")
		dataString = replace(dataString, "\t" => "")
		dataString = replace(dataString, " " => "")
		dataString = String.(split(dataString, "\n"; keepempty = false))

		global SLD_OFFSETS = parse.(Float64, dataString[2:end])

	end

	return nothing

end

function make_probes(qs::Array{Tuple{Float64, Float64, Float64}, 1})

	probes = Array{Array{Complex{Float16}, 3}, 1}()

	x1 = Float16.(repeat(1:Int(2*SAMPLE_RADIUS_PX+1), 1, Int(2*SAMPLE_RADIUS_PX+1), Int(2*SAMPLE_RADIUS_PX+1)))
	x2 = permutedims(x1, [2 1 3])
	x3 = permutedims(x1, [3 2 1])

	for i in 1:length(qs)
		q = qs[i]
		probe = Complex{Float16}.(exp.((x1.*q[1]+x2.*q[2]+x3.*q[3]).*im))
		probes = [probes; [probe]]
		println(i, "/", length(qs))
	end

	return probes::Array{Array{Complex{Float16}, 3}, 1}

end

function get_intensity(sample::Array{Float16, 3}, probe::Array{Complex{Float16}, 3})

	sample = deepcopy(sample)
	sample[sample .== -Inf] .= 0

	intensity = abs(sum(Float64.(sample) .* Complex{Float64}.(probe) .* (1/PX_PER_NM)^3))^2

	return intensity::Float64

end

function measure(samplePath::String)

	load_meta(samplePath)

	sampleVolume = (4 * pi * SAMPLE_RADIUS_NM^3 / 3)

	numMeasur = 3
	for k in 1:numMeasur

		println("\nSetting up measurement ", k, " of ", numMeasur, "...\n")

		global qs = Array{Tuple{Float64, Float64, Float64}, 1}()
		rs = 10 .^ (-1:0.01:1) # q-range [nm^-1]
		for r in rs
			local q = (1, 1, 1)
			while sqrt(sum(q.^2)) > 1
				q = (2*rand()-1, 2*rand()-1, 2*rand()-1)
			end
			q = r .* (q ./ sqrt(sum(q.^2)))
			q = q ./ PX_PER_NM # pixel wave vector
			global qs = [qs; [q]]
		end

		probes = make_probes(qs)

		println("\n...done.")

		intensities = zeros(length(rs))

		n = 1
		while isfile(samplePath*string(n, pad = 5)*".samp")

			sample = deserialize(samplePath*string(n, pad = 5)*".samp")
			sample = sample .+ Float16(SLD_OFFSETS[n])

			println("\nq-value [nm^-1], Intensity [nm^-1]", " # Measurement ", k, " of ", numMeasur, ".")
			for i in 1:length(rs)
				intensities[i] += get_intensity(sample, probes[i])
				local intensity = intensities[i]
				local qValue = sqrt(sum((qs[i].*PX_PER_NM).^2))
				println(qValue, ", ", intensity/sampleVolume, " # "*string(n, pad = 5)*" [", i, "/", length(rs), "]")
			end

			n += 1

		end
		n -= 1

		intensities = intensities ./ (n * sampleVolume)

		outputFileName = "../output/simulated.data"
		io = open(outputFileName, "a")
		if isempty(read(outputFileName, String))
			println(io, "q-Value [nm^-1], Intensity [nm^-1]")
		end

		for i in 1:length(rs)

			local intensity = intensities[i]
			local qValue = sqrt(sum((qs[i].*PX_PER_NM).^2))

			println(io, qValue, ", ", intensity)

		end

		close(io)

	end

	return nothing

end

end # Scatterer

##########################################################################

module Plotter
export make_plot
using ..Auxiliaries, ..Finitizer, Plots

function make_plot()

	gr()

	plt = plot()

	(qValues, intensities) = read_data_from_file("../output/simulated.data")
	plot!(qValues, intensities;
		seriestype = :scatter,
		xaxis = :log,
		yaxis = :log,
		minorgrid = true,
		markershape = :circle,
		markersize = 2,
		markercolor = RGB(0, 0, 1),
		markerstrokecolor = RGB(0, 0, 1),
		markeralpha = 0.5,
		xlabel = "q-Value [nm\u207B\u00B9]",
		ylabel = "Intensity [nm\u207B\u00B9]",
		foreground_color_legend = nothing,
		background_color_legend = nothing,
		lab = "Simulated",
		framestyle = :box)

	(qValues, intensities) = read_data_from_file("../output/finitized.data")
	plot!(qValues, intensities;
		seriestype = :line,
		linewidth = 3,
		linecolor = RGB(0, 1, 0),
		linealpha = 0.5,
		lab = "Calculated")

	savefig(plt, "../output/figure.pdf")

	return nothing

end

end # Plotter

##########################################################################

end # Scatt

max_stars_repo_name: JonPallbo/Scatt

max_issues_repo_name: JonPallbo/Scatt

max_forks_repo_name: JonPallbo/Scatt


text:
import Base: +, setindex!, filter!

# ------------ Parameter store ------------ #

struct ParameterStore
    params::Dict{Address,Any}
    ParameterStore() = new(Dict{Address, Any}())
    ParameterStore(d::Dict{Address, Any}) = new(d)
end
haskey(ps::ParameterStore, addr) = haskey(ps.params, addr)
setindex!(ps::ParameterStore, val, addr) = ps.params[addr] = val

Zygote.@adjoint ParameterStore(params) = ParameterStore(params), store_grad -> (nothing,)

function +(a::ParameterStore, b::ParameterStore)
    params = Dict{Address, Any}()
    for (k, v) in Iterators.flatten((a.params, b.params))
        if !haskey(params, k)
            params[k] = v
        else
            params[k] += v
        end
    end
    ParameterStore(params)
end

# ------------ Backpropagation contexts ------------ #

abstract type BackpropagationContext <: ExecutionContext end

# Learnable parameters
mutable struct ParameterBackpropagateContext{T <: Trace} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::Parameters
    params::ParameterStore
    param_grads::Gradients
end
ParameterBackpropagate(tr::T, init, params) where T <: Trace = ParameterBackpropagateContext(tr, 0.0, init, params, Gradients())
ParameterBackpropagate(tr::T, init, params, param_grads::Gradients) where {T <: Trace, K <: UnconstrainedSelection} = ParameterBackpropagateContext(tr, 0.0, init, params, param_grads)

# Choice sites
mutable struct ChoiceBackpropagateContext{T <: Trace, K <: UnconstrainedSelection} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::Parameters
    params::ParameterStore
    choice_grads::Gradients
    select::K
end
ChoiceBackpropagate(tr::T, init, params, choice_grads) where {T <: Trace, K <: UnconstrainedSelection} = ChoiceBackpropagateContext(tr, 0.0, init, params, choice_grads, UnconstrainedAllSelection())
ChoiceBackpropagate(tr::T, init, params, choice_grads, sel::K) where {T <: Trace, K <: UnconstrainedSelection} = ChoiceBackpropagateContext(tr, 0.0, init, params, choice_grads, sel)

# ------------ Learnable ------------ #

read_parameter(ctx::K, addr::Address) where K <: BackpropagationContext = read_parameter(ctx, ctx.params, addr)
read_parameter(ctx::K, params::ParameterStore, addr::Address) where K <: BackpropagationContext = get_top(ctx.initial_params, addr)

Zygote.@adjoint function read_parameter(ctx, params, addr)
    ret = read_parameter(ctx, params, addr)
    fn = param_grad -> begin
        state_grad = nothing
        params_grad = ParameterStore(Dict{Address, Any}(addr => param_grad))
        (state_grad, params_grad, nothing)
    end
    return ret, fn
end

# ------------ Call sites ------------ #

# Grads for learnable parameters.
simulate_call_pullback(params, param_grads, cl::T, args) where T <: CallSite = cl.ret

Zygote.@adjoint function simulate_call_pullback(params, param_grads, cl::HierarchicalCallSite, args)
    ret = simulate_call_pullback(params, param_grads, cl, args)
    fn = ret_grad -> begin
        arg_grads = accumulate_parameter_gradients!(params, param_grads, cl, ret_grad)
        (nothing, nothing, nothing, arg_grads)
    end
    return ret, fn
end

# Utility.
merge(tp1::Tuple{}, tp2::Tuple{}) = tp1
merge(tp1::Tuple{Nothing}, tp2::Tuple{Nothing}) where T = tp1
merge(tp1::NTuple{N, Float64}, tp2::NTuple{N, Float64}) where N = [tp1[i] + tp2[i] for i in 1 : N]
merge(tp1::Array{Float64}, tp2::NTuple{N, Float64}) where N = [tp1[i] + tp2[i] for i in 1 : N]

Zygote.@adjoint function simulate_call_pullback(params, param_grads, cl::VectorizedCallSite{typeof(markov)}, args)
    ret = simulate_call_pullback(params, param_grads, cl, args)
    fn = ret_grad -> begin
        arg_grads = accumulate_parameter_gradients!(params, param_grads, get_sub(cl, cl.len), ret_grad)
        for i in (cl.len - 1) : -1 : 1
            arg_grads = accumulate_parameter_gradients!(params, param_grads, get_sub(cl, i), arg_grads)
        end
        (nothing, nothing, nothing, arg_grads)
    end
    return ret, fn
end

Zygote.@adjoint function simulate_call_pullback(params, param_grads, cl::VectorizedCallSite{typeof(plate)}, args)
    ret = simulate_call_pullback(params, param_grads, cl, args)
    fn = ret_grad -> begin
        arg_grads = accumulate_parameter_gradients!(params, param_grads, get_sub(cl, cl.len), ret_grad[1])
        for i in 2 : cl.len
            new = accumulate_parameter_gradients!(params, param_grads, get_sub(cl, i), ret_grad[i])
            arg_grads = merge(arg_grads, new)
        end
        (nothing, nothing, nothing, arg_grads)
    end
    return ret, fn
end

# Grads for choices with differentiable logpdfs.
simulate_choice_pullback(params, choice_grads, choice_selection, cl::T, args) where T <: CallSite = cl.ret

Zygote.@adjoint function simulate_choice_pullback(params, choice_grads, choice_selection, cl, args)
    ret = simulate_choice_pullback(params, choice_grads, choice_selection, cl, args)
    fn = ret_grad -> begin
        arg_grads, choice_vals, choice_grads = choice_gradients(params, choice_grads, choice_selection, cl, ret_grad)
        (nothing, nothing, nothing, (choice_vals, choice_grads), arg_grads)
    end
    return ret, fn
end

# ------------ Accumulate gradients ------------ #

function accumulate_parameter_gradients!(initial_params, param_grads, cl::HierarchicalCallSite, ret_grad, scaler::Float64 = 1.0)
    fn = (args, params) -> begin
        ctx = ParameterBackpropagate(cl.trace, initial_params, params, param_grads)
        ret = ctx(cl.fn, args...)
        (ctx.weight, ret)
    end
    blank = ParameterStore()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, ps_grad = back((1.0, ret_grad))
    if !(ps_grad isa Nothing)
        for (addr, grad) in ps_grad.params
            push!(param_grads, addr, scaler .* grad)
        end
    end
    return arg_grads
end

function accumulate_parameter_gradients!(initial_params, param_grads, cl::HierarchicalCallSite, ret_grad::Tuple, scaler::Float64 = 1.0)
    fn = (args, params) -> begin
        ctx = ParameterBackpropagate(cl.trace, initial_params, params, param_grads)
        ret = ctx(cl.fn, args...)
        (ctx.weight, ret)
    end
    blank = ParameterStore()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, ps_grad = back((1.0, ret_grad...))
    if !(ps_grad isa Nothing)
        for (addr, grad) in ps_grad.params
            push!(param_grads, addr, scaler .* grad)
        end
    end
    return arg_grads
end

function accumulate_parameter_gradients!(initial_params, param_grads, cl::VectorizedCallSite{typeof(markov)}, ret_grad, scaler::Float64 = 1.0) where T <: CallSite
    addr = gensym()
    v_params = parameters(addr => initial_params)
    tr = Trace()
    add_call!(tr, addr, cl)
    fn = (args, params) -> begin
        ctx = ParameterBackpropagate(tr, v_params, params, param_grads)
        ret = ctx(markov, addr, cl.fn, cl.len, args...)
        (ctx.weight, ret)
    end
    blank = ParameterStore()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, ps_grad = back((1.0, ret_grad))
    if !(ps_grad isa Nothing)
        for (addr, grad) in ps_grad.params
            push!(param_grads, addr, scaler .* grad)
        end
    end
    return arg_grads
end

function accumulate_parameter_gradients!(initial_params, param_grads, cl::VectorizedCallSite{typeof(plate)}, ret_grad, scaler::Float64 = 1.0) where T <: CallSite
    addr = gensym()
    v_params = parameters(addr => initial_params)
    tr = Trace()
    add_call!(tr, addr, cl)
    fn = (args, params) -> begin
        ctx = ParameterBackpropagate(tr, v_params, params, param_grads)
        ret = ctx(plate, addr, cl.fn, args)
        (ctx.weight, ret)
    end
    blank = ParameterStore()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, ps_grad = back((1.0, ret_grad))
    if !(ps_grad isa Nothing)
        for (addr, grad) in ps_grad.params
            push!(param_grads, addr, scaler .* grad)
        end
    end
    return arg_grads
end

# ------------ Compute choice gradients ------------ #

function filter!(choice_grads, cl::HierarchicalCallSite, grad_tr::NamedTuple, sel::K) where K <: UnconstrainedSelection
    values = ConstrainedHierarchicalSelection()
    for (k, v) in cl.trace.choices
        has_top(sel, k) && begin
            push!(values, k, v.val)
            push!(choice_grads, k, grad_tr.choices[k].val)
        end
    end
    return values
end

function choice_gradients(initial_params, choice_grads, choice_selection, cl, ret_grad)
    call = cl.fn
    fn = (args, trace) -> begin
        ctx = ChoiceBackpropagate(trace, initial_params, ParameterStore(), choice_grads, choice_selection)
        ret = ctx(call, args...)
        (ctx.weight, ret)
    end
    _, back = Zygote.pullback(fn, cl.args, cl.trace)
    arg_grads, grad_ref = back((1.0, ret_grad))
    choice_vals = filter!(choice_grads, cl, grad_ref, choice_selection)
    return arg_grads, choice_vals, choice_grads
end

# ------------ Convenience getters ------------ #

function get_choice_gradients(cl::T, ret_grad) where T <: CallSite
    choice_grads = Gradients()
    choice_selection = UnconstrainedAllSelection()
    _, vals, _ = choice_gradients(Parameters(), choice_grads, choice_selection, cl, ret_grad)
    return vals, choice_grads
end

function get_choice_gradients(params, cl::T, ret_grad) where T <: CallSite
    choice_grads = Gradients()
    choice_selection = UnconstrainedAllSelection()
    _, vals, _ = choice_gradients(params, choice_grads, choice_selection, cl, ret_grad)
    return vals, choice_grads
end

function get_choice_gradients(sel::K, cl::T, ret_grad) where {T <: CallSite, K <: UnconstrainedSelection}
    choice_grads = Gradients()
    _, vals, _ = choice_gradients(Parameters(), choice_grads, sel, cl, ret_grad)
    return vals, choice_grads
end

function get_choice_gradients(sel::K, params, cl::T, ret_grad) where {T <: CallSite, K <: UnconstrainedSelection}
    choice_grads = Gradients()
    _, vals, _ = choice_gradients(params, choice_grads, sel, cl, ret_grad)
    return vals, choice_grads
end

function get_parameter_gradients(params, cl::HierarchicalCallSite, ret_grad, scaler::Float64 = 1.0)
    param_grads = Gradients()
    accumulate_parameter_gradients!(params, param_grads, cl, ret_grad, scaler)
    return param_grads
end

function get_parameter_gradients(params, cl::VectorizedCallSite, ret_grad, scaler::Float64 = 1.0)
    param_grads = Gradients()
    accumulate_parameter_gradients!(params, param_grads, cl, ret_grad, scaler)
    for k in keys(param_grads.tree)
        return param_grads.tree[k]
    end
end

# ------------ train ------------ #

function train(params, fn::Function, args...; opt = ADAM(0.05, (0.9, 0.8)), iters = 1000)
    for i in 1 : iters
        _, cl = simulate(params, fn, args...)
        grads = get_parameter_gradients(params, cl, 1.0)
        params = update_parameters(opt, params, grads)
    end
    return params
end

function train(sel::K, params, fn::Function, args...; opt = ADAM(0.05, (0.9, 0.8)), iters = 1000) where K <: ConstrainedSelection
    for i in 1 : iters
        _, cl, _ = generate(sel, params, fn, args...)
        grads = get_parameter_gradients(params, cl, 1.0)
        params = update_parameters(opt, params, grads)
    end
    return params
end

# ------------ includes ------------ #

include("hierarchical/backpropagate.jl")
include("plate/backpropagate.jl")
include("markov/backpropagate.jl")
include("factor/backpropagate.jl")

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct ParameterBackpropagateContext{T <: Trace} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::Parameters
    params::ParameterStore
    param_grads::Gradients
end
```
`ParameterBackpropagateContext` is used to compute the gradients of parameters with respect to following objective:

Outer constructors:
```julia
ParameterBackpropagate(tr::T, params) where T <: Trace = ParameterBackpropagateContext(tr, 0.0, params, Gradients())
ParameterBackpropagate(tr::T, params, param_grads::Gradients) where {T <: Trace, K <: UnconstrainedSelection} = ParameterBackpropagateContext(tr, 0.0, params, param_grads)
```
""", ParameterBackpropagateContext)

@doc(
"""
```julia
mutable struct ChoiceBackpropagateContext{T <: Trace} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::Parameters
    params::ParameterStore
    param_grads::Gradients
end
```
`ChoiceBackpropagateContext` is used to compute the gradients of choices with respect to following objective:

Outer constructors:
```julia
ChoiceBackpropagate(tr::T, init_params, params, choice_grads) where {T <: Trace, K <: UnconstrainedSelection} = ChoiceBackpropagateContext(tr, 0.0, params, choice_grads, UnconstrainedAllSelection())
ChoiceBackpropagate(tr::T, init_params, params, choice_grads, sel::K) where {T <: Trace, K <: UnconstrainedSelection} = ChoiceBackpropagateContext(tr, 0.0, params, choice_grads, sel)
```
""", ChoiceBackpropagateContext)

@doc(
"""
```julia
gradients = get_choice_gradients(params, cl::T, ret_grad) where T <: CallSite
gradients = get_choice_gradients(cl::T, ret_grad) where T <: CallSite
```

Returns a `Gradients` object which tracks the gradients with respect to the objective of random choices with differentiable `logpdf` in the program.
""", get_choice_gradients)

@doc(
"""
```julia
gradients = get_parameter_gradients(params, cl::T, ret_grad, scaler::Float64 = 1.0) where T <: CallSite
```

Returns a `Gradients` object which tracks the gradients of the objective with respect to parameters in the program.
""", get_parameter_gradients)


text:
include("./Agilent4294A.jl")
include("./Agilent4395A.jl")

"""
    get_frequency_limits(instr)

# Returns
`Tuple{Frequency, Frequency}`: (lower_limit, upper_limit)
"""
get_frequency_limits(obj::Instr{T}) where (T <: ImpedanceAnalyzer) =
    f_query(obj, "STAR?") * u"Hz", f_query(obj, "STOP?") * u"Hz"

"""
    set_frequency_limits(instr, lower_limit, upper_limit)

"""
function set_frequency_limits(
    obj::Instr{T},
    start::Frequency,
    stop::Frequency
) where {T <: ImpedanceAnalyzer}
    start = raw(start)
    stop = raw(stop)
    write(obj, "STAR $start; STOP $stop")
end

"""
    set_num_data_points(instr, num_points)

"""
set_num_data_points(obj::Instr{T}, n) where (T <: ImpedanceAnalyzer) =
    write(obj, "POIN $n")

"""
    get_num_data_points(instr)

"""
get_num_data_points(obj::Instr{T}) where (T <: ImpedanceAnalyzer) =
    i_query(obj, "POIN?")

"""
    get_volt_dc(instr)

"""
get_volt_dc(obj::Instr{T}) where (T <: ImpedanceAnalyzer) =
    f_query(obj, "DCV?") * V

"""
    set_volt_dc(instr, volts)

"""
set_volt_dc(obj::Instr{T}, num::Voltage) where (T <: ImpedanceAnalyzer) =
    write(obj, "DCV $(raw(num))")


text:
function testExpminus(p::Float64, λ::Float64, m::Int64)
  testFunction(BernoulliFactory.expMinus, m, exp(-λ*p), makef(p), λ, GLOBAL_RNG)
end

Random.seed!(192837465)

testExpminus(0.2, 1.0, ntrials)
testExpminus(0.8, 3.0, ntrials)
testExpminus(0.3, 0.5, ntrials)


text:
using BinDeps

vers = "1.5.0"

tagfile = "installed_vers"
target = "libblosc.$(Sys.dlext)"
if !isfile(tagfile) || readchomp(tagfile) != vers
    if OS_NAME == :Windows
        run(download_cmd("http://ab-initio.mit.edu/blosc/libblosc$WORD_SIZE-$vers.dll", target))
    elseif OS_NAME == :Darwin
        run(download_cmd("http://ab-initio.mit.edu/blosc/libblosc$WORD_SIZE-$vers.dylib", target))
    else
        tarball = "c-blosc-$vers.tar.gz"
        srcdir = "c-blosc-$vers/blosc"
        if !isfile(tarball)
            run(download_cmd("https://github.com/Blosc/c-blosc/archive/v$vers.tar.gz", tarball))
        end
        run(unpack_cmd(tarball, ".", ".gz", ".tar"))
        cd(srcdir) do
            println("Compiling libblosc...")
            for f in ("blosc.c", "blosclz.c", "shuffle.c")
                println("   CC $f")
                run(`gcc -fPIC -O3 -msse2 -I. -c $f`)
            end
            println("   LINK libblosc")
            run(`gcc -shared -o ../../$target blosc.o blosclz.o shuffle.o`)
        end
    end
    run(`echo $vers` |> tagfile)
end


text:
using StingerGraphs
using BenchmarkTools
using JLD
using FileIO

function bench(filename)
    s = Stinger()
    fields = [
        max_nv, max_neblock, max_netypes, max_nvtypes, num_insertions,
        num_deletions, num_insertions_last_batch, num_deletions_last_batch
    ]

    getrng = srand(1)

    getbench = @benchmarkable $s[f] samples=10000000  seconds=10 setup=(f=rand($getrng, $fields))

    setrng = srand(2)
    setbench = @benchmarkable $s[f]=val samples=10000000 setup=f,val=rand($setrng, $fields),rand($setrng,0:1000) seconds=10#For a reasonable run time

    info("Running getbench")
    gettrial = run(getbench)
    info("Running setbench")
    settrial = run(setbench)

    jldopen(filename, "w") do f
        write(f, "gettrial", gettrial)
        write(f, "settrial", settrial)
    end

    return gettrial, settrial
end


text:
struct PredictModel{S, F}
    f::F
end

"""
    PredictModel{S}(f::Function)

Create a prediction model for use in a [`BasicParticleFilter`](@ref)

See [`ParticleFilterModel`](@ref) for descriptions of `S` and `f`.
"""
PredictModel{S}(f::F) where {S, F<:Function} = PredictModel{S, F}(f)

function predict!(pm, m::PredictModel, b, u, rng)
    for i in 1:n_particles(b)
        x1 = particle(b, i)
        pm[i] = m.f(x1, u, rng)
    end
end

particle_memory(m::PredictModel{S}) where S = S[]

"""
    ReweightModel(g::Function)

Create a reweighting model for us in a [`BasicParticleFilter`](@ref).

See [`ParticleFilterModel`](@ref) for a description of `g`.
"""
struct ReweightModel{G}
    g::G
end

function reweight!(wm, m::ReweightModel, b, u, pm, y)
    for i in 1:n_particles(b)
        x1 = particle(b, i)
        x2 = pm[i]
        wm[i] = m.g(x1, u, x2, y)
    end
end

struct ParticleFilterModel{S, F, G}
    f::F
    g::G
end

"""
    ParticleFilterModel{S}(f, g)

Create a system model suitable for use in a particle filter. This is a combination prediction/dynamics model and reweighting model.

# Parameters
- `S` is the state type, e.g. `Float64`, `Vector{Float64}`

# Arguments
- `f` is the dynamics function. `xₜ₊₁ = f(xₜ, uₜ, rng)` where `x` is the state, `u` is the control input, and `rng` is a random number generator.
- `g` is the observation weight function. `g(xₜ, uₜ, xₜ₊₁, yₜ₊₁)` returns the likelihood weight of measurement yₜ₊₁ given that the state transitioned from `xₜ` to `xₜ₊₁` when control uₜ was applied. These weights need not be normalized (and, for performance, usually should not be).
"""
function ParticleFilterModel{S}(f::F, g::G) where {S, F<:Function, G<:Function}
    return ParticleFilterModel{S, F, G}(f, g)
end

function predict!(pm, m::ParticleFilterModel{S}, b, u, rng) where S
    predict!(pm, PredictModel{S}(m.f), b, u, rng)
end

function reweight!(wm, m::ParticleFilterModel, b, u, pm, y)
    reweight!(wm, ReweightModel(m.g), b, u, pm, y)
end

particle_memory(m::ParticleFilterModel{S}) where S = S[]


text:
test_image_results = [0, 1, 1, 4, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 3, 1, 2, 2, 2]

# Black pixels
test_image_1 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

# White pixels
test_image_2 = [
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1 ;
1 1 1 1 1]

# Single pixel
test_image_3 = [
0 0 0 0 0 ;
0 0 0 0 0 ;
0 0 1 0 0 ;
0 0 0 0 0 ;
0 0 0 0 0]

# Multiple single pixels
test_image_4 = [
0 0 0 0 0 ;
0 1 0 1 0 ;
0 0 0 0 0 ;
0 1 0 1 0 ;
0 0 0 0 0]

# Checkerboard 1
test_image_5 = [
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1]

# Checkerboard 2
test_image_6 = [
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0 ;
1 0 1 0 1 ;
0 1 0 1 0]

# External contour
test_image_7 = [
0 0 0 0 0 ;
0 1 1 1 0 ;
0 1 0 1 0 ;
0 1 1 1 0 ;
0 0 0 0 0]

# Cut-off external contour
test_image_8 = [
0 0 0 0 0 ;
0 1 1 1 0 ;
0 1 0 1 0 ;
0 1 0 1 0 ;
0 1 0 1 0]

# Multiple external contour
test_image_9 = [
0 0 0 0 0 ;
1 1 1 1 0 ;
0 0 0 0 0 ;
0 1 1 1 1 ;
0 0 0 0 0]

# Single pixel internal contour
test_image_10 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 1 0 1 1 ;
1 1 1 1 1 ;
0 1 1 1 0]

# Internal contour
test_image_11 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 0 0 0 1 ;
1 1 1 1 1 ;
0 1 1 1 0]

# Multiple internal contours
test_image_12 = [
1 0 1 0 1 ;
1 1 1 1 1 ;
1 0 1 0 1 ;
1 1 1 1 1 ;
0 1 0 1 0]

# Labelled internal contours
test_image_13 = [
1 1 1 1 1 ;
1 0 0 1 1 ;
1 0 0 0 1 ;
1 1 0 1 1 ;
1 1 1 1 1]

# Unlabelled internal contours
test_image_14 = [
1 1 1 1 1 ;
1 1 1 1 1 ;
1 0 0 1 1 ;
1 0 0 1 1 ;
1 1 1 1 1]

# No dummy row
test_image_15 = [
0 0 0 0 0 ;
1 1 1 1 1 ;
1 0 1 0 1 ;
1 0 1 0 1 ;
1 1 1 1 1]

# Multiple external and internal contours 1
test_image_16 = [
0 0 0 0 0 0 0 0 0 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 0 0 0 0 1 0;
0 1 0 1 1 1 1 1 1 0;
0 1 1 1 0 0 0 0 0 0;
0 0 0 0 0 1 1 1 1 0;
0 1 1 1 0 1 0 0 1 0;
0 1 0 1 0 1 0 0 1 0;
0 1 1 1 0 1 1 1 1 0;
0 0 0 0 0 0 0 0 0 0]

# Multiple external and internal contours 2
test_image_17 = [
0 0 0 0 0 0 0 0 0 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 1 1 0 0 1 0;
0 1 0 0 0 0 1 0 1 0;
0 1 1 0 0 1 0 0 1 0;
0 1 1 1 1 0 0 1 0 0;
0 0 0 0 0 0 0 0 1 0;
0 1 0 0 0 0 1 1 1 0;
0 1 1 1 1 1 1 1 0 0;
0 0 0 0 0 0 0 0 0 0]

# Multiple external and internal contours 3
test_image_18 = [
0 0 0 0 0 0 0 0 0 0;
0 0 1 1 1 1 1 0 0 0;
0 1 0 1 0 1 1 1 0 0;
0 1 0 0 1 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0;
0 1 1 1 1 0 0 0 0 0;
0 0 1 0 0 0 1 1 0 0;
0 0 0 0 0 1 0 0 1 0;
0 0 0 0 0 1 0 0 1 0;
0 0 0 0 0 0 1 0 1 0]

# Multiple external and internal contours 4
test_image_19 = [
0 0 0 0 1 1 0 0 0 0;
0 0 0 1 0 0 1 0 0 0;
0 0 1 0 0 0 0 1 0 0;
0 1 0 0 1 1 0 0 1 0;
1 0 0 1 0 0 1 0 0 1;
1 0 0 1 0 0 1 0 0 1;
0 1 0 1 0 0 1 0 1 0;
0 1 0 0 1 1 0 0 1 0;
0 0 1 0 0 0 0 1 0 0;
0 0 0 1 1 1 1 0 0 0]

# Multiple external and internal contours 5
test_image_20 = [
0 1 0 0 0 0 0 0 1 0;
0 1 1 1 1 1 1 1 1 0;
0 1 0 0 1 0 0 0 1 0;
0 1 0 0 0 1 0 0 1 0;
0 1 0 1 0 0 1 1 1 0;
0 1 0 0 1 0 0 0 1 0;
0 1 0 1 0 1 0 0 1 0;
0 1 0 1 0 0 1 0 1 0;
0 1 0 0 1 1 1 0 1 0;
0 1 0 0 1 0 0 0 1 0]


text:
# Load Packages
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Geometry
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.AdditiveRungeKuttaMethod
using CLIMA.LinearSolvers
using CLIMA.GeneralizedMinimalResidualSolver
using CLIMA.SubgridScaleParameters
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
using Random
using CLIMA.Atmos: vars_state, vars_aux

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray,)
else
  const ArrayTypes = (Array,)
end

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

# -------------- Problem constants ------------------- #
const (xmin,xmax)      = (0,1000)
const (ymin,ymax)      = (0,400)
const (zmin,zmax)      = (0,1000)
const Ne        = (10,2,10)
const polynomialorder = 4
const dim       = 3
const dt        = 0.1
const timeend   = 10dt
# ------------- Initial condition function ----------- #
"""
@article{doi:10.1175/1520-0469(1993)050<1865:BCEWAS>2.0.CO;2,
author = {Robert, A},
title = {Bubble Convection Experiments with a Semi-implicit Formulation of the Euler Equations},
journal = {Journal of the Atmospheric Sciences},
volume = {50},
number = {13},
pages = {1865-1873},
year = {1993},
doi = {10.1175/1520-0469(1993)050<1865:BCEWAS>2.0.CO;2},
URL = {https://doi.org/10.1175/1520-0469(1993)050<1865:BCEWAS>2.0.CO;2},
eprint = {https://doi.org/10.1175/1520-0469(1993)050<1865:BCEWAS>2.0.CO;2},
}
"""
function Initialise_Rising_Bubble!(state::Vars, aux::Vars, (x1,x2,x3), t)
  FT            = eltype(state)
  R_gas::FT     = R_d
  c_p::FT       = cp_d
  c_v::FT       = cv_d
  γ::FT         = c_p / c_v
  p0::FT        = MSLP

  xc::FT        = 500
  zc::FT        = 260
  r             = sqrt((x1 - xc)^2 + (x3 - zc)^2)
  rc::FT        = 250
  θ_ref::FT     = 303
  Δθ::FT        = 0

  if r <= rc
    Δθ          = FT(1//2)
  end
  #Perturbed state:
  θ            = θ_ref + Δθ # potential temperature
  π_exner      = FT(1) - grav / (c_p * θ) * x3 # exner pressure
  ρ            = p0 / (R_gas * θ) * (π_exner)^ (c_v / R_gas) # density
  P            = p0 * (R_gas * (ρ * θ) / p0) ^(c_p/c_v) # pressure (absolute)
  T            = P / (ρ * R_gas) # temperature
  ρu           = SVector(FT(0),FT(0),FT(0))
  # energy definitions
  e_kin        = FT(0)
  e_pot        = grav * x3
  ρe_tot       = ρ * total_energy(e_kin, e_pot, T)
  state.ρ      = ρ
  state.ρu     = ρu
  state.ρe     = ρe_tot
  state.moisture.ρq_tot = FT(0)
end
# --------------- Driver definition ------------------ #
function run(mpicomm, ArrayType, LinearType,
             topl, dim, Ne, polynomialorder,
             timeend, FT, dt)
  # -------------- Define grid ----------------------------------- #
  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = polynomialorder
                                           )
  # -------------- Define model ---------------------------------- #
  model = AtmosModel(FlatOrientation(),
                     HydrostaticState(IsothermalProfile(FT(T_0)),FT(0)),
                     Vreman{FT}(C_smag),
                     EquilMoist(),
                     NoRadiation(),
                     Gravity(),
                     NoFluxBC(),
                     Initialise_Rising_Bubble!)
  # -------------- Define dgbalancelaw --------------------------- #
  dg = DGModel(model,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  linmodel = LinearType(model)
  lindg = DGModel(linmodel,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty(); auxstate=dg.auxstate)

  Q = init_ode_state(dg, FT(0))

  linearsolver = GeneralizedMinimalResidual(10, Q, sqrt(eps(FT)))
  ark = ARK548L2SA2KennedyCarpenter(dg, lindg, linearsolver, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e
  ArrayType = %s
  FloatType = %s""" eng0 ArrayType FT

  # Set up the information callback (output field dump is via vtk callback: see cbinfo)
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(10, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(ark),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  step = [0]
  cbvtk = GenericCallbacks.EveryXSimulationSteps(3000)  do (init=false)
    mkpath("./vtk-rtb/")
      outprefix = @sprintf("./vtk-rtb/DC_%dD_mpirank%04d_step%04d", dim,
                           MPI.Comm_rank(mpicomm), step[1])
      @debug "doing VTK output" outprefix
      writevtk(outprefix, Q, dg, flattenednames(vars_state(model,FT)), param[1], flattenednames(vars_aux(model,FT)))
      step[1] += 1
      nothing
  end

  solve!(Q, ark; timeend=timeend, callbacks=(cbinfo,cbvtk))
  # End of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, FT(timeend))
  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
engf/eng0
end
# --------------- Test block / Loggers ------------------ #
using Test
let
  MPI.Initialized() || MPI.Init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
      device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end
  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    FloatType = (Float32, Float64)
    for FT in FloatType
      brickrange = (range(FT(xmin); length=Ne[1]+1, stop=xmax),
                    range(FT(ymin); length=Ne[2]+1, stop=ymax),
                    range(FT(zmin); length=Ne[3]+1, stop=zmax))
      topl = StackedBrickTopology(mpicomm, brickrange, periodicity = (false, true, false))
      for LinearType in (AtmosAcousticLinearModel, AtmosAcousticGravityLinearModel)
        engf_eng0 = run(mpicomm, ArrayType, LinearType,
                        topl, dim, Ne, polynomialorder,
                        timeend, FT, dt)
        @test engf_eng0 ≈ FT(0.9999997771981113)
      end
    end
  end
end

#nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# # Assigning concrete function call to a method saves the function, not the result.
# module functions
# function factorial(n)
#     if n>1
#         return n*factorial(n-1)
#     elseif n==1
#         return 1
#     end
# end

# using BenchmarkTools

# fn = factorial(15)

# special() = factorial(15)
# special2() = fn

# @benchmark factorial(15)
# @benchmark special()
# @benchmark special2()

# end

# # define vector interface for custom struct
# module vectorinterface_poc

# using StaticArrays

# struct mysv{Size} <: StaticArray{Tuple{Size}, Bool, 1}
#     v::SVector{Size, Bool}
#     meta::String
# end

# @inline Tuple(m::mysv) = m.v.data
# StaticArrays.@propagate_inbounds function Base.getindex(m::mysv, i::Int)
#     m.v.data[i]
# end

# a = mysv{2}(SVector(false, false), "aha")
# b = SVector(true, true)

# @show b+a
# end

# # traits
# module traits
# using SimpleTraits

# abstract type Superbtype end

# struct maintype <: Superbtype
#     a::Int
# end

# struct otherype <: Superbtype
#     a::Int
#     b::String
# end
# end

# # Comparing speed of rng
# using Random
# using RandomNumbers.Xorshifts
# using BenchmarkTools

# rngm = MersenneTwister(1)
# rngx = Xorshift128Plus(1)
# rngx2 = Xorshift1024Plus(1)

# @benchmark MersenneTwister(1)
# @benchmark Xorshift128Plus(1)
# @benchmark Xorshift1024Plus(1)

# @benchmark rand(rngm)
# @benchmark rand(rngx)
# @benchmark rand(rngx2)

# # Try-catch block scope issues

# function test_trycatch()
#     a = "outer"
#     try
#         a = "try"
#         throw(ErrorException("Oh no!"))
#     catch
#         a = "catch"
#     end

#     a
# end

# test_trycatch()

# # plotting functions
# using Plots

# # f(x) = log.(x)

# x = 1:100
# plot(x, .(log.(x)))
# plot(x, sqrt.(log.(x)))
# plot!(x, sqrt.(log.(x)/10))
# plot!(x, sqrt.(log.(x)/100))

# C,α = (1, 0.1)
# plot(x,C.*x.^α)

# # BSON Dict bug?
# using BSON
# using FileIO

# inner_d = Dict(:v=>0)
# outer_d = Dict(:d=>inner_d)
# save("test.bso", outer_d)
# display(inner_d)
# load("test.bson")

# v = [zeros(1) for i in 1:2]
# D = Dict(:d=>v)
# save("test_v.bson", D)
# BSON.load("test_v.bson")

# JLD issues
using JLD
using StaticArrays

struct MyType{N}
    a::SVector{N, Int64}
end

MyType{2}(SA[1,2])

struct MyType2{T}
    a::SVector{2, T}
end

struct MyType3{T}
    a::T
end

struct MyType4{T}
    b::MyType3{T}
end

MyType4{Int64}(MyType3{Int64}(1))

JLD.save("test.jld", "s", MyType4{Int64}(MyType3{Int64}(1)))
load("test.jld")


text:
using JuAFEM
using MOOSE

using ForwardDiff
using ForwardDiff: Partials, Dual, value, partials, npartials, setindex

using Base.Test

if MPI.Comm_size(MPI.COMM_WORLD) == 1
#    include("test_Node.jl")
#    include("test_Element.jl")
#    include("test_Generation.jl")
#    include("test_Variable.jl")
#    include("test_System.jl")
#    include("test_Solver.jl")
#    include("test_Kernel.jl")
#    include("test_Assembly.jl")
#
#    include("test_Diffusion.jl")
#
#    include("test_DirichletBC.jl")
#
#    include("test_VTKOutput.jl")
#
#    include("test_JuliaDenseImplicitSolver.jl")
#
#    if MOOSE.have_petsc
#       include("test_PetscImplicitSolver.jl")
#    end
#
#    include("test_Mesh.jl")
#    include("test_PerfLog.jl")
    include("test_XDAReader.jl")
end

if MPI.Comm_size(MPI.COMM_WORLD) == 2
#    include("test_parallel_Mesh.jl")
#    include("test_parallel_System.jl")
end

if MOOSE.have_petsc
#   include("test_MetisPartitioner.jl")
end


text:
using Pkg
Pkg.activate(".")

include(joinpath(@__DIR__, "..", "src", "trg.jl"))
include(joinpath(@__DIR__, "..", "src", "2d_classical_ising.jl"))

# Make Ising model MPO
β = 1.1 * βc
d = 2
s = Index(d)
sₕ = addtags(s, "horiz")
sᵥ = addtags(s, "vert")
T = ising_mpo(sₕ, sᵥ, β)

χmax = 20
nsteps = 20
κ, T = trg(T; χmax = χmax, nsteps = nsteps,
           svd_alg = "divide_and_conquer")

κ_exact = exp(-β * ising_free_energy(β))
@show κ, κ_exact
@show abs(κ - κ_exact)


text:
println("===============================================")
println("Testing interface.jl")
println("===============================================")

function _interface_Hertz_f()
    itf = AFM(4.0)
    ramp = RheoTimeData(itf, t = Vector(0:0.1:10), d = Vector(0:0.05:5))
    model=RheoModel(Spring, k=2)
    data = modelpredict(ramp,model)

    # check expected force for d = 1.0
    isapprox(data[itf].f[21], 16/3., atol=tol)
end
@test _interface_Hertz_f()

function _interface_Hertz_d()
    itf = AFM(2)
    ramp = RheoTimeData(itf, t = Vector(0:0.1:10), f = Vector(0:0.05:5))
    model=RheoModel(Spring, k=0.25)
    data = modelpredict(ramp,model)

    # check expected displacement for f = 1.0
    isapprox(data[itf].d[21], 2*(0.75^(2/3)), atol=tol)
end
@test _interface_Hertz_d()

function _interface_Tweezers_f()
    itf = Tweezers(2.0)
    ramp = RheoTimeData(itf, t = Vector(0:0.1:10), d = Vector(0:0.05:5))
    model=RheoModel(Spring, k=2)
    data = modelpredict(ramp,model)

    # check expected force for d = 1.0
    isapprox(data[itf].f[21], 4*6*π, atol=tol)
end
@test _interface_Tweezers_f()

function _interface_Tweezers_d()
    itf = Tweezers(2.0,0) # assumes full slip
    ramp = RheoTimeData(itf, t = Vector(0:0.1:10), f = Vector(0:0.05:5))
    model=RheoModel(Spring, k=0.25)
    data = modelpredict(ramp,model)

    # check expected displacement for f = 1.0
    isapprox(data[itf].d[21], 2/(4*π), atol=tol)
end
@test _interface_Tweezers_d()

function _interface_importcsv()
    fildir = joinpath(@__DIR__, "testdata", "interfacetest.csv")
    itf = Interface(:d, :f, (ϵ,σ)->(d = 2 .* ϵ, f = 2 .* σ ), (d,f)->(ϵ = 0.5 .* d, σ = 0.5 .* f) )
    data = importcsv(fildir, itf, d="d", f="f", t="t")

    # check expected displacement for d = 0.2
    isapprox(data.ϵ[3], 0.1, atol=tol)
end
@test _interface_importcsv()


text:
load(::Day{1}, path) = parse.(Int, (readlines(path)))

# Part 1 – count increases
solve1(::Day{1}, data) = sum(x -> x > 0, diff(data))

# Part 2: count windowed increases ignoring identical summands in window
function solve2(::Day{1}, data; window=3)
    n = length(data)
    return sum(view(data, (1 + window):n) .> view(data, 1:(n - window)))
end

testresult1(::Day{1}) = 7
testresult2(::Day{1}) = 5


text:
module Clines

abstract type Cline end
export Cline

const _DEFAULT_TOL = 1e-10
_TOLERANCE = _DEFAULT_TOL

export set_tolerance, get_tolerance

"""
    set_tolerance(tol)
Set the tolerance to Cline roundoff errors.
"""
function set_tolerance(tol::Real = _DEFAULT_TOL)
    if tol <= 0
        error("Tolerance cannot be set to $tol; must be positive")
    end
    global _TOLERANCE = tol
end
"""
    get_tolerance(tol)
Return the current tolerance to Cline roundoff errors.
"""
function get_tolerance()
    return _TOLERANCE
end

"""
    collinear(a,b,c)
Determine if the three points `a`, `b`, and `c` are collinear.
Here the inputs are finite complex numbers.
"""
function collinear(a::Number, b::Number, c::Number)::Bool
    x1, y1 = reim(a)
    x2, y2 = reim(b)
    x3, y3 = reim(c)

    d = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)

    return abs(d) <= get_tolerance()
end

export collinear

is_huge(z::Number) = isinf(z) || abs(z) >= 1 / get_tolerance()

function Cline(a::Number, b::Number, c::Number)::Cline
    if isnan(a) || isnan(b) || isnan(c)
        error("Cannot construct a Cline with a NaN argument")
    end

    two_inf = "Cannot construct a Cline with two infinite arguments"
    two_dif = "Cannot construct a Cline through a single point"

    # If any of the arguments is ∞ then we make a line

    if is_huge(a)
        if is_huge(b) || is_huge(c)
            error(two_inf)
        end
        b != c || error(two_dif)
        return Line(b, c)
    end

    if is_huge(b)   # know that a is not infinite
        a != c || error(two_dif)
        return Line(a, c)
    end

    if is_huge(c)   # we know a and b are not infinite
        a != b || error(two_dif)
        return Line(a, b)
    end

    # if two of the arguments are equal, we make a line

    if a == b == c     # but not if three are equal
        error(two_dif)
    end

    if a == b
        return Line(a, c)
    end
    if a == c
        return Line(a, b)
    end
    if b == c
        return Line(a, b)
    end

    # if the arguments are collienar, we make a line
    if collinear(a, b, c)
        x, y, z = sort_points(a, b, c)
        return Line(x, z)
    end

    # Otherwise, it's a circle
    return Circle(Complex(a), Complex(b), Complex(c))
end

Cline(a::Number, b::Number) = Cline(a, b, Inf)

import SimpleDrawing: draw

include("Lines.jl")
include("Circles.jl")
include("intersection.jl")
include("angle.jl")

# applying LFT to a Cline

function (F::LFT)(C::Cline)
    pp = three_points(C)
    qq = F.(pp)
    return Cline(qq...)
end

# create a LFT from a pair of Clines

import LinearFractionalTransformations: LFT

"""
    LFT(C::Cline, D::Cline)
Return a linear fractional transformation that maps `C` to `D`. If `D` is
omitted, it is assumed to be the x-axis.
"""
function LFT(C::Cline, D::Cline)::LFT
    F = LFT(C)
    G = LFT(D)
    return inv(G) * F
end

function LFT(C::Circle)::LFT
    return LFT(three_points(C)...)
end

function LFT(L::Line)::LFT
    return LFT(L.a, L.b, Inf)
end

end # module


text:
export ISWAP, SqrtX, SqrtY, SqrtW, singlet_block
export ISWAPGate, SqrtXGate, SqrtYGate, SqrtWGate, CPhaseGate

const CPhaseGate{T} = ControlBlock{<:ShiftGate{T},<:Any}

@const_gate ISWAP = PermMatrix([1,3,2,4], [1,1.0im,1.0im,1])
@const_gate SqrtX = [0.5+0.5im 0.5-0.5im; 0.5-0.5im 0.5+0.5im]
@const_gate SqrtY = [0.5+0.5im -0.5-0.5im; 0.5+0.5im 0.5+0.5im]
# √W is a non-Clifford gate
@const_gate SqrtW = mat(rot((X+Y)/sqrt(2), π/2))

"""
    singlet_block(θ::Real, ϕ::Real)

The circuit block for initialzing a singlet state.
"""
singlet_block() = chain(put(2, 1=>chain(X, H)), control(2, -1, 2=>X))


text:
# This module contains the abstractions for the solver based on the Symbolics package

include("_equilibrium_matrix_conditions.jl")

"""

# Examples
```
using Symbolics
@variables N P T N_A MU_B NP_ALPHA N
elements = ["A", "B", "C"]
phases = ["ALPHA", "BETA", "GAMMA"]
condition_dict = Dict(
   N => 1.0,
   P => 101325.0,
   T => 300.0,
   N_A => 0.5,
   MU_B => -10000,
   NP_ALPHA => 0.8,
)
unpack_indices(elements, phases, condition_dict)
```
"""
function unpack_indices(elements, phases, conditions_keys; validate_phases=false)
    elements = sort(elements)
    POTENTIALS = ("P", "T",)

    fixed_chempot_symbols = []
    fixed_chempot_elements = []
    fixed_pot_symbols = []
    fixed_pot_strings = []
    fixed_phase_symbols = []
    fixed_phase_names = []
    condition_row_symbols = []
    # TODO: move validation outside this function
    for cond in conditions_keys
        str_cond = string(cond)
        if startswith(str_cond, "MU_")
            el = str_cond[4:end]
            if el ∉ elements
                throw("Element $el in condition $str_cond is not in the elements $elements")
            end
            push!(fixed_chempot_symbols, cond)
            push!(fixed_chempot_elements, el)
        elseif str_cond in POTENTIALS
            push!(fixed_pot_symbols, cond)
            push!(fixed_pot_strings, str_cond)
        elseif startswith(str_cond, "ℵ_")
            phase_name = str_cond[5:end]  # TODO: not sure why it needs to be 5
            # handle case if unpack_indices is used for a single phase Δy
            if validate_phases
                if phase_name ∉ phases
                    throw("Phase $phase_name in condition $str_cond is not in the phases $phases")
                end
            end
            push!(fixed_phase_symbols, cond)
            push!(fixed_phase_names, phase_name)
        elseif startswith(str_cond, "N_") | startswith(str_cond, "X_") | (str_cond in ("N", ))
            # TODO: validate N_ and X_ conditions that they contain elements
            push!(condition_row_symbols, cond)
        else
            throw("Unknown condition $cond")
        end
    end

    # The fixed conditions are all in sorted order. Now we have to determine the
    # indices of the free conditions
    free_chempot_idxs = [A for A in 1:length(elements) if elements[A] ∉ fixed_chempot_elements]
    free_pot_idxs = [pot for pot in 1:length(POTENTIALS) if POTENTIALS[pot] ∉ fixed_pot_strings]
    free_phase_idxs = [α for α in 1:length(phases) if phases[α] ∉ fixed_phase_names]

    return (
        condition_row_symbols,
        fixed_chempot_symbols, free_chempot_idxs,
        fixed_pot_symbols, free_pot_idxs,
        fixed_phase_symbols, free_phase_idxs,
    )

end

function get_solution_parts(phase_records, elements, conditions_keys)
    conditions_keys = collect(conditions_keys)
    elements = sort(elements)
    phases = [prx.phase_name for prx in phase_records]
    idxs = unpack_indices(elements, phases, conditions_keys)
    condition_row_symbols = first(idxs)
    fixed_free_soln_terms = Base.tail(idxs)

    # construct rows and right-hand-side
    free_chempots = fixed_free_soln_terms[2]
    free_pots = fixed_free_soln_terms[4]
    free_phases = fixed_free_soln_terms[6]

    neqns = length(phase_records) + length(condition_row_symbols)
    soln_size = length(free_chempots) + length(free_pots) + length(free_phases)
    @assert neqns == soln_size "The number of matrix equations ($neqns) does not match the number of solution terms ($soln_size)"

    A = Matrix{Num}(undef, neqns, soln_size)
    b = Array{Num}(undef, soln_size)

    # Get the equations
    # One row for each stable phase:
    for phase_idx in 1:length(phase_records)
        row, rhs = get_stable_phase_row_rhs(phase_records, phase_idx, fixed_free_soln_terms...)
        A[phase_idx, :] = row
        b[phase_idx] = rhs
    end
    # One row for each condition equation
    row_offset = length(phase_records)
    for cond_idx in 1:length(condition_row_symbols)
        cond = condition_row_symbols[cond_idx]
        row, rhs = cond_row_rhs(cond, elements, phases, phase_records, fixed_free_soln_terms)
        A[row_offset+cond_idx, :] = row
        b[row_offset+cond_idx] = rhs
    end

    return A, b
end

function get_solution(phase_records, elements, conditions_keys)
    A, b = get_solution_parts(phase_records, elements, conditions_keys)
    soln = A \ b
    return soln
end

@doc raw"""
    get_Delta_y_mat

Notationally, Δy is a vector of length(phase_record.site_fractions) that
updates the site fractions. However in this case, we need to plug in the
results from the solution vector that we do not have symbolic variables for. To
resolve that, we'll design the site fractions to be a matrix that can be used
to take the dot product with the solution vector (padded with a prefixed one to
handle the c_iG term). The usage is therefore

# Examples
```
delta_y_M = get_Delta_y_mat(prx, ["A", "B"], [T, P, N_A, N_B])
# ... compute solution
soln = [3271.04833019, 7271.04833015, 1e-16]
delta_y = delta_y_M * vcat(1, soln...)
```

# Equation
```math
\sum_i \Delta y_i^\alpha = \sum_i \left(c_{iG}^\alpha + \sum_A c_{iA}^\alpha \mu_A + \sum_\mathrm{Pot} c_{i\mathrm{Pot}}^\alpha \Delta \mathrm{Pot} \right)
```
"""
function get_Delta_y_mat(phase_record, elements, conditions_keys)
    conditions_keys = collect(conditions_keys)
    elements = sort(elements)
    idxs = unpack_indices(elements, [phase_record.phase_name], conditions_keys; validate_phases=false)
    fixed_free_soln_terms = Base.tail(idxs)
    num_free_chempots = length(fixed_free_soln_terms[2])
    num_free_pots = length(fixed_free_soln_terms[4])

    Δys = Matrix{Num}(undef, length(phase_record.site_fractions), 1+num_free_chempots+num_free_pots)
    for i in 1:length(phase_record.site_fractions)
        Δys[i, 1] = c_iG(phase_record, i)
        offset = 1
        for A in 1:num_free_chempots
            Δys[i, offset+A] = c_iA(phase_record, i, A)
        end
        offset += num_free_chempots
        for pot in 1:num_free_pots
            Δys[i, offset+pot] = c_iPot(phase_record, i, pot)
        end
    end
    return Δys
end

function get_subs_dict(compsets, conditions_dict)
    d = Dict{Num, Num}()
    for compset in compsets
        for (Y_sym, Y_num) in zip(compset.phase_rec.site_fractions, compset.Y)
            d[Y_sym] = Num(Y_num)
        end
        d[compset.phase_rec.ℵ] = compset.ℵ
    end
    for cond in keys(conditions_dict)
        d[cond] = conditions_dict[cond]
    end
    return d
end

function isphasecond(cond::Num)
    str_cond = string(cond)
    return startswith(str_cond, "ℵ_")
end

function vectorize_inputs(compsets::Vector{CompSet}, free_potentials::OrderedDict{Num,Float64}, conditions::OrderedDict{Num,Float64})
    compset_conds = [vcat(cs.phase_rec.site_fractions, [cs.phase_rec.ℵ]) for cs in compsets]
    return vcat(collect(keys(free_potentials)), [ky for ky in keys(conditions) if !isphasecond(ky)], compset_conds...)
end

function vectorize_values(compsets::Vector{CompSet}, free_potentials::OrderedDict{Num,Float64}, conditions::OrderedDict{Num,Float64})
    compset_conds = [vcat(cs.Y, [cs.ℵ]) for cs in compsets]
    return vcat(collect(values(free_potentials)), [val for (ky, val) in conditions if !isphasecond(ky)], compset_conds...)
end

function build_callable(expr, inputs)
    # The function at index 1 is a callable with a return, the second one is for in-place modification
    f = eval(build_function(expr, inputs)[1])
    return f
end

function build_delta_y_callables(delta_y_matrices::Vector{Matrix{Num}}, soln::Vector{Num}, inp::Vector{Num}, num_free_phases::Int)
    soln_potentials = vcat([1], soln[1:end-num_free_phases])
    funcs = Vector{Function}(undef, length(delta_y_matrices))
    for i in 1:length(delta_y_matrices)
        delta_y_expr = delta_y_matrices[i] * soln_potentials
        # The function at index 1 is a callable with a return, the second one is for in-place modification
        funcs[i] = build_callable(delta_y_expr, inp)
    end
    return funcs
end

# TODO: test that phase amount updating works properly. Most I've seen so far
# keep Δℵ close to zero (single phase, but I think phase amount should still)
# change because the mass condition RHS shouldn't be satisfied, that is:
# (`(N_A - Ñ_A) != 0`)
function solve_and_update(compsets::Vector{CompSet}, free_potentials::OrderedDict{Num,Float64}, conditions::OrderedDict{Num,Float64}, soln_func::Function, delta_y_funcs::Vector{Function}, free_phase_idxs::Vector{Int}; step_size=1.0, doprint=false)
    num_free_phases = length(free_phase_idxs)
    x = vectorize_values(compsets, free_potentials, conditions)
    if doprint
        println(x)
    end
    soln = soln_func(x)
    # Extract chemical potentials
    num_free_chempots = length(soln) - num_free_phases - length(free_potentials)
    chempots = soln[1:num_free_chempots]
    if doprint
        println("equilibrium_soln: ", soln)
    end
    # Update ΔPot
    # TODO: assumes free potentials are in the same order as the potentials in the solution vector, is there a way to enforce this structurally without incurring overhead time?
    for (i, (ky, val)) in enumerate(free_potentials)
        # offset for number of free chemical potentials
        free_potentials[ky] = val + soln[num_free_chempots+i]
    end
    # Update Δℵ
    for β in 1:num_free_phases
        α = free_phase_idxs[β]
        Δℵ = soln[end-num_free_phases+β]
        if doprint
            println("$(compsets[α].phase_rec.phase_name): Δℵ=$(Δℵ)")
        end
        compsets[α].ℵ = max(compsets[α].ℵ+Δℵ, 0.0)
    end
    # Update Δy
    for α in 1:length(compsets)
        Δy = delta_y_funcs[α](x)
        if doprint
            println("$(compsets[α].phase_rec.phase_name): Δy=$(Δy) (step size=$step_size)")
        end
        compsets[α].Y += step_size*Δy
    end
    return chempots
end


text:
"""
    PreprocessedDataCausalityTest <: MetaCausalityTest

A causality test where the data is subjected to some preprocessing before
applying the test.
"""
abstract type PreprocessedDataCausalityTest{CT} <: MetaCausalityTest{CT} end

export PreprocessedDataCausalityTest

include("BinnedDataCausalityTest.jl")
include("ConstrainedTest.jl")
include("InterpolateBinTest.jl")

include("causality_BinnedDataCausalityTest.jl")
include("causality_ConstrainedTest.jl")
include("causality_InterpolateBinTest.jl")


text:
include("scaling_nthreads_target.jl")
using Setfield: @set!

if abspath(PROGRAM_FILE) == @__FILE__
    let benchmarks = copy(BENCHMARKS)
        @set! benchmarks[:sort].paramaxes.datasize =
            [100_000, 1_000_000, 10_000_000, 100_000_000]
        @set! benchmarks[:foreach_symmetrize].paramaxes.datasize = [1000, 10_000]

        main(benchmark_definitions = benchmarks, tags = ["datasize"])
    end
end


text:
function test_noc(bridged_mock, F, S, n)
    @test MOI.get(bridged_mock, MOI.NumberOfConstraints{F, S}()) == n
    @test length(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{F, S}())) == n
    @test ((F, S) in MOI.get(bridged_mock, MOI.ListOfConstraints())) == !iszero(n)
end

# Test deletion of bridge
function test_delete_bridge(
    m::MOIB.AbstractBridgeOptimizer, ci::MOI.ConstraintIndex{F, S}, nvars::Int,
    nocs::Tuple; used_bridges = 1, num_bridged = 1) where {F, S}
    function num_bridges()
        return count(bridge -> bridge !== nothing, m.bridges) +
            length(m.single_variable_constraints)
    end
    start_num_bridges = num_bridges()
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    test_noc(m, F, S, num_bridged)
    for noc in nocs
        test_noc(m, noc...)
    end
    @test MOI.is_valid(m, ci)
    MOI.delete(m, ci)
    @test_throws MOI.InvalidIndex{typeof(ci)} MOI.delete(m, ci)
    try
        MOI.delete(m, ci)
    catch err
        @test err.index == ci
    end
    @test !MOI.is_valid(m, ci)
    @test num_bridges() == start_num_bridges - used_bridges
    test_noc(m, F, S, num_bridged - 1)
    # As the bridge has been removed, if the constraints it has created where not removed, it wouldn't be there to decrease this counter anymore
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    for noc in nocs
        test_noc(m, noc...)
    end
end