text:
using Test
using JSON
using Random
import Jupyter2Pluto: PlutoCell, JupyterMarkdownCell, JupyterCodeCell, generate_code

const markdown_arr = [
    "# Data structures\n",
    "\n",
    "Once we start working with many pieces of data at once, it will be convenient for us to store data in structures like arrays or dictionaries (rather than just relying on variables).
\n",
    "\n",
    "Types of data structures covered:\n",
    "1. Tuples\n",
    "2. Dictionaries\n",
    "3. Arrays\n",
    "\n",
    "
\n",
    "As an overview, tuples and arrays are both ordered sequences of elements (so we can index into them). Dictionaries and arrays are both mutable.\n",
    "We'll explain this more below!"
]

##TODO: fix remove \n from the last jupyter cell
const code_arr = ["# First, restore fibonacci\n", "fibonacci[1] = 1\n", "fibonacci"]
Random.seed!(1234)

function test_jupyter2pluto()
    #@info join(markdown_arr)
    jmark_cell = JupyterMarkdownCell(join(markdown_arr))
    pmc = PlutoCell(jmark_cell)
    @testset "test pluto markdown codegen" begin
        @test string(pmc) == """# ╔═╡ 196f2941-2d58-45ba-9f13-43a2532b2fa8
md\"\"\"
# Data structures

Once we start working with many pieces of data at once, it will be convenient for us to store data in structures like arrays or dictionaries (rather than just relying on variables).

Types of data structures covered:
1. Tuples
2. Dictionaries
3. Arrays

As an overview, tuples and arrays are both ordered sequences of elements (so we can index into them). Dictionaries and arrays are both mutable.
We'll explain this more below!\"\"\"

"""
    end
    pcc = PlutoCell(JupyterCodeCell(1, code_arr))
    @testset "Pluto code block codegen" begin
        @test string(pcc) == """# ╔═╡ bd85187e-0531-4a3e-9fea-713204a818a2
begin
    # First, restore fibonacci
    fibonacci[1] = 1
    fibonacci
end

"""
    end
    @test string(generate_code(PlutoCell[pmc, pcc])) == """### A Pluto.jl notebook ###

# ╔═╡ 196f2941-2d58-45ba-9f13-43a2532b2fa8
md\"\"\"
# Data structures

Once we start working with many pieces of data at once, it will be convenient for us to store data in structures like arrays or dictionaries (rather than just relying on variables).

Types of data structures covered:
1. Tuples
2. Dictionaries
3. Arrays

As an overview, tuples and arrays are both ordered sequences of elements (so we can index into them). Dictionaries and arrays are both mutable.
We'll explain this more below!\"\"\"

# ╔═╡ bd85187e-0531-4a3e-9fea-713204a818a2
begin
    # First, restore fibonacci
    fibonacci[1] = 1
    fibonacci
end

# ╔═╡ Cell order:
# ╟─196f2941-2d58-45ba-9f13-43a2532b2fa8
# ╠═bd85187e-0531-4a3e-9fea-713204a818a2
"""
end

test_jupyter2pluto()


text:
using Multigraphs: Multigraph
import ZXCalculus
using ZXCalculus: ZXDiagram, SpiderType, add_spider!, Phase

function ZXCalculus.ZXDiagram(nin::Integer, nout::Integer)
    N = nin + nout
    mg = Multigraph(nin + nout)
    st = vcat([SpiderType.In for _ = 1:nin], [SpiderType.Out for _ = 1:nout])
    ps = [Phase(0//1) for _ = 1:N]
    zxd = ZXDiagram(mg, st, ps)
    sort!(zxd.inputs)
    sort!(zxd.outputs)
    return zxd
end

function ZXCalculus.ZXDiagram(q::Tait)
    simplify!(q, Rule(:genus_fusion))
    nin = length(q.inputs)
    nout = length(q.outputs)
    zxd = ZXDiagram(nin, nout)
    v_map = Dict{Int, Int}()
    for i = 1:nin
        v_map[q.inputs[i]] = i
    end
    for i = 1:nout
        v_map[q.outputs[i]] = i + nin
    end

    for v in vertices(q)
        haskey(v_map, v) && continue
        is_genus(q, v) && continue
        v_map[v] = add_spider!(zxd, SpiderType.Z)
    end

    he_map = Dict{Int, Int}()
    for he in half_edges(q)
        he_twin = twin(q, he)
        s = src(q, he)
        d = dst(q, he)
        if !haskey(he_map, he)
            if haskey(q.quon_params, he)
                p = quon_param(q, he)
                if is_genus(q, s) || is_genus(q, d)
                    v_z = v_map[is_genus(q, s) ? d : s]
                    if real(p.param) == 0
                        if !is_parallel(p)
                            zxd.ps[v_z] += imag(p.param) / pi
                            he_map[he] = v_z
                            he_map[he_twin] = v_z
                        else
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p.param) / pi), [v_z])
                            he_map[he] = v_x
                            he_map[he_twin] = v_x
                        end
                    else    # the parameter is not pure imag number
                        p1 = QuonParam{QuonComplex}(real(p.param) + pi/2*im, is_parallel(p))
                        p2 = QuonParam{QuonComplex}(imag(p.param) - pi/2*im, is_parallel(p))
                        p1 = change_direction(p1)
                        if !is_parallel(p1)
                            zxd.ps[v_z] += imag(p1.param) / pi
                            he_map[he] = v_z
                            he_map[he_twin] = v_z
                        else
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p1.param) / pi), [v_z])
                            he_map[he] = v_x
                            he_map[he_twin] = v_x
                        end
                        if !is_parallel(p2)
                            zxd.ps[v_z] += imag(p2.param) / pi
                            he_map[he] = v_z
                            he_map[he_twin] = v_z
                        else
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p2.param) / pi), [v_z])
                            he_map[he] = v_x
                            he_map[he_twin] = v_x
                        end
                    end
                else
                    if real(p.param) == 0
                        if is_parallel(p)
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p.param) / pi), [v_map[s], v_map[d]])
                            he_map[he] = v_x
                            he_map[he_twin] = v_x
                        else
                            v_x = add_spider!(zxd, SpiderType.X, Phase(0//1), [v_map[s], v_map[d]])
                            v_z = add_spider!(zxd, SpiderType.Z, Phase(imag(p.param) / pi), [v_x])
                            he_map[he] = v_z
                            he_map[he_twin] = v_z
                        end
                    else    # the parameter is not pure imag number
                        p1 = QuonParam{QuonComplex}(real(p.param) + pi/2*im, is_parallel(p))
                        p2 = QuonParam{QuonComplex}(imag(p.param) - pi/2*im, is_parallel(p))
                        p1 = change_direction(p1)
                        if is_parallel(p1)
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p1.param) / pi), [v_map[s], v_map[d]])
                            v_z = add_spider!(zxd, SpiderType.Z, Phase(imag(p2.param) / pi), [v_x])
                            he_map[he] = v_x
                            he_map[he_twin] = v_x
                        else
                            v_x = add_spider!(zxd, SpiderType.X, Phase(imag(p2.param) / pi), [v_map[s], v_map[d]])
                            v_z = add_spider!(zxd, SpiderType.Z, Phase(imag(p1.param) / pi), [v_x])
                            he_map[he] = v_z
                            he_map[he_twin] = v_z
                        end
                    end
                end
            else
                if !(is_genus(q, s) || is_genus(q, d))
                    ZXCalculus.add_edge!(zxd, v_map[s], v_map[d])
                    he_map[he] = v_map[s]
                    he_map[he_twin] = v_map[d]
                end
            end
        end
    end
    return zxd
end

ext: jl


text:
function mdl_shaft_ode!(device_states,
        output_ode,
        f0::Float64,
        device::DynGenerator{M, SingleMass, A, TG, P})  where {M <: Machine,
                                                               A <: AVR,
                                                               TG <: TurbineGov,
                                                               P <: PSS}

    #Obtain references
    ω_ref = get_ω_ref(device)

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.shaft]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    δ = internal_states[1]
    ω = internal_states[2]

    #Obtain inner variables for component
    τe = device.inner_vars[τe_var]
    τm = device.inner_vars[τm_var]

    #Get parameters
    H = device.shaft.H
    D = device.shaft.D

    #Compute 2 states ODEs
    output_ode[local_ix[1]] = 2*π*f0*(ω-ω_ref)                    #15.5
    output_ode[local_ix[2]] = (1/(2*H))*(τm - τe - D*(ω-ω_ref))   #15.5

    return
end

function mdl_shaft_ode!(device_states,
                        output_ode,
                        f0::Float64,
                        device::DynGenerator{M, FiveMassShaft, A, TG, P})  where {M <: Machine,
                                                                                A <: AVR,
                                                                                TG <: TurbineGov,
                                                                                P <: PSS}

    #Obtain references
    ω_ref = get_ω_ref(device)

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.shaft]

    #Define internal states for component
    internal_states = @view device_states[local_ix]
    δ = internal_states[1]
    ω = internal_states[2]
    δ_hp = internal_states[3]
    ω_hp = internal_states[4]
    δ_ip = internal_states[5]
    ω_ip = internal_states[6]
    δ_lp = internal_states[7]
    ω_lp = internal_states[8]
    δ_ex = internal_states[9]
    ω_ex = internal_states[10]

    #Obtain inner variables for component
    τe = device.inner_vars[τe_var]
    τm = device.inner_vars[τm_var]

    #Get parameters
    H = device.shaft.H
    H_hp = device.shaft.H_hp
    H_ip = device.shaft.H_ip
    H_lp = device.shaft.H_lp
    H_ex = device.shaft.H_ex
    D = device.shaft.D
    D_hp = device.shaft.D_hp
    D_ip = device.shaft.D_ip
    D_lp = device.shaft.D_lp
    D_ex = device.shaft.D_ex
    D_12 = device.shaft.D_12
    D_23 = device.shaft.D_23
    D_34 = device.shaft.D_34
    D_45 = device.shaft.D_45
    K_hp = device.shaft.K_hp
    K_ip = device.shaft.K_ip
    K_lp = device.shaft.K_lp
    K_ex = device.shaft.K_ex

    #Compute 10 states ODEs #15.51
    output_ode[local_ix[1]] = 2.0*π*f0*(ω-ω_ref)

    output_ode[local_ix[2]] = (1.0/(2.0*H))*(-τe - D*(ω-ω_ref) - D_34*(ω-ω_lp)
                              -D_45*(ω - ω_ex) + K_lp*(δ_lp - δ)
                              +K_ex*(δ_ex - δ))

    output_ode[local_ix[3]] = 2.0*π*f0*(ω_hp - ω_ref)

    output_ode[local_ix[4]] =  (1.0/(2.0*H_hp))*(τm - D_hp*(ω_hp - ω_ref)
                                - D_12*(ω_hp - ω_ip) +K_hp*(δ_ip - δ_hp))

    output_ode[local_ix[5]] = 2*π*f0*(ω_ip - ω_ref)

    output_ode[local_ix[6]] =  (1.0/(2*H_ip))*(-D_ip*(ω_ip - ω_ref)
                                - D_12*(ω_ip - ω_hp) -D_23*(ω_ip - ω_lp)
                                + K_hp*(δ_hp - δ_ip) +K_ip*(δ_lp - δ_ip))

    output_ode[local_ix[7]] = 2.0*π*f0*(ω_lp - ω_ref)

    output_ode[local_ix[8]] =  (1.0/(2.0*H_lp))*(-D_lp*(ω_lp - ω_ref)
                                - D_23*(ω_lp - ω_ip) -D_34*(ω_lp - ω)
                                + K_ip*(δ_ip - δ_lp) +K_lp*(δ-δ_lp))

    output_ode[local_ix[9]] = 2*π*f0*(ω_ex - ω_ref)

    output_ode[local_ix[10]] = (1.0/(2.0*H_ex))*(-D_ex*(ω_ex - ω_ref)
                                - D_45*(ω_ex - ω) +K_ex*(δ-δ_ex))

    return
end


text:
using PhotoOrganizer

dry_run = false
rm_src = false
dst_root="/home/hertz/Pictures/Pictures"

src_dirs = String[]
#push!(src_dirs, "/run/user/1000/gvfs/mtp:host=%5Busb%3A002%2C009%5D/Samsung SD card/CameraZOOM")
#push!(src_dirs, "/run/user/1000/gvfs/mtp:host=%5Busb%3A002%2C009%5D/Samsung SD card/DCIM/Camera")
#push!(src_dirs, "/home/hertz/Documents/Glen/backups/Google/Takeout/Google Photos")
#push!(src_dirs, "/home/hertz/Documents/Erin/backups/Google/Takeout/Google Photos")
#src_dirs = readlines(`/home/hertz/bin/ls_phone_backup_dirs.jl`)

src_dirs = ["/home/hertz/Documents.local/Pictures"]
src_root = "/home/hertz/Documents/backups/S7"
src_root = "/run/user/1000/gvfs/mtp:host=%5Busb%3A003%2C004%5D"
user = "eh"
if user == "gh"
	src_root = "/var/run/user/1000/gvfs/mtp:host=%5Busb%3A002%2C009%5D"
	src_dirs = [
		    "$src_root/Card/DCIM/Camera",
		    "$src_root/Phone/DCIM/Camera",
		    "$src_root/Phone/DCIM/PhotoScan",
		    "$src_root/Phone/DCIM/Screenshots",
		    "$src_root/Phone/DCIM/Video Editor",
		    "$src_root/Phone/Movies/Instagram",
		    "$src_root/Phone/Pictures",
		    "$src_root/Phone/Snapchat",
		    "$src_root/Phone/Snapseed",
		    "$src_root/Phone/Studio",
		    "$src_root/Phone/Telegram",
		    "$src_root/Phone/WhatsApp/Media"]
elseif user == "eh"
	src_root = "/var/run/user/1000/gvfs/mtp:host=%5Busb%3A002%2C011%5D/Internal shared storage"
	src_dirs = [
		    "$src_root/CameraZOOM",
		    "$src_root/DCIM",
		    "$src_root/Movies",
		    "$src_root/Pictures",
		    "$src_root/Phone/DCIM/Screenshots",
		    "$src_root/Phone/DCIM/Video Editor",
		    "$src_root/Phone/Movies/Instagram",
		    "$src_root/SlowMotion",
		    "$src_root/Snapchat",
		    "$src_root/Telegram",
		    "$src_root/WhatsApp/Media"]
end
if length(src_dirs) > 0
    organize_photos(src_dirs, dst_root, rm_src, dry_run)
else
    warn("No directories found to backup")
end


text:
using PkgSRA
using Plots, Random, OrdinaryDiffEq, Statistics
using Distributions, Turing, StatsPlots, DSP
using BSON: @save
pyplot()
Random.seed!(11)

#####
##### Define the controlled dynamical system
#####
include("../scripts/paper_settings.jl")
include(EXAMPLE_FOLDERNAME*"example_lorenz.jl")

# Define the piecewise forcing function
num_jumps = 10;
    U_min = -10;
    U_max = 10
    my_U_func_time = generate_U_func_piecewise(ts, (3, length(ts)),
                                        num_jumps = num_jumps,
                                        U_min=U_min, U_max=U_max,
                                        F_dim=2)
#####
##### Produce data
#####
sol = solve_lorenz_system(my_U_func_time)
dat = Array(sol)
numerical_grad = numerical_derivative(dat, ts)

# Get true derivative
true_grad = zeros(size(dat))
for i in 1:size(dat,2)
    true_grad[:,i] = lorenz_system(true_grad[:,i], dat[:,i], p, [0])
end
# True control signal
U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time(t, dat)
end

plot(U_true')
    title!("True control signal")

#####
##### Calculate distribution of residuals
#####
# Calculate posterior distribution of parameters
chain, train_ind = generate_chain(dat,
                        numerical_grad, lorenz_grad_residual)
# Generate test trajectories from the posterior
vars = [:ρ, :σ, :β]
sample_ind = 1:length(ts)
sample_trajectories, sample_noise = sample_posterior_grad(chain,
            dat, sample_ind, vars, lorenz_system)
# Calculate the residuals
sample_trajectory_mean = transpose(drop_all_1dims(mean(sample_trajectories, dims=1)))
    sample_trajectory_noise = mean(sample_noise)

# Align the signals
dat_grad = numerical_grad[:, sample_ind]

residual = dat_grad .- sample_trajectory_mean
ctr_guess = process_residual(residual, sample_trajectory_noise)

# plot(residual, lw=2, ribbon=mean(sample_noise), fillalpha=0.5,
#             linecolor=COLOR_DICT["intrinsic"], legend=false);
#     title!("Residual", titlefontsize=24)
#     plot!(U_true[save_ind,sample_ind][this_plot_ind], label="True controller")

#####
##### Save everything; plotted in a different script
#####

# Intermediate data
fname = DAT_FOLDERNAME*"dat_piecewise_controller_intermediate.bson";
@save fname chain sol sample_trajectories sample_noise sample_ind numerical_grad

# Original ODE and controller variables
fname = DAT_FOLDERNAME*"dat_piecewise_controller_ode_vars.bson";
@save fname dat dat_grad true_grad U_true

# Bayesian variables
fname = DAT_FOLDERNAME*"dat_piecewise_controller_ctr_vars.bson";
@save fname sample_trajectory_mean sample_trajectory_noise residual ctr_guess


text:
"""
```julia
ExtraState{T,U} <: IterationState{T}
```

Extends `DefaultState` to track extra information each iteration. The field
`other` can be of any type and is simply a storage facility for other data we
want to preserve between iterations. It is not used to check convergence or in
any other place inside IterationManagers.jl. It is simply a convenience for
users who want to keep track of more than just explicit state on each iteration
"""
type ExtraState{T,U} <: IterationState{T}
    default::DefaultState{T}
    other::U
end

ExtraState{T,U}(v::T, other::U) = ExtraState(DefaultState(v), other)

# just forward these methods on to the underlying default state
function Base.show{T,U}(io::IO, es::ExtraState{T,U})
    print(io,
    """
    ExtraState{$T,$U}:
        - default: $(stringmime(MIME"text/plain"(), es.default))
        - other  : $(stringmime(MIME"text/plain"(), es.other))
    """)
end
display_iter(io::IO, es::ExtraState, prefix) =
    display_iter(io, es.default, prefix)
num_iter(es::ExtraState) = num_iter(es.default)
Base.norm(es::ExtraState) = norm(es.default)

# this method calls the default's update! and update the `other` field
function update!{T,U}(es::ExtraState{T,U}, v::T, new_other::U;
                      by::Function=default_by)
    update!(es.default, v)
    es.other = new_other
    nothing
end

# need a special version of managed_iteration!? for this state to handle other
function managed_iteration!{T<:AbstractArray}(f!::Base.Callable,
                                              mgr::IterationManager,
                                              dest::T,
                                              es::ExtraState{T};
                                              by::Base.Callable=default_by)
    pre_hook(mgr, es)

    while !(finished(mgr, es))
        # this function updates dest in place and returns other
        other = f!(dest, es.prev)
        update!(es, dest, other; by=by)
        iter_hook(mgr, es)
    end

    post_hook(mgr, es)
    es
end

function managed_iteration{T}(f::Base.Callable, mgr::IterationManager,
                              es::ExtraState{T};
                              by=default_by)
    pre_hook(mgr, es)

    while !(finished(mgr, es))
        v, other = f(es.default.prev)
        v::T
        update!(es, v, other; by=by)
        iter_hook(mgr, es)
    end

    post_hook(mgr, es)
    es
end


text:
#=
Created on Saturday 3 August 2019
Last update: Thursday 8 August 2019

@author: Michiel Stock
michielfmstock@gmail.com

Some methods to simulate Kronecker graphs.

Leskovec, J., Chakrabarti, D., Kleinberg, J., Faloutsos, C., & Ghahramani,
Z. (2008). Kronecker graphs: an approach to modeling networks.
Journal of Machine Learning Research, 11, 985–1042.
Retrieved from https://cs.stanford.edu/~jure/pubs/kronecker-jmlr10.pdf
=#

using SparseArrays: spzeros
using StatsBase: sample, Weights

"""
    isprob(A::AbstractArray)

Test if a matrix can be interpeted as a probability matrix, i.e., all elements
are between 0 and 1.
"""
isprob(A::AbstractArray) = all(0 .<= A .<= 1)

"""
    isprob(K::AbstractKroneckerProduct)

Test if a Kronecker product can be interpeted as a probability matrix,
i.e., all elements are between 0 and 1.
"""
function isprob(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return isprob(A) && isprob(B)
end

"""
    naivesample(P::AbstractKroneckerProduct)

Sample a Kronecker graph from a probabilistic Kronecker product P using the
naive method. This method has a time complexity in the size of the Kronecker
product (but is still light in memory use). Consider using `fastsample`.
"""
function naivesample(P::AbstractKroneckerProduct)
    @assert isprob(P) throw(DomainError(
                            "All values of K should be between 0 and 1"))
    G = spzeros(Bool, size(P)...)
    for I in CartesianIndices(P)
        if P[I] > rand()  # QUESTION: is this the most efficient way?
            @inbounds G[I] = true
        end
    end
    return G
end

"""
    sampleindices(A::AbstractMatrix, s::Int)

Samples the indices from an `AbstractMatrix`. Probability of sampling indices is
proportional to the size of the corresponding value. Does not do any checks on A.
"""
sampleindices(A::AbstractMatrix, s::Int) = Tuple.(sample(CartesianIndices(A),
                                                    Weights(vec(A), sum(A)), s,
                                                    replace=true))

"""
sampleindices(K::AbstractKroneckerProduct, s::Int)

Samples the indices from an `AbstractKroneckerProduct`. Probability of
sampling indices is proportional to the size of the corresponding value.
Does not do any checks on A.
"""
function sampleindices(K::AbstractKroneckerProduct, s::Int)
    A, B = getmatrices(K)
    p, q = size(B)
    indicesA = sampleindices(A, s)
    indicesB = sampleindices(B, s)
    indices = similar(indicesA)
    for (o, (Ia, Ib)) in enumerate(zip(indicesA, indicesB))
        (i, j), (k, l) = Ia, Ib
        @inbounds indices[o] = ((i - 1) * p + k, (j - 1) * q + l)
    end
    return indices
end

"""
    fastsample(P::AbstractKroneckerProduct)

Uses the heuristic sampling from Leskovec et al. (2008) to sample a large
Kronecker graph.
"""
function fastsample(P::AbstractKroneckerProduct)
    @assert isprob(P) throw(DomainError(
                            "All values of K should be between 0 and 1"))
    G = spzeros(Bool, size(P)...)
    n = Int(round(sum(P)))  # expected number of edges
    for (i, j) in sampleindices(P, n)
        @inbounds G[i,j] = true
    end
    return G
end


text:
# Source: https://www.johndcook.com/blog/2020/02/08/arenstorf-orbit/

module Arenstorf

using Plots

const mmoon = 0.012277471
const mearth = 1-mmoon
const p0 = [0.994, 0]
const v0 = [0, -2.001585106]

function derivative(t, x)
    D1 = ((x[1]+mmoon)*(x[1]+mmoon) + x[2]*x[2])
    D1 *= sqrt(D1)
    D2 = ((x[1]-mearth)*(x[1]-mearth) + x[2]*x[2])
    D2 *= sqrt(D2)
    return [
        x[3],
        x[4],
        x[1] + 2*x[4] - mearth*(x[1]+mmoon)/D1 - mmoon*(x[1]-mearth)/D2,
        x[2] - 2*x[3] - mearth*x[2]/D1 - mmoon*x[2]/D2
    ]
end

function vderivative(t, x)
    r = zeros(size(x))
    for i in 1:size(x, 2)
        r[:,i] = derivative(t, x[:,i])
    end
    return r
end

function init_plot()
    p = plot(
        cos,
        sin,
        0, 2π,
        label = "moon",
        color = :black,
        style = :dot,
        title = "Arenstorf orbit",
        aspect_ratio = :equal,
        framestyle = :box,
        xlims = (-1.4, 1.2),
        legend = :bottomright,
    )
    scatter!((0, 0), label = "", color = :black, markersize = 1.5)
    return p
end

end

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
# SPLITTING DATA SETS

# Helper function for partitioning in the non-stratified case
function _partition(rows, fractions, ::Nothing)
    # container for the row selections (head:tail)
    n_splits = length(fractions) + 1
    heads    = zeros(Int, n_splits)
    tails    = zeros(Int, n_splits)
    n_rows   = length(rows)
    head     = 1
    for (i, p) in enumerate(fractions)
        n = round(Int, p * n_rows)
        iszero(n) && (@warn "A split has only one element."; n = 1)
        # update tail
        tail = head + n - 1
        # store
        heads[i] = head
        tails[i] = tail
        # update head
        head = tail + 1
    end
    if head > n_rows
        @warn "Last vector in the split has only one element."
        head = n_rows
    end
    heads[end] = head
    tails[end] = n_rows
    return tuple((rows[h:t] for (h, t) in zip(heads, tails))...)
end

_make_numerical(v::AbstractVector) =
    throw(ArgumentError("`stratify` must have `Count`, `Continuous` "*
                        "or `Finite` element scitpye. Consider "*
                        "`coerce(stratify, Finite)`. "))
_make_numerical(v::AbstractVector{<:Union{Missing,Real}}) = v
_make_numerical(v::AbstractVector{<:Union{Missing,CategoricalValue}}) =
                                int.(v)

# Helper function for partitioning in the stratified case
function _partition(rows, fractions, raw_stratify::AbstractVector)
    stratify = _make_numerical(raw_stratify)
    length(stratify) == length(rows) ||
        throw(ArgumentError("The stratification vector must "*
                            "have as many entries as " *
                            "the rows to partition."))
    uv    = unique(stratify)
    # construct table (n_classes * idx_of_that_class)
    # NOTE use of '===' is important to handle missing.
    idxs  = [[i for i in rows if stratify[rows[i]] === v] for v in uv]

    # number of occurences of each class and proportions
    nidxs = length.(idxs)
    props = length.(idxs) ./ length(rows)

    n_splits = length(fractions) + 1
    n_rows   = length(rows)

    ns_props = round.(Int, n_rows * fractions * props')
    ns_props = vcat(ns_props, nidxs' .- sum(ns_props, dims=1))

    # warn if anything is >= 1
    if !all(e -> e > 1, ns_props)
        @warn "Some splits have a single or no representative of some class."
    end
    # container for the rows
    split_rows = []
    heads      = ones(Int, length(uv))
    for r in 1:size(ns_props, 1)
        tails   = heads .+ ns_props[r, :] .- 1
        # take chunks of the indices corresponding to the current fraction
        indices = vcat((idxs[i][heads[i]:tails[i]] for i in eachindex(uv))...)
        # rearrange by order of appearance
        indices = sort(indices)
        push!(split_rows, rows[indices])
        heads .= tails .+ 1
    end
    if !all(sl -> sl > 1, length.(split_rows))
        @warn "Some splits have a single or no representative of some class."
    end
    return tuple(split_rows...)
end

const ERR_PARTITION_UNSUPPORTED = ArgumentError(
    "Function `partition` only supports "*
    "AbstractVector, AbstractMatrix or containers implementing the "*
    "Tables interface.")
const ERR_PARTITION_DIMENSION_MISMATCH = DimensionMismatch(
    "Expected a tuple of objects with a common length. ")

__nrows(X) = Tables.istable(X) ? nrows(X) : throw(ERR_PARTITION_UNSUPPORTED)
__nrows(X::Union{AbstractMatrix,AbstractVector}) = nrows(X)

"""
    partition(X, fractions...;
              shuffle=nothing,
              rng=Random.GLOBAL_RNG,
              stratify=nothing,
              multi=false)

Splits the vector, matrix or table `X` into a tuple of objects of the
same type, whose vertical concatenation is `X`. The number of rows in
each component of the return value is determined by the
corresponding `fractions` of `length(nrows(X))`, where valid fractions
are floats between 0 and 1 whose sum is less than one. The last
fraction is not provided, as it is inferred from the preceding ones.

For "synchronized" partitioning of multiple objects, use the
`multi=true` option described below.

    julia> partition(1:1000, 0.8)
    ([1,...,800], [801,...,1000])

    julia> partition(1:1000, 0.2, 0.7)
    ([1,...,200], [201,...,900], [901,...,1000])

    julia> partition(reshape(1:10, 5, 2), 0.2, 0.4)
    ([1 6], [2 7; 3 8], [4 9; 5 10])

    X, y = make_blobs() # a table and vector
    Xtrain, Xtest = partition(X, 0.8, stratify=y)

    (Xtrain, Xtest), (ytrain, ytest) = partition((X, y), 0.8, rng=123, multi=true)

## Keywords

* `shuffle=nothing`: if set to `true`, shuffles the rows before taking
  fractions.

* `rng=Random.GLOBAL_RNG`: specifies the random number generator to be
  used, can be an integer seed. If specified, and `shuffle ===
  nothing` is interpreted as true.

* `stratify=nothing`: if a vector is specified, the partition will
  match the stratification of the given vector. In that case,
  `shuffle` cannot be `false`.

* `multi=false`: if `true` then `X` is expected to be a `tuple` of
  objects sharing a common length, which are each partitioned
  separately using the same specified `fractions` *and* the same row
  shuffling. Returns a tuple of partitions (a tuple of tuples).

"""
function partition(X, fractions::Real...;
                   shuffle::Union{Nothing,Bool}=nothing,
                   rng=Random.GLOBAL_RNG,
                   stratify::Union{Nothing,AbstractVector}=nothing,
                   multi=false)

    # check the fractions
    if !all(e -> 0 < e < 1, fractions) || sum(fractions) >= 1
        throw(DomainError(fractions,
                          "Fractions must be in (0, 1) with sum < 1."))
    end

    # determinen `n_rows`:
    if X isa Tuple && multi
        isempty(X) && return tuple(fill((), length(fractions) + 1)...)
        x = first(X)
        n_rows = __nrows(x)
        all(X[2:end]) do x
            nrows(x) === n_rows
        end || throw(ERR_PARTITION_DIMENSION_MISMATCH)
    else
        n_rows = __nrows(X)
    end

    # check the rng & adjust shuffling
    if rng isa Integer
        rng = MersenneTwister(rng)
    end
    if rng != Random.GLOBAL_RNG && shuffle === nothing
        shuffle = true
    end
    rows = collect(1:n_rows)
    shuffle !== nothing && shuffle && shuffle!(rng, rows)

    # determine the partition of `rows`:
    row_partition = _partition(rows, collect(fractions), stratify)
    return _partition(X, row_partition)
end

function _partition(X, row_partition)
#    _X = Tables.istable(X) ? X : collect(X)
    return tuple((selectrows(X, p) for p in row_partition)...)
end

_partition(X::Tuple, row_partition) =
    map(x->_partition(x, row_partition), X)

# # UNPACK

"""
    unpack(table, f1, f2, ... fk;
           wrap_singles=false,
           shuffle=false,
           rng::Union{AbstractRNG,Int,Nothing}=nothing,
           coerce_options...)

Horizontally split any Tables.jl compatible `table` into smaller
tables or vectors by making column selections determined by the
predicates `f1`, `f2`, ..., `fk`. Selection from the column names is
without replacement. A *predicate* is any object `f` such that
`f(name)` is `true` or `false` for each column `name::Symbol` of
`table`.

Returns a tuple of tables/vectors with length one greater than the
number of supplied predicates, with the last component including all
previously unselected columns.

```
julia> table = DataFrame(x=[1,2], y=['a', 'b'], z=[10.0, 20.0], w=["A", "B"])
2×4 DataFrame
 Row │ x      y     z        w
     │ Int64  Char  Float64  String
─────┼──────────────────────────────
   1 │     1  a        10.0  A
   2 │     2  b        20.0  B

Z, XY, W = unpack(table, ==(:z), !=(:w))
julia> Z
2-element Vector{Float64}:
 10.0
 20.0

julia> XY
2×2 DataFrame
 Row │ x      y
     │ Int64  Char
─────┼─────────────
   1 │     1  a
   2 │     2  b

julia> W  # the column(s) left over
2-element Vector{String}:
 "A"
 "B"
```

Whenever a returned table contains a single column, it is converted to
a vector unless `wrap_singles=true`.

If `coerce_options` are specified then `table` is first replaced
with `coerce(table, coerce_options)`. See
[`ScientificTypes.coerce`](@ref) for details.

If `shuffle=true` then the rows of `table` are first shuffled, using
the global RNG, unless `rng` is specified; if `rng` is an integer, it
specifies the seed of an automatically generated Mersenne twister. If
`rng` is specified then `shuffle=true` is implicit.

"""
function unpack(X, predicates...;
                wrap_singles=false,
                shuffle=nothing,
                rng=nothing, pairs...)

    # add a final predicate to unpack all remaining columns into to
    # the last return value:
    predicates = (predicates..., _ -> true)

    shuffle, rng = shuffle_and_rng(shuffle, rng)

    shuffle && (X = selectrows(X, Random.shuffle(rng, 1:nrows(X))))

    if isempty(pairs)
        Xfixed = X
    else
        Xfixed = coerce(X, pairs...)
    end

    unpacked = Any[]
    names_left = schema(Xfixed).names |> collect

    for c in predicates
        names = filter(c, names_left)
        filter!(!in(names), names_left)
        length(names) == 1 && !wrap_singles && (names = names[1])
        push!(unpacked, selectcols(Xfixed, names))
    end
    return Tuple(unpacked)
end

## RESTRICTING TO A FOLD

struct FoldRestrictor{i,N}
    f::NTuple{N,Vector{Int}}
end
(r::FoldRestrictor{i})(X) where i = selectrows(X, (r.f)[i])

"""
    restrict(X, folds, i)

The restriction of `X`, a vector, matrix or table, to the `i`th fold
of `folds`, where `folds` is a tuple of vectors of row indices.

The method is curried, so that `restrict(folds, i)` is the operator
on data defined by `restrict(folds, i)(X) = restrict(X, folds, i)`.

### Example

    folds = ([1, 2], [3, 4, 5],  [6,])
    restrict([:x1, :x2, :x3, :x4, :x5, :x6], folds, 2) # [:x3, :x4, :x5]

See also [`corestrict`](@ref)

"""
restrict(f::NTuple{N}, i) where N = FoldRestrictor{i,N}(f)
restrict(X, f, i) = restrict(f, i)(X)

## RESTRICTING TO A FOLD COMPLEMENT

"""
    complement(folds, i)

The complement of the `i`th fold of `folds` in the concatenation of
all elements of `folds`. Here `folds` is a vector or tuple of integer
vectors, typically representing row indices or a vector, matrix or
table.

    complement(([1,2], [3,], [4, 5]), 2) # [1 ,2, 4, 5]

"""
complement(f, i) = reduce(vcat, collect(f)[Not(i)])

struct FoldComplementRestrictor{i,N}
    f::NTuple{N,Vector{Int}}
end
(r::FoldComplementRestrictor{i})(X) where i =
    selectrows(X, complement(r.f, i))

"""
    corestrict(X, folds, i)

The restriction of `X`, a vector, matrix or table, to the *complement*
of the `i`th fold of `folds`, where `folds` is a tuple of vectors of
row indices.

The method is curried, so that `corestrict(folds, i)` is the operator
on data defined by `corestrict(folds, i)(X) = corestrict(X, folds, i)`.

### Example

    folds = ([1, 2], [3, 4, 5],  [6,])
    corestrict([:x1, :x2, :x3, :x4, :x5, :x6], folds, 2) # [:x1, :x2, :x6]

"""
corestrict(f::NTuple{N}, i) where N = FoldComplementRestrictor{i,N}(f)
corestrict(X, f, i) = corestrict(f, i)(X)

## to be replaced (not used anywhere):
## ACCESSORS FOR JULIA NDSPARSE ARRAYS (N=2)

# nrows(::Val{:sparse}, X) = maximum([r[1] for r in keys(X)])
# function select(::Val{:sparse}, X, r::Integer, c::Symbol)
#     try
#         X[r,c][1]
#     catch exception
#         exception isa KeyError || throw(exception)
#         missing
#     end
# end
# select(::Val{:sparse}, X, r::AbstractVector{<:Integer}, c::Symbol) = [select(X, s, c) for s in r]
# select(::Val{:sparse}, X, ::Colon, c::Symbol) = [select(X, s, c) for s in 1:nrows(X)]
# selectrows(::Val{:sparse}, X, r::Integer) = X[r:r,:]
# selectrows(::Val{:sparse}, X, r) = X[r,:]
# selectcols(::Val{:sparse}, X, c::Symbol) = select(X, :, c)
# selectcols(::Val{:sparse}, X, c::AbstractVector{Symbol}) = X[:,sort(c)]
# selectcols(::Val{:sparse}, X, ::Colon) = X
# select(::Val{:sparse}, X, r::Integer, c::AbstractVector{Symbol}) = X[r,sort(c)]
# select(::Val{:sparse}, X, r::Integer, ::Colon) = X[r,:]
# select(::Val{:sparse}, X, r, c) = X[r,sort(c)]

## TRANSFORMING BETWEEN CATEGORICAL ELEMENTS AND RAW VALUES

MLJModelInterface.transform(
    e::Union{CategoricalArray,CategoricalValue,CategoricalPool},
    arg) = CategoricalDistributions.transform(e, arg)

## SKIPPING MISSING AND NAN: skipinvalid

_isnan(x) = false
_isnan(x::Number) = isnan(x)

skipnan(x) = Iterators.filter(!_isnan, x)

"""
    skipinvalid(itr)

Return an iterator over the elements in `itr` skipping `missing` and
`NaN` values. Behaviour is similar to [`skipmissing`](@ref).

    skipinvalid(A, B)

For vectors `A` and `B` of the same length, return a tuple of vectors
`(A[mask], B[mask])` where `mask[i]` is `true` if and only if `A[i]`
and `B[i]` are both valid (non-`missing` and non-`NaN`). Can also
called on other iterators of matching length, such as arrays, but
always returns a vector. Does not remove `Missing` from the element
types if present in the original iterators.

"""
skipinvalid(v) = v |> skipmissing |> skipnan

isinvalid(x) = ismissing(x) || _isnan(x)

function skipinvalid(yhat, y)
    mask = .!(isinvalid.(yhat) .| isinvalid.(y))
    return yhat[mask], y[mask]
end

# TODO: refactor balanced accuracy to get rid of these:

function _skipinvalid(yhat, y, w::Arr)
    mask = .!(isinvalid.(yhat) .| isinvalid.(y))
    return yhat[mask], y[mask], w[mask]
end

function _skipinvalid(yhat, y, w::Union{Nothing,AbstractDict})
    mask = .!(isinvalid.(yhat) .| isinvalid.(y))
    return yhat[mask], y[mask], w
end


text:
import GuillotineModels.Utilities: save_and_relax!, save_and_fix!, restore!
import JuMP
import JuMP: @variable, @constraint, @objective
import GLPK
import MathOptInterface
const MOI = MathOptInterface
import Suppressor: @suppress_out

function test_relax_continuous()
	# Not much to do, just solve an LP and check if the solution keeps being
	# the same after the relax.
	m = JuMP.Model(GLPK.Optimizer)
	@variable(m, 0.5 <= c <= 2)
	@variable(m, -2 <= d <= 0.5)
	# Get the optimal value of the maximization LP.
	@objective(m, Max, c + d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6
	# Get the optimal value of the maximization LP (relaxed).
	svc_c = save_and_relax!(c)
	svc_d = save_and_relax!(d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6
	# Get the optimal value of the maximization LP (restored).
	restore!.([c, d], [svc_c, svc_d])
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6
end

function test_relax_integer()
	m = JuMP.Model(GLPK.Optimizer)
	@variable(m, 0.5 <= c <= 2, Int)
	@variable(m, -2 <= d <= 0.5, Int)
	# Get the optimal value of the maximization MIP.
	@objective(m, Max, c + d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.0 rtol=1e-6 atol=1e-6
	# Get the optimal value of the maximization LP.
	svc_c = save_and_relax!(c)
	svc_d = save_and_relax!(d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6
	# Get the optimal value of the minimization LP.
	@objective(m, Min, c + d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ -1.5 rtol=1e-6 atol=1e-6
	# Get the optimal value of the restored MIP (but for minimization now).
	restore!.([c, d], [svc_c, svc_d])
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ -1.0 rtol=1e-6 atol=1e-6
end

function test_relax_binary()
	# Given a model with three binary variables, it gives a 2.0 objective value
	# because two variables may only assume value zero.
	m = JuMP.Model(GLPK.Optimizer)
	@variable(m, a, Bin)
	@variable(m, b, Bin)
	@variable(m, 0.5 <= c <= 2, Bin)
	@variable(m, -2 <= d <= 0.5, Bin)
	@objective(m, Max, a + b + c + d)
	@constraint(m, b <= 0.5)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.0 rtol=1e-6 atol=1e-6

	# If variable `a` is relaxed, the implicit bound at 1 should yet
	# be respected.
	svc_a = save_and_relax!(a)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.0 rtol=1e-6 atol=1e-6

	# If variable `b` is relaxed, the constraint should keep it at 0.5 anyway.
	svc_b = save_and_relax!(b)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6

	# If variable `c` is relaxed, the implicit bound at 1 should be respected,
	# even if there is a explicit bound at 2.
	svc_c = save_and_relax!(c)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.5 rtol=1e-6 atol=1e-6

	# If variable `d` is relaxed, the implicit bound at 0.5 should be respected.
	svc_d = save_and_relax!(d)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 3.0 rtol=1e-6 atol=1e-6

	# If the variables are restores to their original state, the same initial
	# objective value should reappear.
	restore!.([a, b, c, d], [svc_a, svc_b, svc_c, svc_d])
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 2.0 rtol=1e-6 atol=1e-6

	# If the problem minimizes, the implicit bound at c should prevent it
	# becoming zero (instead it stops at 0.5). None should dip below zero.
	@objective(m, Min, a + b + c + d)
	save_and_relax!.([a, b, c, d])
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 0.5 rtol=1e-6 atol=1e-6
	@test JuMP.value(a) ≈ 0.0 rtol=1e-6 atol=1e-6
	@test JuMP.value(b) ≈ 0.0 rtol=1e-6 atol=1e-6
	@test JuMP.value(c) ≈ 0.5 rtol=1e-6 atol=1e-6
	@test JuMP.value(d) ≈ 0.0 rtol=1e-6 atol=1e-6
end

function test_relax()
	@testset "Relax Binary" begin test_relax_binary() end
	@testset "Relax Integer" begin test_relax_integer() end
	@testset "Relax Continuous" begin test_relax_continuous() end
end

function test_fix()
	m = JuMP.Model(GLPK.Optimizer)
	@variable(m, 10 <= a <= 20)
	@variable(m, b, Bin)
	@variable(m, 30 <= c <= 40, Int)
	@objective(m, Max, a + b + c)
	@constraint(m, b >= 0.1)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 61.0 rtol=1e-6 atol=1e-6
	# Fix a continuous variable to a valid value.
	svc_a = save_and_fix!(a, 15.5)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 56.5 rtol=1e-6 atol=1e-6
	# Fix the binary variable to a value that would be valid by its binary
	# bounds, but there is a constraint in conflict.
	svc_b = save_and_fix!(b, 0.0)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.NO_SOLUTION
	# Restore the binary variable.
	restore!(b, svc_b)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 56.5 rtol=1e-6 atol=1e-6
	# Fix an integer variable to a valid value.
	svc_c = save_and_fix!(c, 35)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 51.5 rtol=1e-6 atol=1e-6
	# Fix a integer variable to a value that is not integer.
	svc_c2 = save_and_fix!(c, 35.5)
	@suppress_out JuMP.optimize!(m) # GLPK outputs some warnings if not suppressed
	@test JuMP.primal_status(m) == MOI.NO_SOLUTION
	# Restore the `c` variable from its first save_and_fix! not the last.
	restore!(c, svc_c)
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ 56.5 rtol=1e-6 atol=1e-6
end

test_relax()
test_fix()


text:
"""
    neighbor_fpidx(IDs, ps, fpg)

global id
local rank

"""
function neighbor_fpidx(IDs, ps, fpg)
    # id-th cell, fd-th face, jd-th point
    id, fd, jd = IDs

    # ending point ids of a face
    if fd == 1
        pids = [ps.cellid[id, 1], ps.cellid[id, 2]]
    elseif fd == 2
        pids = [ps.cellid[id, 2], ps.cellid[id, 3]]
    elseif fd == 3
        pids = [ps.cellid[id, 3], ps.cellid[id, 1]]
    end

    # global face index
    faceids = ps.cellFaces[id, :]

    function get_faceid()
        for i in eachindex(faceids)
            if sort(pids) == sort(ps.facePoints[faceids[i], :])
                return faceids[i]
            end
        end

        @warn "no face id found"
    end
    faceid = get_faceid()

    # neighbor cell id
    neighbor_cid = setdiff(ps.faceCells[faceid, :], id)[1]

    # in case of boundary cell
    if neighbor_cid <= 0
        return neighbor_cid, -1, -1
    end

    # face rank in neighbor cell
    if ps.cellid[neighbor_cid, 1] ∉ ps.facePoints[faceid, :]
        neighbor_frk = 2
    elseif ps.cellid[neighbor_cid, 2] ∉ ps.facePoints[faceid, :]
        neighbor_frk = 3
    elseif ps.cellid[neighbor_cid, 3] ∉ ps.facePoints[faceid, :]
        neighbor_frk = 1
    end

    # point rank in neighbor cell
    neighbor_nrk1 =
        findall(x -> x == fpg[id, fd, jd, 1], fpg[neighbor_cid, neighbor_frk, :, 1])
    neighbor_nrk2 =
        findall(x -> x == fpg[id, fd, jd, 2], fpg[neighbor_cid, neighbor_frk, :, 2])
    neighbor_nrk = intersect(neighbor_nrk1, neighbor_nrk2)[1]

    return neighbor_cid, neighbor_frk, neighbor_nrk
end


text:
using RobotZoo
using RobotDynamics
using Test
using StaticArrays
using BenchmarkTools
using Rotations

@testset "allocations" begin
    include("models_test.jl")
end

@testset "LinearModels" begin
    include("linear_models_test.jl")
end


text:
abstract type AccessMode end
abstract type RecordMarkerType end
abstract type Record<:IO end

struct SequentialAccess{RT<:RecordMarkerType} <: AccessMode
   recmrktyp :: RT
end

Base.show(io::IO, a::SequentialAccess) =
   print(io, "sequential-access, ", a.recmrktyp)

struct DirectAccess <: AccessMode
   reclen :: Int64
end

Base.show(io::IO, a::DirectAccess) =
   print(io, "direct-access, $(a.reclen)-byte records")

struct WithoutSubrecords{T} <: RecordMarkerType
end

const RECMRK8B = WithoutSubrecords{Int64}()

Base.show(io::IO, ::WithoutSubrecords{T}) where {T} =
   print(io, "$(sizeof(T))-byte record markers, no subrecords")

struct WithSubrecords <: RecordMarkerType
   max_subrecord_length :: Int32
end

const max_subrecord_length = 2^31-9
const RECMRK4B = WithSubrecords(max_subrecord_length)

const RECMRKDEF = RECMRK4B

Base.show(io::IO, rt::WithSubrecords) =
   print(io, "4-byte record markers, subrecords of max $(rt.max_subrecord_length) bytes")

struct Conversion{R,W}
   onread  :: R
   onwrite :: W
   name    :: String
end

const converts = [
   Conversion( identity, identity, "native"        ),
   Conversion( ntoh    , hton    , "big-endian"    ),
   Conversion( ltoh    , htol    , "little-endian" )
   ]

const NOCONV = typeof(converts[1])

function get_convert(name::String)
   for conv in converts
      if conv.name == name
         return conv
      end
   end
   fthrow("unknown convert method \"$name\"")
end


text:
using TaskGraphs

## set up the environment
using GraphUtils

indicator_grid = IndicatorGrid(zeros(Bool,11,5))                           # initialize a 10x10 environment with no obstacles
# indicator_grid[4,4] = true                                                                # obstacle at position 4,4
indicator_grid[1:9,2] .= true
indicator_grid[1:9,4] .= true                                                          # obstacle covering positions 6:7,6:7
env = construct_factory_env_from_indicator_grid(indicator_grid) # initialize the factory environment
@show env.vtx_map                                                                     # visualize the mapping from grid position to vertex indices
# exit()
## set up the environment
# vtx_grid = initialize_dense_vtx_grid(7,5) # 4 x 4 grid world
#  1   2   3   4	5
#  5   6   7   8
#  9  10  11  12
# 13  14  15  16
# env = construct_factory_env_from_vtx_grid(vtx_grid)
# add_obstacle!(env)
## Define the initial conditions of the robots
robot_ics = [
    ROBOT_AT(1,27), # robot 1 starts at vertex 2
    ROBOT_AT(2,36), # robot 2 starts at vertex 9
    ROBOT_AT(3,29),
]

## Define the manufacturing project
spec = ProjectSpec()
# set initial conditions of "raw materials"
set_initial_condition!(spec,OBJECT_AT(1,27))
set_initial_condition!(spec,OBJECT_AT(2,36))
set_initial_condition!(spec,OBJECT_AT(3,29))

# set_initial_condition!(spec,OBJECT_AT(5,7))
# set_initial_condition!(spec,OBJECT_AT(6,2))
# set_initial_condition!(spec,OBJECT_AT(7,12))

# define the operations that need to take place
# define the operations that need to take place

op1 = Operation(Δt=0)
op2 = Operation(Δt=0)
op3 = Operation(Δt=0)

# inputs
set_precondition!(op1,OBJECT_AT(1,26))
set_precondition!(op2,OBJECT_AT(2,20))
set_precondition!(op3,OBJECT_AT(3,14))

# outputs
set_postcondition!(op1,OBJECT_AT(4,23))
set_postcondition!(op2,OBJECT_AT(5,17))
set_postcondition!(op3,OBJECT_AT(6,11))

add_operation!(spec,op1)
add_operation!(spec,op2)
add_operation!(spec,op3)

# add a terminal operation with no outputs
op4 = Operation(Δt=0)
op5 = Operation(Δt=0)
op6 = Operation(Δt=0)

# outputs
set_precondition!(op4,OBJECT_AT(4,2))
set_precondition!(op5,OBJECT_AT(5,5))
set_precondition!(op6,OBJECT_AT(6,8))

add_operation!(spec,op4)
add_operation!(spec,op5)
add_operation!(spec,op6)

## define solver
solver = NBSSolver()
# finalize problem construction (the solver is passed as an argument here
# because it determines what cost model is used for the problem)
prob = pctapf_problem(solver,spec,env,robot_ics)
# solve the problem
solution, cost = solve!(solver,prob)
println(env.vtx_map)
println(solution)
println(cost)
# check if the problem was solved to optimality
@show feasible_status(solver)
@show optimal_status(solver)


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module RLESCAS

export include_visualize
using Reexport

const DIR = dirname(@__FILE__)

include("config/config_ACASX_GM.jl") #defineSim
@reexport using .Config_ACASX_GM

#Config AdaptiveStressTest
include("config/config_ast.jl") #defineAST
@reexport using .ConfigAST

#Config MCTS solver
include("config/config_mcts.jl") #defineMCTS
@reexport using .ConfigMCTS

#Config MCBest solver
include("defines/mcbest.jl") #defineSim
@reexport using .MCBest
include("config/config_mcbest.jl") #defineMCBest
@reexport using .ConfigMCBest

include("defines/define_save.jl") #trajSave, trajLoad and helpers
@reexport using .DefineSave
include("defines/define_log.jl") #SimLog
@reexport using .DefineLog
include("defines/save_types.jl") #ComputeInfo
@reexport using .SaveTypes
include("helpers/save_helpers.jl")
@reexport using .SaveHelpers

include("helpers/label270.jl")
@reexport using .Label270
include("helpers/add_supplementary.jl") #add label270
@reexport using .AddSupplementary
include("tools/label270_to_text.jl")
@reexport using .Label270_To_Text
include("tools/summarize.jl")
@reexport using .Summarize
include("converters/log_to_csv.jl")
@reexport using .Log_To_CSV
include("converters/log_to_scripted.jl")
@reexport using .Log_To_Scripted
include("converters/log_to_waypoints.jl")
@reexport using .Log_To_Waypoints

include("helpers/plot_nmacs.jl")
@reexport using .PlotNMACs

include("visualize/visualize.jl")

function include_visualize()
  @eval @reexport using .Visualize
end

include("trajsave/trajSave_common.jl")
@reexport using .TrajSaveCommon
include("helpers/fill_to_max_time.jl")
@reexport using .Fill_To_Max_Time
include("trajsave/trajSave_replay.jl")
@reexport using .TrajSaveReplay
include("trajsave/post_process.jl")
@reexport using .PostProcess
include("trajsave/trajSave_once.jl")
@reexport using .TrajSaveOnce
include("trajsave/trajSave_mcbest.jl")
@reexport using .TrajSaveMCBest
include("trajsave/trajSave_mcts.jl")
@reexport using .TrajSaveMCTS

include("tools/nmac_stats.jl")
@reexport using .NMACStats
#include("tools/json_update.jl")
#@reexport using .JsonUpdate
include("tools/recursive_plot.jl")
@reexport using .RecursivePlot
include("tools/collage.jl")
@reexport using .Collage
include("tools/sortbylikelihood.jl")
@reexport using .SortByLikelihood

include("metrics/metrics.jl")
@reexport using .Metrics

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
module HelFEM
using ReadOnlyArrays
using CxxWrap: CxxPtr

export PolynomialBasis, RadialBasis, FEMBasis,
    boundaries, nelements, radial_integral

include("helfem_cxxwrap.jl")
include("polynomialbasis.jl")
include("radialbasis.jl")
include("fembasis.jl")

invh(S) = helfem.invh(S, false)

## Quadrature

"""
    chebyshev(nquad) -> (xs::Vector, ws::Vector)

Return HelFEM's modified Gauss-Chebyshev quadrature points and weights for numerical
integration of functions on ``x \\in [-1, 1]``.
"""
function chebyshev(nquad)
    xs, ws = HelFEM.helfem.ArmaVector(), HelFEM.helfem.ArmaVector()
    HelFEM.helfem.chebyshev(nquad, xs, ws)
    (xs = collect(xs), ws = collect(ws))
end

"""
    lobatto(nquad) -> (xs::Vector, ws::Vector)

Return HelFEM's Gauss-Lobatto quadrature points and weights for numerical on
``x \\in [-1, 1]``.
"""
function lobatto(nquad)
    xs, ws = HelFEM.helfem.ArmaVector(), HelFEM.helfem.ArmaVector()
    HelFEM.helfem.lobatto(nquad, xs, ws)
    (xs = collect(xs), ws = collect(ws))
end

end # module


text:
"""
    module AtomicMod

Simple module for atom type.
"""
module AtomicMod

using ..Utility:str_w_spaces
using ..BandTools:band_energy
using ..BandTools:smearing_energy

"""
    struct atom

Hold basic atomic information

- `name::String`
- `Z::Int64` Atomic number
- `row::Float64` Periodic table row
- `col::Float64` Periodic table col
- `mass::Float64` Mass in amu
- `nval::Float64` Number of valence electrons in TB calculation.
- `nsemicore::Int64` Number of semicore electrons. Depends on pseudopotential choice
- `nwan::Int64` Number of TB orbitals
- `orbitals::Array{Symbol,1}` Names of the orbitals, like `[:s, :p]`
- `total_energy::Float64` DFT total energy, depends on pseudopotentials.
- `eigs::Dict` orbital eigenvalues of isolated non-spin-polarized nuetral atom.
- `energy_offset::Float64` energy to add to TB calculation to make isolated atoms have zero energy.
- `U::Float64` U value for Ewald correction

"""
struct atom
    name::String
    Z::Int64
    row::Float64
    col::Float64
    mass::Float64
    nval::Float64
    nsemicore::Int64
    nwan::Int64
    orbitals::Array{Symbol,1}
    total_energy::Float64
    eigs::Dict
    energy_offset::Float64
    U::Float64
end

#printing
Base.show(io::IO, a::atom) = begin

    println(io,str_w_spaces([a.name, "Z=", a.Z, "Mass=",a.mass, "RC=",a.row, a.col]))
    println(io,str_w_spaces(["Nval=", a.nval, "Nsemicore=", a.nsemicore, "Orbitals=",a.orbitals]))
    println(io,"U = ",a.U)
    println(io,"total energy=", a.total_energy)
    for i in a.orbitals
        println(i, "  " , a.eigs[i])
    end
end

"""
    function makeatom(name, Z, row, col, mass, nval, nsemicore, orbitals, etot, eigs, vac_potential=0.0, U=0.0)

Constructor for atom.
"""
function makeatom(name, Z, row, col, mass, nval, nsemicore, orbitals, etot, eigs, vac_potential=0.0, U=0.0)
    #vac potential in ryd

    orb2 = map(x->convert(Symbol, x), orbitals)

    nwan = 0

    for o in orb2
        if o == :s
            nwan += 2
        elseif o == :p
            nwan += 6
        elseif o == :d
            nwan += 10
        elseif o == :f
            nwan += 14
        else
            exit("oribtal entry incorrect, good values are s p d f: ", orbitals, orb2)
        end
    end

    convert_ev_ryd = 1.0/13.605693122
    EIGS = zeros(Int64(nwan/2), 1)
    c=0
    for (i,o) in enumerate(orb2)
        if o == :s
            c+=1
            EIGS[c,1] = eigs[i]*convert_ev_ryd
        elseif o == :p
            c += 3
            EIGS[c-2,1] = eigs[i]*convert_ev_ryd
            EIGS[c-1,1] = eigs[i]*convert_ev_ryd
            EIGS[c  ,1] = eigs[i]*convert_ev_ryd
        elseif o == :d
            c += 5
#            EIGS[c-5,1] = eigs[i]*convert_ev_ryd
            EIGS[c-4,1] = eigs[i]*convert_ev_ryd
            EIGS[c-3,1] = eigs[i]*convert_ev_ryd
            EIGS[c-2,1] = eigs[i]*convert_ev_ryd
            EIGS[c-1,1] = eigs[i]*convert_ev_ryd
            EIGS[c  ,1] = eigs[i]*convert_ev_ryd
        elseif o == :f
            c += 7
#            EIGS[c-7,1] = eigs[i]*convert_ev_ryd
            EIGS[c-6,1] = eigs[i]*convert_ev_ryd
            EIGS[c-5,1] = eigs[i]*convert_ev_ryd
            EIGS[c-4,1] = eigs[i]*convert_ev_ryd
            EIGS[c-3,1] = eigs[i]*convert_ev_ryd
            EIGS[c-2,1] = eigs[i]*convert_ev_ryd
            EIGS[c-1,1] = eigs[i]*convert_ev_ryd
            EIGS[c  ,1] = eigs[i]*convert_ev_ryd
        end
    end

    EIGS = EIGS .- vac_potential

#    println("EIGS after subtraction")
#    println(EIGS)

    band_energy0, efermi = band_energy(EIGS, [1.0], nval, returnef=true)
    smear_energy = smearing_energy(EIGS, [1.0], efermi)

#    if name=="Ta"
#        println("name $name")
#        println("EIGS before subtraction")
#        println(EIGS)
#        println("band_energy0 ", band_energy0, " ", efermi)
#    end

#    shift = (0.0 - band_energy0)/nval

    energy_offset = -(band_energy0+smear_energy)  #this is the atomic energy, no spin, after vac correction

    if name == "X" || name == "Xa"
        energy_offset = 0.0
    end

#    println("$name energy_offset $energy_offset    $band_energy0 $smear_energy  efermi $efermi")

    d = Dict()

    for (o, i) in zip(orb2, eigs)
#        d[o] = (i * convert_ev_ryd + shift)
        d[o] = (i * convert_ev_ryd - vac_potential)

#        println("d ", d[o], " ", o)
    end

#    band_energy_new = band_energy(EIGS .- band_energy0 / nval , [1.0], nval)#
#
#    println("Loading $name, new band energy $band_energy_new, ", d[:s])

    return atom(name, Z, row, col, mass, nval, nsemicore, nwan, orb2, etot, d, energy_offset, U)

end

end #end module


text:
"""
Define your agents to be a subtype of AbstractAgent. Your agent type has to have the following fields: `id`, `pos`

e.g.

```
mutable struct MyAgent <: AbstractAgent
  id::Integer
  pos::Tuple{Integer, Integer, Integer}
end
```

Agents should have an `id` and a `pos` (position) field. If your space is a grid, the position should accept a `Tuple{Integer, Integer, Integer}` or a `Tuple{Integer, Integer}` representing x, y, and optionally z coordinates.
"""
abstract type AbstractAgent end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the precession computation of IAU-2006.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Wallace, P. T., Capitaine, N (2006). Precession-nutation procedures
#       consistent with IAU 2006 resolutions. Astronomy & Astrophysics.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/transformations/iau2006/precession.jl
# =================================================

# Function precession_iau2006
# ---------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Numerical example at Appendix of [1]
#
# According to this example, using:
#
#   UTC    = 2006 January 15, 21h 24.375
#   JD_TT  = 2400000.5+53750.892855138888889
#   JD_UT1 = 2400000.5+53750.892104561342593
#
# one gets:
#
#   Ψ_a = +304.359364139   arcsec
#   ω_a = +84381.404629617 arcsec
#   χ_a = +0.628998164     arcsec
#
################################################################################

@testset "Function precession_iau2006" begin
    JD_TT  = 2400000.5+53750.892855138888889

    r2a = 180/π*3600

    Ψ_a, ω_a, χ_a = precession_iau2006(JD_TT)

    @test Ψ_a*r2a ≈   +304.359364139 atol = 1e-9
    @test ω_a*r2a ≈ +84381.404629617 atol = 1e-9
    @test χ_a*r2a ≈     +0.628998164 atol = 1e-9
end


text:
# Tests for our auditing infrastructure

@testset "Auditor - cppfilt" begin
    # We take some known x86_64-linux-gnu symbols and pass them through c++filt
    mangled_symbol_names = [
        "_ZNKSt7__cxx1110_List_baseIiSaIiEE13_M_node_countEv",
        "_ZNKSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEE6lengthEv@@GLIBCXX_3.4.21",
        "_Z10my_listlenNSt7__cxx114listIiSaIiEEE",
        "_ZNKSt7__cxx114listIiSaIiEE4sizeEv",
    ]
    unmangled_symbol_names = BinaryBuilder.cppfilt(mangled_symbol_names, Linux(:x86_64))
    @test all(unmangled_symbol_names .== [
        "std::__cxx11::_List_base<int, std::allocator<int> >::_M_node_count() const",
        "std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >::length() const@@GLIBCXX_3.4.21",
        "my_listlen(std::__cxx11::list<int, std::allocator<int> >)",
        "std::__cxx11::list<int, std::allocator<int> >::size() const",
    ])
end

@testset "Auditor - ISA tests" begin
    mktempdir() do build_path
        products = Product[
            ExecutableProduct("main_sse", :main_sse),
            ExecutableProduct("main_avx", :main_avx),
            ExecutableProduct("main_avx2", :main_avx2),
        ]

        build_output_meta = nothing
        @test_logs (:warn, r"sandybridge") (:warn, r"haswell") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "isa_tests",
                v"1.0.0",
                [DirectorySource(build_tests_dir)],
                # Build the test suite, install the binaries into our prefix's `bin`
                raw"""
                cd ${WORKSPACE}/srcdir/isa_tests
                make -j${nproc} install
                install_license /usr/include/ltdl.h
                """,
                # Build for our platform
                [platform],
                # Ensure our executable products are built
                products,
                # No dependencies
                Dependency[];
                # We need to build with very recent GCC so that we can emit AVX2
                preferred_gcc_version=v"8",
            )
        end

        # Extract our platform's build
        @test haskey(build_output_meta, platform)
        tarball_path, tarball_hash = build_output_meta[platform][1:2]
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        prefix = Prefix(testdir)

        # Run ISA tests
        for (product, true_isa) in zip(products, (:core2, :sandybridge, :haswell))
            readmeta(locate(product, prefix)) do oh
                detected_isa = BinaryBuilder.analyze_instruction_set(oh, platform; verbose=true)
                @test detected_isa == true_isa
            end
        end
    end
end

@testset "Auditor - cxxabi selection" begin
    libcxxstringabi_test = LibraryProduct("libcxxstringabi_test", :libcxxstringabi_test)

    # Factor the autobuild() out
    function do_build(build_path, script, platform, gcc_version)
        autobuild(
            build_path,
            "libcxxstringabi_test",
            v"1.0.0",
            # Copy in the build_tests sources
            [DirectorySource(build_tests_dir)],
            script,
            # Build for this platform
            [platform],
            # The products we expect to be build
            [libcxxstringabi_test],
            # No depenedencies
            Dependency[];
            preferred_gcc_version=gcc_version
        )
    end

    for platform in (Linux(:x86_64; compiler_abi=CompilerABI(;cxxstring_abi=:cxx03)),
                     Linux(:x86_64; compiler_abi=CompilerABI(;cxxstring_abi=:cxx11)))
        # Look across multiple gcc versions; there can be tricksy interactions here
        for gcc_version in (v"4", v"6", v"9")
            # Do each build within a separate temporary directory
            mktempdir() do build_path
                script = raw"""
                    cd ${WORKSPACE}/srcdir/cxxstringabi_tests
                    make install
                    install_license /usr/share/licenses/libuv/LICENSE
                """
                build_output_meta = do_build(build_path, script, platform, gcc_version)
                # Extract our platform's build
                @test haskey(build_output_meta, platform)
                tarball_path, tarball_hash = build_output_meta[platform][1:2]
                @test isfile(tarball_path)

                # Unpack it somewhere else
                @test verify(tarball_path, tarball_hash)
                testdir = joinpath(build_path, "testdir")
                mkdir(testdir)
                unpack(tarball_path, testdir)
                prefix = Prefix(testdir)

                # Ensure that the library detects as the correct cxxstring_abi:
                readmeta(locate(libcxxstringabi_test, prefix)) do oh
                    detected_cxxstring_abi = BinaryBuilder.detect_cxxstring_abi(oh, platform)
                    @test detected_cxxstring_abi == cxxstring_abi(platform)
                end

                # Explicitly test cxx string abi mismatches
                if gcc_version > v"4"
                    script = """
                        mkdir -p \${libdir}
                        /opt/\${target}/bin/\${target}-g++ -fPIC \\
                            -D_GLIBCXX_USE_CXX11_ABI=$(cxxstring_abi(platform) == :cxx03 ? "1" : "0") \\
                            -o \${libdir}/libcxxstringabi_test.\${dlext} \\
                            -shared \${WORKSPACE}/srcdir/cxxstringabi_tests/lib.cc
                        install_license /usr/share/licenses/libuv/LICENSE
                    """
                    @test_logs (:warn, r"ignoring our choice of compiler") match_mode=:any begin
                        do_build(build_path, script, platform, gcc_version)
                    end
                end
            end
        end
    end

    # Explicitly test not setting a cxx string abi at all
    script = raw"""
        cd ${WORKSPACE}/srcdir/cxxstringabi_tests
        make install
        install_license /usr/share/licenses/libuv/LICENSE
    """
    platform = Linux(:x86_64)
    mktempdir() do build_path
        @test_logs (:warn, r"contains std::string values") match_mode=:any begin
            do_build(build_path, script, platform, v"6")
        end
    end
end

@testset "Auditor - .dll moving" begin
    for platform in [Windows(:x86_64)]
        mktempdir() do build_path
            build_output_meta = nothing
            @test_logs (:warn, r"lib/libfoo.dll should be in `bin`") (:warn, r"Simple buildsystem detected") match_mode=:any begin
                build_output_meta = autobuild(
                    build_path,
                    "dll_moving",
                    v"1.0.0",
                    GitSource[],
                    # Install a .dll into lib
                    raw"""
                    mkdir -p ${prefix}/lib
                    cc -o ${prefix}/lib/libfoo.${dlext} -shared /usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
                    install_license /usr/include/ltdl.h
                    """,
                    # Build for our platform
                    [platform],
                    # Ensure our executable products are built
                    Product[LibraryProduct("libfoo", :libfoo)],
                    # No dependencies
                    Dependency[]
                )
            end

            @test haskey(build_output_meta, platform)
            tarball_path, tarball_hash = build_output_meta[platform][1:2]
            @test isfile(tarball_path)

            # Test that `libfoo.dll` gets moved to `bin` if it's a windows
            contents = list_tarball_files(tarball_path)
            @test "bin/libfoo.$(dlext(platform))" in contents
        end
    end
end

@testset "Auditor - .dylib identity mismatch" begin
    mktempdir() do build_path
        no_id = LibraryProduct("no_id", :no_id)
        abs_id = LibraryProduct("abs_id", :wrong_id)
        wrong_id = LibraryProduct("wrong_id", :wrong_id)
        right_id = LibraryProduct("right_id", :wrong_id)
        platform = MacOS()

        build_output_meta = autobuild(
            build_path,
            "dll_moving",
            v"1.0.0",
            FileSource[],
            # Intsall a .dll into lib
            raw"""
            mkdir -p ${prefix}/lib
            SRC=/usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
            cc -o ${libdir}/no_id.${dlext} -shared $SRC
            cc -o ${libdir}/abs_id.${dlext} -Wl,-install_name,${libdir}/abs_id.${dlext} -shared $SRC
            cc -o ${libdir}/wrong_id.${dlext} -Wl,-install_name,@rpath/totally_different.${dlext} -shared $SRC
            cc -o ${libdir}/right_id.${dlext} -Wl,-install_name,@rpath/right_id.${dlext} -shared $SRC
            install_license /usr/include/ltdl.h
            """,
            # Build for MacOS
            [platform],
            # Ensure our executable products are built
            Product[no_id, abs_id, wrong_id, right_id],
            # No dependencies
            Dependency[],
        )

        # Extract our platform's build
        @test haskey(build_output_meta, platform)
        tarball_path, tarball_hash = build_output_meta[platform][1:2]
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        prefix = Prefix(testdir)

        # Helper to extract the dylib id of a path
        function get_dylib_id(path)
            return readmeta(path) do oh
                dylib_id_lcs = [lc for lc in MachOLoadCmds(oh) if isa(lc, MachOIdDylibCmd)]
                @test !isempty(dylib_id_lcs)
                return dylib_name(first(dylib_id_lcs))
            end
        end

        # Locate the build products within the prefix, ensure that all the dylib ID's
        # now match the pattern `@rpath/$(basename(p))`
        no_id_path = locate(no_id, prefix; platform=platform)
        abs_id_path = locate(abs_id, prefix; platform=platform)
        right_id_path = locate(right_id, prefix; platform=platform)
        for p in (no_id_path, abs_id_path, right_id_path)
            @test any(startswith.(p, libdirs(prefix)))
            @test get_dylib_id(p) == "@rpath/$(basename(p))"
        end

        # Only if it already has an `@rpath/`-ified ID, it doesn't get touched.
        wrong_id_path = locate(wrong_id, prefix; platform=platform)
        @test any(startswith.(wrong_id_path, libdirs(prefix)))
        @test get_dylib_id(wrong_id_path) == "@rpath/totally_different.dylib"
    end
end

@testset "Auditor - absolute paths" begin
    mktempdir() do build_path
        sharedir = joinpath(realpath(build_path), "share")
        mkpath(sharedir)
        open(joinpath(sharedir, "foo.conf"), "w") do f
            write(f, "share_dir = \"$sharedir\"")
        end

        # Test that `audit()` warns about an absolute path within the prefix
        @test_logs (:warn, r"share/foo.conf contains an absolute path") match_mode=:any begin
            BinaryBuilder.audit(Prefix(build_path); verbose=true)
        end
    end
end

@testset "Auditor - broken symlinks" begin
    mktempdir() do build_path
        bindir = joinpath(realpath(build_path), "bin")
        mkpath(bindir)
        # Test both broken and working (but external) symlinks
        symlink("../../artifacts/1a2b3/lib/libzmq.dll.a", joinpath(bindir, "libzmq.dll.a"))
        # The following symlinks shouldn't raise a warning
        symlink("/bin/bash", joinpath(bindir, "bash.exe"))
        symlink("libfoo.so.1.2.3", joinpath(bindir, "libfoo.so"))

        # Test that `audit()` warns about broken symlinks
        @test_logs (:warn, r"Broken symlink: bin/libzmq.dll.a") match_mode=:any begin
            BinaryBuilder.warn_deadlinks(build_path)
        end
    end
end

@testset "Auditor - gcc version" begin
    # These tests assume our gcc version is concrete (e.g. that Julia is linked against libgfortran)
    our_libgfortran_version = libgfortran_version(compiler_abi(platform))
    @test our_libgfortran_version != nothing

    mktempdir() do build_path
        hello_world = ExecutableProduct("hello_world_fortran", :hello_world_fortran)
        build_output_meta = @test_logs (:warn, r"CompilerSupportLibraries_jll") match_mode=:any begin
            autobuild(
                build_path,
                "hello_fortran",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build the test suite, install the binaries into our prefix's `bin`
                raw"""
                    # Build fortran hello world
                    make -j${nproc} -sC /usr/share/testsuite/fortran/hello_world install
                    # Install fake license just to silence the warning
                    install_license /usr/share/licenses/libuv/LICENSE
                    """,
                # Build for our platform
                [platform],
                #
                Product[hello_world],
                # No dependencies
                Dependency[];
            )
        end

        # Extract our platform's build, run the hello_world tests:
        output_meta = select_platform(build_output_meta, platform)
        @test output_meta != nothing
        tarball_path, tarball_hash = output_meta[1:2]

        # Ensure the build products were created
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)

        # Attempt to run the executable, we expect it to work since it's our platform:
        hello_world_path = locate(hello_world, Prefix(testdir); platform=platform)
        with_libgfortran() do
            @test strip(String(read(`$hello_world_path`))) == "Hello, World!"
        end

        # If we audit the testdir, pretending that we're trying to build an ABI-agnostic
        # tarball, make sure it warns us about it.
        @test_logs (:warn, r"links to libgfortran!") match_mode=:any begin
            BinaryBuilder.audit(Prefix(testdir); platform=BinaryBuilder.abi_agnostic(platform), autofix=false)
        end
    end
end

@testset "Auditor - soname matching" begin
    mktempdir() do build_path
        build_output_meta = nothing
        linux_platform = Linux(:x86_64)
        @test_logs (:info, r"creating link to libfoo\.so\.1\.0\.0") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "soname_matching",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build the library only with the versioned name
                raw"""
                mkdir -p ${prefix}/lib
                cc -o ${prefix}/lib/libfoo.${dlext}.1.0.0 -fPIC -shared /usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
                # Set the soname to a non-existing file
                patchelf --set-soname libfoo.so ${prefix}/lib/libfoo.${dlext}.1.0.0
                """,
                # Build for Linux
                [linux_platform],
                # Ensure our executable products are built
                [LibraryProduct("libfoo", :libfoo)],
                # No dependencies
                Dependency[];
                autofix = true,
                verbose = true,
                require_license = false
            )
        end
        # Extract our platform's build
        @test haskey(build_output_meta, linux_platform)
        tarball_path, tarball_hash = build_output_meta[linux_platform][1:2]
        # Ensure the build products were created
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        @test readlink(joinpath(testdir, "lib", "libfoo.so")) == "libfoo.so.1.0.0"
    end
end


text:
# Methods that need to be implemented:
#   •  starting_address(::AbstractHamiltonian) - not needed
"""
    G2Correlator(d::Int,c=:cross) <: AbstractHamiltonian{ComplexF64}

Two-body correlation operator representing the density-density
correlation at distance `d` of a two component system
in a momentum-space Fock-state basis.
It returns a `Complex` value.

Correlation across two components:
```math
\\hat{G}^{(2)}(d) = \\frac{1}{M}\\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} b^†_{p}  b_q a_r δ_{s+p,q+r}
```
Correlation within a single component:
```math
\\hat{G}^{(2)}(d) = \\frac{1}{M}\\sum_{spqr=1}^M e^{-id(p-q)2π/M} a^†_{s} a^†_{p}  a_q a_r δ_{s+p,q+r}
```
# Arguments
- `d::Integer`: the distance between two particles.
- `c`: possible instructions: `:cross`: default instruction, computing correlation between particles across two components;
  `:first`: computing correlation between particles within the first component;
  `:second`: computing correlation between particles within the second component.
  These are the only defined instructions, using anything else will produce errors.

# To use on a one-component system

For a system with only one component, e.g. with `BoseFS`, the second argument `c` is irrelevant
and can be any of the above instructions, one could simply skip this argument and let it be the default value.

# See also

* [`BoseHubbardMom1D2C`](@ref)
* [`BoseFS2C`](@ref)
* [`AbstractHamiltonian`](@ref)

"""
struct G2Correlator{C} <: AbstractHamiltonian{ComplexF64}
    d::Int

    function G2Correlator(d,c=:cross)
        if c == :first
            return new{1}(d)
        elseif c == :second
            return new{2}(d)
        elseif c == :cross
            return new{3}(d)
        else
            error("Unknown instruction for G2Correlator!")
        end
    end
end

function Base.show(io::IO, g::G2Correlator{C}) where C
    if C == 1
        print(io, "G2Correlator($(g.d),:first)")
    elseif C == 2
        print(io, "G2Correlator($(g.d),:second)")
    elseif C == 3
        print(io, "G2Correlator($(g.d),:cross)")
    end
end

num_offdiagonals(g::G2Correlator{1},add::BoseFS2C) = num_offdiagonals(g, add.bsa)
num_offdiagonals(g::G2Correlator{2},add::BoseFS2C) = num_offdiagonals(g, add.bsb)

function num_offdiagonals(g::G2Correlator{3}, add::BoseFS2C)
    m = num_modes(add)
    sa = num_occupied_modes(add.bsa)
    sb = num_occupied_modes(add.bsb)
    return sa*(m-1)*sb
    # number of excitations that can be made
end

function num_offdiagonals(g::G2Correlator, add::BoseFS)
    m = num_modes(add)
    singlies, doublies = num_singly_doubly_occupied_sites(add)
    return singlies*(singlies-1)*(m - 2) + doublies*(m - 1)
end

"""
    diagonal_element(g::G2Correlator, add::BoseFS2C{NA,NB,M,AA,AB})

The diagonal element in [`G2Correlator`](@ref), where `(p-q)=0`, hence
it becomes
```math
\\frac{1}{M}\\sum_{k,p=1}^M a^†_{k} b^†_{p}  b_p a_k .
```
# See also

* [`G2Correlator`](@ref)
"""
diagonal_element(g::G2Correlator{1}, add::BoseFS2C) = diagonal_element(g, add.bsa)
diagonal_element(g::G2Correlator{2}, add::BoseFS2C) = diagonal_element(g, add.bsb)

function diagonal_element(g::G2Correlator{3}, add::BoseFS2C{NA,NB,M,AA,AB}) where {NA,NB,M,AA,AB}
    onrep_a = onr(add.bsa)
    onrep_b = onr(add.bsb)
    gd = 0
    for p in 1:M
        iszero(onrep_b[p]) && continue
        for k in 1:M
            gd += onrep_a[k]*onrep_b[p] # b†_p b_p a†_k a_k
        end
    end
    return ComplexF64(gd/M)
end

function diagonal_element(g::G2Correlator, add::BoseFS{N,M,A}) where {N,M,A}
    onrep = onr(add)
    gd = 0
    for p in 1:M
        iszero(onrep[p]) && continue
        for k in 1:M
            gd += onrep[k]*onrep[p] # a†_p a_p a†_k a_k
        end
    end
    return ComplexF64(gd/M)
end

function get_offdiagonal(g::G2Correlator{1}, add::A, chosen)::Tuple{A,ComplexF64} where A<:BoseFS2C
    new_bsa, elem = get_offdiagonal(g, add.bsa, chosen)
    return A(new_bsa,add.bsb), elem
end

function get_offdiagonal(g::G2Correlator{2}, add::A, chosen)::Tuple{A,ComplexF64} where A<:BoseFS2C
    new_bsb, elem = get_offdiagonal(g, add.bsb, chosen)
    return A(add.bsa, new_bsb), elem
end

function get_offdiagonal(
    g::G2Correlator{3},
    add::A,
    chosen,
    ma=OccupiedModeMap(add.bsa),
    mb=OccupiedModeMap(add.bsb),
)::Tuple{A,ComplexF64} where {A<:BoseFS2C}

    m = num_modes(add)
    new_bsa, new_bsb, gamma, _, _, Δp = momentum_transfer_excitation(
        add.bsa, add.bsb, chosen, ma, mb
    )
    gd = exp(-im*g.d*Δp*2π/m)*gamma
    return A(new_bsa, new_bsb), ComplexF64(gd/m)
end

function get_offdiagonal(
    g::G2Correlator,
    add::A,
    chosen,
)::Tuple{A,ComplexF64} where {A<:BoseFS}
    M = num_modes(add)
    new_add, gamma, Δp = momentum_transfer_excitation(add, chosen, OccupiedModeMap(add))
    gd = exp(-im*g.d*Δp*2π/M)*gamma
    return new_add, ComplexF64(gd/M)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 773a7b60-5a6c-11eb-0a88-1de623c4a412
begin
	using Revise
	using Pkg
	Pkg.activate("..")
	Pkg.instantiate()
	using ACCESS
	using CSV
	using StatsPlots
	using Measurements
	using CCDReduction
	using PlutoUI
	using BenchmarkTools
	using Glob, FITSIO, DataFramesMeta
	using Statistics
end

# ╔═╡ 7a8509df-2c66-456c-94ff-42b52fabc91b
using PhysicalConstants

# ╔═╡ edbe9545-e948-4cdb-94c9-eea816b8cb7e
md"""
# 🪐 ACCESS Notebook
"""

# ╔═╡ 24519e4d-5cfc-4ff7-acec-632088d84757
TableOfContents()

# ╔═╡ a457660a-7167-11eb-0eb6-adcef9cf0bb9
md"""## $(@bind run_data_info CheckBox()) Data info

View header information
"""

# ╔═╡ c452dfb6-7167-11eb-3b93-dba356a52ec3
if run_data_info
	const DATA_DIR = "/home/mango/data/WASP107/ut200302"
	df = fitscollection(
		DATA_DIR;
		recursive=false,
		exclude=r"^((?!ift[0-9]{4}c1.fits).)*$",
		exclude_key=("", "COMMENT"),
	)

	# Select header items to show
	cols = [
		"FILENAME",
		"OBJECT",
		"DISPERSR",
		"SLITMASK",
		"SPEED",
		"FILTER",
		"EXPTIME",
		"BINNING",
		"UT-DATE",
		"UT-TIME",
		"RA",
		"DEC",
		"EQUINOX",
		"EPOCH",
		"AIRMASS",
		"OBSERVER",
	]
	night_log = select(df, cols)

	#CSV.write("/home/mango/Desktop/test.csv", night_log; delim="    ")
end

# ╔═╡ a4eadc9c-cb5d-4c76-a313-3e207182f70c
if run_data_info
	with_terminal() do
		for df in groupby(night_log, "OBJECT")
			println(
				unique(df[!, "OBJECT"]),
				" ",
				length(df[!, "OBJECT"]),
				" ",
				df[begin, "FILENAME"],
				"-",
				df[end, "FILENAME"],
				" ",
				df[begin, "UT-DATE"],
				" ",
				df[begin, "UT-TIME"],
				" ",
				df[begin, "BINNING"],
			)
			for col in ["SLITMASK", "DISPERSR", "FILTER", "EXPTIME", "BINNING"]
				println(col, unique(df[!, col]))
			end
			println()
		end
	end
end

# ╔═╡ d4214693-34a5-44c3-9ba5-898b703e9503
# begin
# 	# 724 ± 7, HAT-P-26
# 	# 136 ± 7, c01
# 	# 880 ± 7, c03
# 	# 510 ± 7, c05
# 	spec_2ds = []
# 	for df in eachrow(df_sci)
# 		get_2D_spectra!(spec_2ds, df; x_cen=510)
# 	end
# end

# ╔═╡ 56819d41-5e06-463b-b927-969eac0c6da9
# begin
# 	top_5s = []

# 	for spec_2d in spec_2ds
# 		top_5 = partialsort(vec(spec_2d), 1:5; rev=true)
# 		push!(top_5s, top_5)
# 	end
# end

# ╔═╡ 0b96bbab-fd2d-4cea-967b-b19293c9ac82
# let
# 	title = "c05"
# 	p = scatter(
# 		ylims = (0, 65_000),
# 		xlabel = "Science frame",
# 		ylabel = "Top 5 counts",
# 		legend = false,
# 		title = title,
# 	)

# 	for (i, top_5) in enumerate(top_5s)
# 		scatter!(p, fill(i, 5), top_5s[i], c=pal[4])
# 	end

# 	savefig(p, "/home/mango/Desktop/$(title)_top5.png")

# 	p
# end

# ╔═╡ 9cb16d97-6116-4c94-a1ec-356f1a97ee15
# function get_2D_spectra!(spec_2ds, df; x_cen=0, app=7)
# 	FITS(df.path) do f
# 		data = read(f[1])[x_cen-app÷2 : x_cen+app÷2, :]
# 		push!(spec_2ds, data)
# 	end
# end

# ╔═╡ 5031a1ce-7803-11eb-15b8-4f3d867409b4
md"""
### $(@bind run_init_spec CheckBox()) Intitial Extracted spectra
From spec.fits files
"""

# ╔═╡ b93d53f6-7809-11eb-1924-47c439b506a9
@bind extr_spec_key Select([
	"Wavelength",
	"Simple extracted object spectrum",
	"Simple extracted flat spectrum",
	"Pixel sensitivity (obtained by the flat)",
	"Simple extracted object spectrum/pixel sensitivity",
	"Sky flag (0 = note uneven sky, 1 = probably uneven profile,
		2 = Sky_Base failed)",
	"Optimally extracted object spectrum",
	"Optimally extracted object spectrum/pixel sensitivity",
])

# ╔═╡ 12f69690-7815-11eb-2208-3563a75de8e9
md"""
### $(@bind run_trace CheckBox()) Trace
View spectra traces
"""

# ╔═╡ f9cbabb3-7a98-46fd-adee-0ea2a1e53284
md"""## $(@bind run_GPT_WLC CheckBox()) GPT WLC

Visualize GPT data stored in `GPT_DET_WLC`
"""

# ╔═╡ c34212e6-780a-11eb-035b-d3d3b4750145
if run_init_spec
	df_spec = fitscollection(
		"/home/mango/data/data_reductions/WASP50/$data_dir";
		recursive=false,
		#exclude=r"^((?!ift[0-9]{4}c4.fits).)*$",
		#exclude_key=("", "COMMENT"),
	)

	# p = plot()
	# plot!(
	# 	p,
	# 	extracted_spectra_med,
	# 	ribbon = σ,
	# 	#label = label,
	# 	legend = :topright,
	# )
end

# ╔═╡ 06b86f86-780e-11eb-1e2c-27150bbb2140
if run_init_spec
	spec_index = Dict(
		"Wavelength" => 1,
		"Simple extracted object spectrum" => 2,
		"Simple extracted flat spectrum" => 3,
		"Pixel sensitivity (obtained by the flat)" => 4,
		"Simple extracted object spectrum/pixel sensitivity" => 5,
		"Sky flag (0 = note uneven sky, 1 = probably uneven profile,
		2 = Sky_Base failed)" => 6,
		"Optimally extracted object spectrum" => 7,
		"Optimally extracted object spectrum/pixel sensitivity" => 8,
	)[extr_spec_key]
	# imgs = arrays(df_spec) |> collect
	# d = imgs[1][:, 2, :]
	# extracted_spectra_med = median(d, dims=1) |> vec
	# σ = std(d, dims=1) |> vec

	extracted_spec_med = map(arrays(df_spec)) do arr
		median(arr[:, spec_index, :], dims=1) |> vec
	end
end;

# ╔═╡ 20af7a2e-780e-11eb-112b-e9382852c19f
#plot(extracted_spec_med, label=reshape(df_spec.name, 1, :), legend=:outertopright)
if run_init_spec
	let
		plotly()
		names = map(df_spec.name) do arr
					split(arr, "_spec")[1]
				end

		plot(
			extracted_spec_med,
			label=reshape(names, 1, :),
			legend = :outertopleft,
			palette = :Paired,
			title = "$data_dir: $extr_spec_key",
		)
	end
end

# ╔═╡ cb691d70-7814-11eb-1e86-5d4d598a3ace
if run_trace
	DATA_RED = "/home/mango/data/data_reductions/WASP50/$data_dir"
	XX = load_pickle("$DATA_RED/XX.pkl")
    YY = load_pickle("$DATA_RED/YY.pkl")

	md"""
	**Number of traces to show for:**
	$(@bind trace_key Select(sort(collect([Pair(key,key) for key in keys(XX)]))))
	$(@bind num_traces Slider(2:7, show_value=true))
	"""
end

# ╔═╡ 48714914-7815-11eb-388b-2d69581df140
if run_trace
	gr()
	trace_idxs = round.(Int, range(1, length(XX[trace_key]), length=num_traces))
	plot(
		XX[trace_key][trace_idxs],
		YY[trace_key][trace_idxs],
		label = reshape(trace_idxs, 1, :) .- 1,
		title = trace_key,
		legend = :topleft,
		legend_title = "index",
		palette = palette(:diverging_isoluminant_cjo_70_c25_n256, num_traces),
	)
end

# ╔═╡ 2a9ead8e-619d-11eb-15ee-43f8d97c4f19
md"""## $(@bind run_retrievals CheckBox()) `exoretrievals`

Summarize Bayesian evidences returned by `exoretrievals`
"""

# ╔═╡ 4ffdfbc0-619d-11eb-1506-65e1a4efabb7
if run_retrievals
	const BASE_DIR = "data_retrievals/spot_lower_bound"
	const model_types = (
		clear = "HATP23_E1_NoHet_FitP0_NoClouds_NoHaze_fitR0",
		haze = "HATP23_E1_NoHet_FitP0_NoClouds_Haze_fitR0",
		spot = "HATP23_E1_Het_FitP0_NoClouds_NoHaze_fitR0",
		spot_haze = "HATP23_E1_Het_FitP0_NoClouds_Haze_fitR0",
	)

	retrieval_data = load_retrieval_data(BASE_DIR, model_types)

	Z = get_evidences(retrieval_data)

	min_Z, min_Z_loc = findmin(Z)

	ΔlnZ = Z .- min_Z
end

# ╔═╡ 84b6b120-61a9-11eb-2bc8-df73e542a46d
if run_retrievals
	species = names(Z, 1)
	models = names(Z, 2)
	ΔlnZ_mat = hcat(Symbol.(species), ΔlnZ.array)
	min_model = species[min_Z_loc[1]], models[min_Z_loc[2]]
	md"Lowest evidence: $(min_model)"
end

# ╔═╡ 4ebb0012-61a9-11eb-1513-678e52371305
if run_retrievals
	md"""
	### Copy paste ``\LaTeX`` table
	"""
end

# ╔═╡ 14d5c3b2-61a8-11eb-349f-cb24dc8b57ea
if run_retrievals
	with_terminal() do
		lap(ΔlnZ_mat)
	end
end

# ╔═╡ 78561026-61b6-11eb-2100-153023cf94cf
if run_retrievals
	md"""
	### Evidences graph
	"""
end

# ╔═╡ 63378968-61b1-11eb-2ca4-d35d300a7b68
if run_retrievals
	groupedbar(
		species,
		ΔlnZ;
		label = reshape(models, 1, :),
		lw = 0,
		xlabel = "Species",
		ylabel = "Relative change in log-evidence",
		legend = :outertopright,
		title = BASE_DIR,
		fmt = :png,
		dpi = 250,
	)
end

# ╔═╡ a7f5ccb5-3a63-4197-91cf-942d894f58ba
md"""
## Packages
"""

# ╔═╡ 0022d8c7-cafa-4964-a30b-35669ec6c285
import DarkMode

# ╔═╡ c76e0088-2a1e-4084-9d8b-88d5d4564bc7
DarkMode.enable()

# ╔═╡ 110457da-61bf-11eb-1a02-916d60cb2096
begin
	const pal = palette(["#f5d300", "#d55e00", "#08f7fe", "#0173b2", "#029e73"])
	theme(:dark)
	default(
		titlefont = "Lato",
		guidefont = "Lato",
		markerstrokewidth = 0,
		palette = pal,
		dpi = 250,
		fmt = :png,
	)
end

# ╔═╡ 4903d454-7811-11eb-1cc2-4d1fa33da6e1
@bind data_dir Select(readdir("/home/mango/data/data_reductions/WASP50"))

# ╔═╡ 1061134b-a53c-4abb-96b5-b95def13f02b
if run_GPT_WLC
	data_dir = "/home/mango/data/data_detrending/WASP50b"
	out_folder = "$data_dir/out_l/WASP50/w50_131219/white-light"

# 	# Load pickle
# 	gpts_pkl_to_dict = load_pickle(

# )
# 	# Show keys
# 	println(keys(gpts_pkl_to_dict))

# 	# Load specified params
# 	const gpts_params = ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"]
# 	gpts_params_dict = sub_dict(
# 		gpts_pkl_to_dict, gpts_params
# 	)

# 	# Add desired offsets
# 	gpts_params_dict[:t0] .-= 2.45485e6; gpts_params_dict

# 	# Load into DataFrame
# 	df_gpts = DataFrame(gpts_params_dict)
end

# ╔═╡ Cell order:
# ╟─edbe9545-e948-4cdb-94c9-eea816b8cb7e
# ╠═24519e4d-5cfc-4ff7-acec-632088d84757
# ╟─a457660a-7167-11eb-0eb6-adcef9cf0bb9
# ╠═c452dfb6-7167-11eb-3b93-dba356a52ec3
# ╠═a4eadc9c-cb5d-4c76-a313-3e207182f70c
# ╠═d4214693-34a5-44c3-9ba5-898b703e9503
# ╠═56819d41-5e06-463b-b927-969eac0c6da9
# ╠═0b96bbab-fd2d-4cea-967b-b19293c9ac82
# ╠═9cb16d97-6116-4c94-a1ec-356f1a97ee15
# ╟─5031a1ce-7803-11eb-15b8-4f3d867409b4
# ╠═4903d454-7811-11eb-1cc2-4d1fa33da6e1
# ╟─b93d53f6-7809-11eb-1924-47c439b506a9
# ╠═c34212e6-780a-11eb-035b-d3d3b4750145
# ╠═06b86f86-780e-11eb-1e2c-27150bbb2140
# ╠═20af7a2e-780e-11eb-112b-e9382852c19f
# ╟─12f69690-7815-11eb-2208-3563a75de8e9
# ╠═cb691d70-7814-11eb-1e86-5d4d598a3ace
# ╠═48714914-7815-11eb-388b-2d69581df140
# ╟─f9cbabb3-7a98-46fd-adee-0ea2a1e53284
# ╠═1061134b-a53c-4abb-96b5-b95def13f02b
# ╟─2a9ead8e-619d-11eb-15ee-43f8d97c4f19
# ╠═84b6b120-61a9-11eb-2bc8-df73e542a46d
# ╠═4ffdfbc0-619d-11eb-1506-65e1a4efabb7
# ╠═4ebb0012-61a9-11eb-1513-678e52371305
# ╠═14d5c3b2-61a8-11eb-349f-cb24dc8b57ea
# ╠═78561026-61b6-11eb-2100-153023cf94cf
# ╠═63378968-61b1-11eb-2ca4-d35d300a7b68
# ╟─a7f5ccb5-3a63-4197-91cf-942d894f58ba
# ╠═773a7b60-5a6c-11eb-0a88-1de623c4a412
# ╠═0022d8c7-cafa-4964-a30b-35669ec6c285
# ╠═c76e0088-2a1e-4084-9d8b-88d5d4564bc7
# ╠═110457da-61bf-11eb-1a02-916d60cb2096
# ╠═7a8509df-2c66-456c-94ff-42b52fabc91b


text:
using BenchmarkTools, LinearAlgebra, Plots
gr()
import JuliaBLAS: addmul!, Block

function bench(nn=30)
    time2gflops(mnk, time) = 2 * mnk^3 / time

    BLAS.set_num_threads(1)
    mnks = 48*(1:2:nn)
    obtimes = zeros(length(mnks)); jbtimes = zeros(length(mnks))

    for (i,siz) in enumerate(mnks)
        A,B,C = (zeros(siz,siz) for i in 1:3)
        bj = @benchmarkable addmul!($C,$A,$B) samples=20 time_tolerance=0.1 seconds=0.8
        jbtimes[i] = run(bj) |> minimum |> time
        bo = @benchmarkable mul!($C,$A,$B) samples=20 time_tolerance=0.1 seconds=0.8
        obtimes[i] = run(bo) |> minimum |> time
    end

    jflops = time2gflops.(mnks, jbtimes)
    oflops = time2gflops.(mnks, obtimes)
    plot(mnks, jflops, lab="JuliaBLAS")
    plot!(mnks, oflops, lab="OpenBLAS", ylabel="GFLOPS", xlabel="M=N=K", legend=:bottomright, dpi=400, ylims=(0,60), yticks=0:5:50)
end
plt = bench()
savefig(plt, "bench1.png")


text:
# ---
# title: 139. Word Break
# id: problem139
# author: zhwang
# date: 2022-02-18
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/word-break/description/>
# hidden: true
# ---
#
# Given a **non-empty** string _s_ and a dictionary _wordDict_ containing a list
# of **non-empty** words, determine if _s_ can be segmented into a space-
# separated sequence of one or more dictionary words.
#
# **Note:**
#
#   * The same word in the dictionary may be reused multiple times in the segmentation.
#   * You may assume the dictionary does not contain duplicate words.
#
# **Example 1:**
#
#
#
#     Input: s = "leetcode", wordDict = ["leet", "code"]
#     Output: true
#     Explanation: Return true because "leetcode" can be segmented as "leet code".
#
#
# **Example 2:**
#
#
#
#     Input: s = "applepenapple", wordDict = ["apple", "pen"]
#     Output: true
#     Explanation: Return true because "applepenapple" can be segmented as "apple pen apple".
#                  Note that you are allowed to reuse a dictionary word.
#
#
# **Example 3:**
#
#
#
#     Input: s = "catsandog", wordDict = ["cats", "dog", "sand", "and", "cat"]
#     Output: false
#
#
#
## @lc code=start
using LeetCode

## BFS
function word_break_bfs(s::String, word_dict::Vector{String})::Bool
    n, word_dict = length(s), sort!(word_dict; by=x -> length(x))
    valids, valid_pos = fill(false, n), [0]
    while !isempty(valid_pos)
        pos = popfirst!(valid_pos)
        for word in word_dict
            (new_pos = pos + length(word)) > n && break
            new_pos == n && s[(pos + 1):new_pos] == word && return true
            if !valids[new_pos] && s[(pos + 1):new_pos] == word
                valids[new_pos] = true
                push!(valid_pos, new_pos)
            end
        end
    end
    return false
end

## Dynamic Programming
function word_break(s::String, word_dict::Vector{String})::Bool
    lens = sort!(unique!(length.(word_dict)))
    dp = append!([true], fill(false, length(s)))
    for i in eachindex(dp)
        for len in lens
            i > len || break
            dp[i - len] && s[(i - len):(i - 1)] ∈ word_dict && (dp[i] = true)
        end
    end
    return last(dp)
end
## @lc code=end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using ArraysOfArrays, Distributions, PDMats

@testset "distribution_density" begin
    mvt = @inferred MvTDist(1.5, PDMat([2.0 1.0; 1.0 3.0]))
    mvdd = @inferred BAT.DistributionDensity(mvt)

    @testset "properties" begin
        @test typeof(mvdd) <: AbstractDensity
        @test parent(mvdd) == mvt
        @test nparams(mvdd) == 2
    end

    @testset "BAT.density_logval" begin
        @test (@inferred BAT.density_logval(mvdd, [0.0, 0.0])) ≈ -2.64259602
    end
end