text:
using AstroImages
#] pkg> add https://github.com/JuliaAstro/AstroImages.jl

using FITSIO
using LinearAlgebra
using NaNMath; nm = NaNMath
using OpenCL
#using Plots
using Makie
using Random
using Statistics
using Wavelets

#pyplot()

dir = "/home/chris/ダウンロード"
#file = "e20121211_0010500001_dp_sf_st_mos.fits"
file = "ALMA01030862.fits"
#file = "ALMA01157085.fits"
#file = "Cygnus_sp46_vs-150_ve100_dv0.50_CN_Tmb.fits.gz" #good for testing noisy datasets

#dir = "/home/chris/NAO/NRO/SF/"

fitspath = string(dir, "/", file)

#fitspath = "/home/chris/NAO/SubaruWebQL/FITSCACHE/SUPM139C52520E922500.fits"
#fitspath = "/home/chris/projects/fits_web_ql/FITSCACHE/d61977e0-fada-3cee-8198-31c00ea64ab3.fits"

println(ARGS, "\targc:", size(ARGS)[1])

if size(ARGS)[1] > 0
    fitspath = ARGS[1]
end

println("FITS file:", fitspath)

#fits = load(fitspath)
#img = AstroImage(fitspath)
#plot(img)

fits = FITS(fitspath)
println(fits)
println(fits[1])

data = @time read(fits[1])
dim = minimum(size(data))
println(size(data), " the smaller dim:", dim)
range = @time nm.extrema(data)
println("min,max:", range)
med = @time median(data)
println("median:", med)

#=
wt = wavelet(WT.cdf97, WT.Lifting)
# a full wavelet transform
#... array must be square/cube ...
sqdata = @view data[1:dim,1:dim]
#@time dwt!(sqdata, wt) #no support for ::GLS Float32 in-place transform
xt = @time dwt(sqdata, wt)
println(size(xt))

# wavelet packet transform
# execution errors (GSL only supports double!!!)
wt = wavelet(WT.db2, WT.Lifting)
L = 4
@time wpt(sqdata, wt, L)
=#

function rbf_gradient_julia(_x1, _x2, _y, _data, _e, c1, c2, p0, p1, p2, w, _grad_c1, _grad_c2, _grad_p0, _grad_p1, _grad_p2, _grad_w)
    for index = 1:length(_data)
        x1 = _x1[index]
        x2 = _x2[index]

        tmp = w[NCLUST + 1]
        grad_c1 = zeros(NCLUST)
        grad_c2 = zeros(NCLUST)
        grad_p0 = zeros(NCLUST)
        grad_p1 = zeros(NCLUST)
        grad_p2 = zeros(NCLUST)
        grad_w = zeros(NCLUST + 1)
        grad_w[NCLUST + 1] = 1.0

        for i = 1:NCLUST
            a = exp(p0[i])
            b = p1[i]
            c = exp(p2[i])

            tmp1 = (x1 - c1[i]);
            tmp2 = (x2 - c2[i]);
            dist = a * tmp1 * tmp1 - 2.0 * b * tmp1 * tmp2 + c * tmp2 * tmp2;
            act = exp(-dist);
            tmp += w[i] * act;

            #gradients
            grad_w[i] = act;
            grad_c1[i] = 2.0 * w[i] * act * (a * tmp1 - b * tmp2)
           	grad_c2[i] = 2.0 * w[i] * act * (c * tmp2 - b * tmp1)
           	grad_p0[i] = w[i] * act * a * ( - tmp1 * tmp1 )
           	grad_p1[i] = 2.0 * w[i] * act * tmp1 * tmp2
           	grad_p2[i] = w[i] * act * c * ( -tmp2 * tmp2)
        end

        e = tmp - _data[index];
        _y[index] = tmp;
        _e[index] = e;

        #gradients
        for i = 1:NCLUST + 1
            _grad_w[i] += e * grad_w[i]
        end

        for i = 1:NCLUST
            _grad_c1[i] += e * grad_c1[i]
            _grad_c2[i] += e * grad_c2[i]
            _grad_p0[i] += e * grad_p0[i]
            _grad_p1[i] += e * grad_p1[i]
            _grad_p2[i] += e * grad_p2[i]
        end

        #println("index: $(index), x1: $(x1), x2: $(x2), y: $(tmp), e: $(e)")
    end
end

dims = size(data)
width = dims[1]
height = dims[2]
depth = dims[3]
capacity = width * height

XCLUST = min(Int(round(width / 8)), 32)#/16
YCLUST = min(Int(round(height / 8)), 32)#/16
NCLUST = XCLUST * YCLUST
NITER = 500

println("width : ", width, "\theight : ", height, "\tdepth : ", depth, "\tcapacity : ", capacity)
println("XCLUST : ", XCLUST, "\tYCLUST : ", YCLUST, "\tNCLUST : ", NCLUST)

device, ctx, queue = cl.create_compute_context()
println(device)

compression_code = open("rbf.cl") do file
    #define NCLUST dynamically
    "#define NCLUST $(NCLUST)\n" * read(file, String)
end

program = cl.Program(ctx, source = compression_code) |> cl.build!
rbf_gradient = cl.Kernel(program, "rbf_gradient")
rbf_compute = cl.Kernel(program, "rbf_compute")

#scene = Scene(resolution = (500, 500))
#center!(scene)

for frame = Int(round(depth / 2)):Int(round(depth / 2))#1:depth
#for frame = 1:5
    sub = view(data, :, :, frame, 1)
    println("frame : ", frame, "\tdims: ", size(sub))
    (frame_min, frame_max) = @time nm.extrema(sub)

    x1 = Float32[]
    x2 = Float32[]
    d = Float32[]
    y = Float32[]
    e = Float32[]
    count = 0

    x1test = Float32[]
    x2test = Float32[]

    for iy = 1:height
        for ix = 1:width
            tmp = sub[ix,iy]
            if !isnan(tmp)
                count = count + 1
                pixel = log(0.5 + (tmp - frame_min) / (frame_max - frame_min))
                push!(d, pixel)
                push!(x1, (ix - 1) / (width - 1))   #[0,1]
                push!(x2, (iy - 1) / (height - 1))  #[0,1]
                push!(y, 0)
                push!(e, 0)
            end
            push!(x1test, (ix - 1) / (width - 1))   #[0,1]
            push!(x2test, (iy - 1) / (height - 1))  #[0,1]
        end
    end

    (x1min, x1max) = nm.extrema(x1)
    (x2min, x2max) = nm.extrema(x2)

    println("count : $(count), $(x1min), $(x1max), $(x2min), $(x2max)")

    #training data OK

    #parameter initialisation
    c1 = shuffle(x1)[1:NCLUST]
    c2 = shuffle(x2)[1:NCLUST]

    #=
    c1 = Float32[]
    c2 = Float32[]

    for i = 1:XCLUST
        for j = 1:YCLUST
            push!(c1, (i - 1) / (XCLUST - 1))
            push!(c2, (j - 1) / (YCLUST - 1))
        end
    end
    =#

    scene = Scene(resolution = (1500, 1500))
    center!(scene)
    scatter!(scene, c1, c2, markersize = 1 / NCLUST)
    display(scene)
    #gui()

    σX = 0.1 / (XCLUST - 1)
    σY = 0.1 / (YCLUST - 1)
    σX = randn(NCLUST) * (0.01 * σX) .+ σX
    σY = randn(NCLUST) * (0.01 * σY) .+ σY
    θ = 2 * π * rand(NCLUST)

    a = 0.5 * cos.(θ) .* cos.(θ) ./ (σX .* σX) .+ 0.5 * sin.(θ) .* sin.(θ) ./ (σY .* σY)
    b = -0.25 * sin.(2.0 * θ) ./ (σX .* σX) .+ 0.25 * sin.(2.0 * θ) ./ (σY .* σY)
    c = 0.5 * sin.(θ) .* sin.(θ) ./ (σX .* σX) .+ 0.5 * cos.(θ) .* cos.(θ) ./ (σY .* σY)

    p0 = log.(a)
    p1 = b
    p2 = log.(c)

    w = randn(Float32, NCLUST + 1)

    #learning rates etc.
    η₊ = 1.2
    η₋ = 0.5
    ΔMin = 1e-5
    ΔMax = 1e-1
    #ΔMin = 1e-7#1e-5
    #ΔMax = 1e-3#1e-1
    d₀ = ΔMin

    #previous gradients
    grad_c1_prev = zeros(Float32, NCLUST)
    grad_c2_prev = zeros(Float32, NCLUST)
    grad_p0_prev = zeros(Float32, NCLUST)
    grad_p1_prev = zeros(Float32, NCLUST)
    grad_p2_prev = zeros(Float32, NCLUST)
    grad_w_prev = zeros(Float32, NCLUST + 1)

    #changes
    Δc1 = fill(Float32(d₀), NCLUST)
    Δc2 = fill(Float32(d₀), NCLUST)
    Δp0 = fill(Float32(d₀), NCLUST)
    Δp1 = fill(Float32(d₀), NCLUST)
    Δp2 = fill(Float32(d₀), NCLUST)
    Δw = fill(Float32(d₀), NCLUST + 1)

    #OpenCL buffers
    data_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = d)
    x1_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = x1)
    x2_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = x2)
    y_buff = cl.Buffer(Float32, ctx, :w, length(y))
    e_buff = cl.Buffer(Float32, ctx, :w, length(e))

    x1test_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = x1test)
    x2test_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = x2test)
    ytest_buff = cl.Buffer(Float32, ctx, :w, length(x1test))

    #println("w (before):", w)

    for iter = 1:NITER
    #gradients
        grad_c1 = zeros(Float32, NCLUST)
        grad_c2 = zeros(Float32, NCLUST)
        grad_p0 = zeros(Float32, NCLUST)
        grad_p1 = zeros(Float32, NCLUST)
        grad_p2 = zeros(Float32, NCLUST)
        grad_w = zeros(Float32, NCLUST + 1)

    #gradient buffers
    #=
        grad_c1_buff = cl.Buffer(Float32, ctx, :w, length(grad_c1))
        grad_c2_buff = cl.Buffer(Float32, ctx, :w, length(grad_c2))
        grad_p0_buff = cl.Buffer(Float32, ctx, :w, length(grad_p0))
        grad_p1_buff = cl.Buffer(Float32, ctx, :w, length(grad_p1))
        grad_p2_buff = cl.Buffer(Float32, ctx, :w, length(grad_p2))
        grad_w_buff = cl.Buffer(Float32, ctx, :w, length(grad_w))
    =#
        grad_c1_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_c1)
        grad_c2_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_c2)
        grad_p0_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_p0)
        grad_p1_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_p1)
        grad_p2_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_p2)
        grad_w_buff = cl.Buffer(Float32, ctx, (:w, :copy), hostbuf = grad_w)

    #parameter buffers
        c1_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = c1)
        c2_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = c2)
        p0_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = Float32.(p0))
        p1_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = Float32.(p1))
        p2_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = Float32.(p2))
        w_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf = w)

    #execute a forward pass
        @time queue(rbf_gradient, size(d), nothing, x1_buff, x2_buff, y_buff, data_buff, e_buff, c1_buff, c2_buff, p0_buff, p1_buff, p2_buff, w_buff, grad_c1_buff, grad_c2_buff, grad_p0_buff, grad_p1_buff, grad_p2_buff, grad_w_buff)

        y = cl.read(queue, y_buff)
        e = cl.read(queue, e_buff)
        grad_c1 = cl.read(queue, grad_c1_buff)
        grad_c2 = cl.read(queue, grad_c2_buff)
        grad_p0 = cl.read(queue, grad_p0_buff)
        grad_p1 = cl.read(queue, grad_p1_buff)
        grad_p2 = cl.read(queue, grad_p2_buff)
        grad_w = cl.read(queue, grad_w_buff)

        println("frame $(frame) ==> GPU batch training iteration: $(iter), error: ", norm(e))

        #@time rbf_gradient_julia(x1, x2, y, d, e, c1, c2, p0, p1, p2, w, grad_c1, grad_c2, grad_p0, grad_p1, grad_p2, grad_w)
        #println("frame $(frame) ==> CPU batch training iteration: $(iter), error: ", norm(e))

    #update parameters
        #w
        for i = 1:NCLUST + 1
            if grad_w_prev[i] * grad_w[i] > 0.0
                Δw[i] = min(ΔMax, η₊ * Δw[i])
            end

            if grad_w_prev[i] * grad_w[i] < 0.0
                Δw[i] = max(ΔMin, η₋ * Δw[i])
                grad_w[i] = 0.0
            end

            w[i] = w[i] - sign(grad_w[i]) * Δw[i]
            grad_w_prev[i] = grad_w[i] ;
        end

        #c1, c2, p0, p1, p2
        for i = 1:NCLUST
            #c1
            if grad_c1_prev[i] * grad_c1[i] > 0.0
                Δc1[i] = min(ΔMax, η₊ * Δc1[i])
            end

            if grad_c1_prev[i] * grad_c1[i] < 0.0
                Δc1[i] = max(ΔMin, η₋ * Δc1[i])
                grad_c1[i] = 0.0
            end

            c1[i] = c1[i] - sign(grad_c1[i]) * Δc1[i]
            grad_c1_prev[i] = grad_c1[i] ;

            #c2
            if grad_c2_prev[i] * grad_c2[i] > 0.0
                Δc2[i] = min(ΔMax, η₊ * Δc2[i])
            end

            if grad_c2_prev[i] * grad_c2[i] < 0.0
                Δc2[i] = max(ΔMin, η₋ * Δc2[i])
                grad_c2[i] = 0.0
            end

            c2[i] = c2[i] - sign(grad_c2[i]) * Δc2[i]
            grad_c2_prev[i] = grad_c2[i] ;

            #p0
            if grad_p0_prev[i] * grad_p0[i] > 0.0
                Δp0[i] = min(ΔMax, η₊ * Δp0[i])
            end

            if grad_p0_prev[i] * grad_p0[i] < 0.0
                Δp0[i] = max(ΔMin, η₋ * Δp0[i])
                grad_p0[i] = 0.0
            end

            p0[i] = p0[i] - sign(grad_p0[i]) * Δp0[i]
            grad_p0_prev[i] = grad_p0[i] ;

            #p1
            if grad_p1_prev[i] * grad_p1[i] > 0.0
                Δp1[i] = min(ΔMax, η₊ * Δp1[i])
            end

            if grad_p1_prev[i] * grad_p1[i] < 0.0
                Δp1[i] = max(ΔMin, η₋ * Δp1[i])
                grad_p1[i] = 0.0
            end

            p1[i] = p1[i] - sign(grad_p1[i]) * Δp1[i]
            grad_p1_prev[i] = grad_p1[i] ;

            #p2
            if grad_p2_prev[i] * grad_p2[i] > 0.0
                Δp2[i] = min(ΔMax, η₊ * Δp2[i])
            end

            if grad_p2_prev[i] * grad_p2[i] < 0.0
                Δp2[i] = max(ΔMin, η₋ * Δp2[i])
                grad_p2[i] = 0.0
            end

            p2[i] = p2[i] - sign(grad_p2[i]) * Δp2[i]
            grad_p2_prev[i] = grad_p2[i] ;
        end

        #plot(x1, x2, ocl_y)
        #density(x1, x2, ocl_y)
        #scatter(c1, c2)
        #gui()

        #scene = Scene(resolution = (500, 500))
        #center!(scene)
        #scatter!(scene, c1, c2, markersize = 1 / NCLUST)
        #AbstractPlotting.force_update!()
        scene = Scene(resolution = (1500, 1500))
        center!(scene)

        if count < capacity
            #scatter!(scene, c1, c2, markersize = 1 / NCLUST)

            @time queue(rbf_compute, capacity, nothing, x1test_buff, x2test_buff, ytest_buff, c1_buff, c2_buff, p0_buff, p1_buff, p2_buff, w_buff)
            y = cl.read(queue, ytest_buff)
            img = reshape(y, width, height)
            heatmap!(scene, img)
        else
            img = reshape(y, width, height)
            heatmap!(scene, img)
        end

        display(scene)

        #scatter!(scene, x1, x2, ocl_y, markersize = 1 / NCLUST)

    #validate in Julia
    #=
        @time rbf_gradient_julia(x1, x2, y, d, e, c1, c2, p0, p1, p2, w, grad_c1, grad_c2, grad_p0, grad_p1, grad_p2, grad_w)

        if isapprox(norm(y - ocl_y) / length(d), zero(Float32), atol = 1e-8)
            println("y: Success!")
        else
            println("y: Norm should be ≈ 0.0f, not ", norm(y - ocl_y) / length(d))
        end

        if isapprox(norm(grad_c1 - ocl_grad_c1) / length(grad_c1), zero(Float32), atol = 1e-8)
            println("grad_c1: Success!")
        else
            println("grad_c1: Norm should be ≈ 0.0f, not ", norm(grad_c1 - ocl_grad_c1) / length(grad_c1))
        end

        if isapprox(norm(grad_c2 - ocl_grad_c2) / length(grad_c2), zero(Float32), atol = 1e-8)
            println("grad_c2: Success!")
        else
            println("grad_c2: Norm should be ≈ 0.0f, not ", norm(grad_c2 - ocl_grad_c2) / length(grad_c2))
        end

        if isapprox(norm(grad_p0 - ocl_grad_p0) / length(grad_p0), zero(Float32), atol = 1e-8)
            println("grad_p0: Success!")
        else
            println("grad_p0: Norm should be ≈ 0.0f, not ", norm(grad_p0 - ocl_grad_p0) / length(grad_p0))
        end

        if isapprox(norm(grad_p1 - ocl_grad_p1) / length(grad_p1), zero(Float32), atol = 1e-8)
            println("grad_p1: Success!")
        else
            println("grad_p1: Norm should be ≈ 0.0f, not ", norm(grad_p1 - ocl_grad_p1) / length(grad_p1))
        end

        if isapprox(norm(grad_p2 - ocl_grad_p2) / length(grad_p2), zero(Float32), atol = 1e-8)
            println("grad_p2: Success!")
        else
            println("grad_p2: Norm should be ≈ 0.0f, not ", norm(grad_p2 - ocl_grad_p2) / length(grad_p2))
        end

        if isapprox(norm(grad_w - ocl_grad_w) / length(grad_w), zero(Float32), atol = 1e-8)
            println("grad_w: Success!")
        else
            println("grad_w: Norm should be ≈ 0.0f, not ", norm(grad_w - ocl_grad_w) / length(grad_w))
        end
    =#
    end

    #println("w (after):", w)
    #println("θ (after):", θ)
end


text:
############################################################################################
# Types and structs - Resampling Methods

############################################################################################
"""
$(TYPEDEF)

Container that determines if jittering has to be applied

# Fields
$(TYPEDFIELDS)
"""
struct JitterTune{B<:UpdateBool}
    "Bolean if adaptive jittering is applied. If false, always max steps are jittered."
    adaption::B
    "Number of jitttering steps performed while jittering"
    Nsteps::Iterator
    "Jittering threshold for correlation of parameter."
    threshold::Float64
    "Minimum amount of jittering steps."
    min::Int64
    "Maximum amount of jittering steps."
    max::Int64
    function JitterTune(
        adaption::B,
        threshold::Float64,
        min::Integer=1,
        max::Integer=30
    ) where {B<:UpdateBool}
        ArgCheck.@argcheck 0 < min <= max
        ArgCheck.@argcheck 0.0 < threshold <= 1.0 "threshold needs to be positive"
        return new{B}(adaption, Iterator(0), threshold, min, max)
    end
end

############################################################################################
function update!(tune::JitterTune)
    init!(tune.Nsteps, 0)
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Check if we can stop jittering step based on correlation ρ.

# Examples
```julia
```

"""
function jitter!(tune::JitterTune, adaption::UpdateTrue, ρ::T) where {T<:AbstractFloat}
    ## Add step to Nsteps
    update!(tune.Nsteps)
    ## Check if jittering has to be continued
    #!NOTE: Need to start with >= tune.max and < tune.min conditions, and then with ρ
    if tune.Nsteps.current >= tune.max
        return false
    elseif tune.Nsteps.current < tune.min
        return true
    #!NOTE: if not above or below boundary, check if correlation threshold is fullfiled
    # If correlation smaller than threshold, stop jittering
    elseif ρ < tune.threshold
        return false
    else
        return true
    end
end

function jitter!(tune::JitterTune, adaption::UpdateFalse, ρ::T) where {T<:AbstractFloat}
    ## Add step to Nsteps
    update!(tune.Nsteps)
    ## Continue jittering until max steps are reached
    if tune.Nsteps.current >= tune.max
        return false
    else
        return true
    end
end

function jitter!(tune::JitterTune, ρ::T) where {T<:AbstractFloat}
    return jitter!(tune, tune.adaption, ρ)
end

############################################################################################
#export
export JitterTune, jitter!, update!

ext: jl


text:
#
# This is a script to generate some results for the AABI-2020 submission. It will produce
# a couple of graphs that will both provide an explanation of the experimental set up, and
# show an approximate posterior mean (because why not :shrug:).
#
# This uses my reduced scripting format that only includes a single experiment block /
# results analysis / plotting block. There's no need for the usual first block because there
# are only two things generated. It's a bit of a mush, but I was pushed for time.
#

println("pseudo_point_synthetic_posterior_plots")

using Revise

using Pkg
Pkg.activate(".");
Pkg.resolve();
Pkg.instantiate();

using ColorTypes
using DataFrames
using BenchmarkTools
using DrWatson
using FileIO
using ImageMagick
using Images
using LaTeXStrings
using LinearAlgebra
using PGFPlotsX
using Plots
using Random
using Stheno
using TemporalGPs
using TexTables

# Configure Plots.
pgfplotsx();
Plots.default(
    tickfontsize=nothing,
    legendfontsize=nothing,
);

# Use default document font-size if nothing provided. This is a hack.
function Plots.pgfx_font(fontsize::Nothing, thickness_scaling = 1, font = "\\selectfont")
    return string("{", font, "}")
end

BLAS.set_num_threads(1)

const exp_dir_name = "pseudo_point_synthetic_posterior_plots"

# A kernel used throughout the experiments.
function separable_kernel(σ::Real, λ_space::Real, λ_time::Real)
    k_space = σ^2 * stretch(EQ(), λ_space)
    k_time = stretch(Matern52(), λ_time)
    return TemporalGPs.Separable(k_space, k_time)
end

# Ensure that output directory for results exists.
const data_dir = joinpath(datadir(), exp_dir_name)
mkpath(data_dir)
mkpath(joinpath(data_dir, "meta"))

# Ensure that results directory exists.
const results_dir = joinpath(projectdir(), "results", exp_dir_name)
mkpath(results_dir)

# Locations to save the plots to.
save_dirs = [
    results_dir,
    "/Users/willtebbutt/Dropbox/University/PhD/first_year_report/chapter_3/pseudo_points"
]
foreach(mkpath, save_dirs)

# Run experiments.
let

    let
        # Set seed for reproducibility.
        rng = MersenneTwister(123456)

        # Extract experiment parameters.
        θ = (σ = 0.92, λ_space = 1.1, λ_time = 0.83, S_min = 0.1)
        T = 100
        N_space_per_time = 50
        N_missing_per_time = 5
        M_per_time = 10

        # Specify input locations. We place these irregularly through time.
        t = sort(rand(rng, T)) .* 50
        t = RegularSpacing(0.0, 0.2, T)
        r = collect(range(0.0, 5.0; length=N_space_per_time))
        r_pred = collect(range(minimum(r), maximum(r); length=100))
        x = TemporalGPs.RectilinearGrid(r, t);
        x_naive_full = collect(x)

        # Specify inducing inputs. They are regularly-spaced between 0 and 10.
        z_r = collect(range(0.0, maximum(r); length=M_per_time))
        z_naive = collect(TemporalGPs.RectilinearGrid(z_r, t))

        # Specify observation noise.
        S_full = [Diagonal(rand(rng, N_space_per_time) .+ θ.S_min) for _ in 1:T]
        S_full_diag = Diagonal(vcat(diag.(S_full)...))

        # Generate entire data set exactly / naively.
        k = separable_kernel(θ.σ, θ.λ_space, θ.λ_time)
        f = GP(k, GPC())
        fx = to_sde(f)(x, S_full_diag)

        # Generate entire data set using state-space model.
        y = rand(rng, TemporalGPs.build_lgssm(fx))

        # Remove a couple of data points at each time point directly from
        # x_naive and y_naive.
        missing_idxs = map(t -> randperm(rng, N_space_per_time)[1:N_missing_per_time], 1:T)
        x_naive_vecs = map(1:T) do t
            x_naive_full[(t - 1) * N_space_per_time + 1:t * N_space_per_time]
        end
        x_naive_missing_vecs = map(eachindex(missing_idxs)) do n
            present_idx = deleteat!(collect(1:N_space_per_time), sort(missing_idxs[n]))
            return x_naive_vecs[n][present_idx]
        end
        y_naive_missing_vecs = map(eachindex(missing_idxs)) do n
            present_idx = deleteat!(collect(1:N_space_per_time), sort(missing_idxs[n]))
            return y[n][present_idx]
        end
        S_naive_missing_vecs = map(eachindex(missing_idxs)) do n
            present_idx = deleteat!(collect(1:N_space_per_time), sort(missing_idxs[n]))
            return diag(S_full[n])[present_idx]
        end

        x_naive = vcat(x_naive_missing_vecs...)
        y_naive = vcat(y_naive_missing_vecs...)
        S_naive = vcat(S_naive_missing_vecs...)

        # Modify S where there is missing data.
        S = map(eachindex(missing_idxs)) do n
            s = copy(diag(S_full[n]))
            s[missing_idxs[n]] .= 1e9
            return Diagonal(s)
        end
        S_diag = vcat(diag.(S)...)

        # Perform approximate inference using the state-space approach.
        post_marginals_sde = TemporalGPs.approx_posterior_marginals(
            dtc, to_sde(f)(x, S_diag), vcat(y...), z_r, r_pred,
        )

        # Perform approximate inference naively.
        fx = f(x_naive, S_naive)
        f_post_approx = f | Stheno.PseudoObs(fx ← y_naive, f(z_naive, 1e-9))
        x_pred_naive = collect(TemporalGPs.RectilinearGrid(r_pred, t))
        post_marginals_naive = marginals(f_post_approx(x_pred_naive))

        approx_post_mean_sde = reshape(mean.(post_marginals_sde), length(r_pred), :);
        approx_post_mean_naive = reshape(mean.(post_marginals_naive), length(r_pred), :);
        @assert maximum(abs.(approx_post_mean_sde .- approx_post_mean_naive)) < 1e-7

        approx_post_std_sde = reshape(std.(post_marginals_sde), length(r_pred), :);
        approx_post_std_naive = reshape(std.(post_marginals_naive), length(r_pred), :);
        @assert maximum(abs.(approx_post_std_sde .- approx_post_std_naive)) < 1e-7

        # Generate plot. This should be roughly half a page width.
        width_inches = 2.8
        width_px = width_inches * Plots.DPI
        height_px = width_px

        plt_post = plot(
            xlabel="time",
            ylabel="space",
            legend=:bottomright,
            legendfontsize=9,
            background_color_legend=RGBA(1, 1, 1, 1),
            size=(width_px, height_px),
            grid = (0, 0),
            colorbar=nothing,
            xlims=(minimum(t), maximum(t)),
            ylim=(minimum(r_pred), maximum(r_pred)),
        );

        scatter!(plt_post, fill(0.0, length(z_r)), z_r;
            color=:black,
            label="Pseudo-Inputs",
        );
        scatter!(plt_post, last.(x_naive), first.(x_naive);
            markersize=0.1,
            color=:black,
            label="Observation Inputs",
        )
        savefig(plt_post, joinpath(first(save_dirs), "post_rectlinear.tikz"));

        # This image will have to be manually added onto the plot after it's been made.
        mean_min = minimum(approx_post_mean_sde);
        mean_max = maximum(approx_post_mean_sde);
        img = colorview(Gray, (approx_post_mean_sde .- mean_min) ./ (mean_max - mean_min));
        FileIO.save(joinpath(first(save_dirs), "post_rectlinear.jpg"), img);
    end

    let
        # Set seed for reproducibility.
        rng = MersenneTwister(123456)

        # Extract experiment parameters.
        θ = (σ = 0.92, λ_space = 1.1, λ_time = 0.83, S_min = 0.1)
        T = 100
        N_space_per_time = 10
        M_per_time = 10

        # Specify input locations. We place these irregularly through time.
        t = sort(rand(rng, T)) .* 50;
        t = RegularSpacing(0.0, 0.3, T);
        r = [rand(rng, N_space_per_time) * 10 for _ in 1:T];
        r_pred = collect(range(0.0, 10.0; length=100))
        x = TemporalGPs.RegularInTime(t, r);
        x_naive = collect(x);

        # Specify inducing inputs. They are regularly-spaced between 0 and 10.
        z_r = collect(range(0.0, 10.0; length=M_per_time));
        z_naive = collect(TemporalGPs.RectilinearGrid(z_r, t));

        # Specify observation noise.
        S = [Diagonal(rand(rng, N_space_per_time) .+ θ.S_min) for _ in 1:T];
        S_naive = vcat(diag.(S)...);

        # Generate data exactly / naively.
        k = separable_kernel(θ.σ, θ.λ_space, θ.λ_time);
        f = GP(k, GPC());

        # Hack together generation.
        N_per_slice = 5_000;
        ys_naive = map(1:div(length(x_naive), N_per_slice)) do j
            idx = N_per_slice * (j - 1) + 1:min(N_per_slice * j, length(x_naive))
            return rand(rng, f(x_naive[idx], S_naive[idx]))
        end;
        ys_naive = if rem(length(x_naive), N_per_slice) == 0
            ys_naive
        else
            idx = length(x_naive) - rem(length(x_naive), N_per_slice) + 1:length(x_naive)
            vcat(
                ys_naive,
                rand(rng, f(x_naive[idx], S_naive[idx]))
            )
        end;
        y_naive = vcat(ys_naive...);
        y = TemporalGPs.match_to(y_naive, x);

        # Perform approximate inference using the state-space approach.
        post_marginals_sde = TemporalGPs.approx_posterior_marginals(
            dtc, to_sde(f)(x, S_naive), y_naive, z_r, r_pred,
        )

        # Perform approximate inference naively.
        fx = f(x_naive, S_naive);
        f_post_approx = f | Stheno.PseudoObs(fx ← y_naive, f(z_naive, 1e-9));
        x_pred_naive = collect(TemporalGPs.RectilinearGrid(r_pred, t));
        post_marginals_naive = marginals(f_post_approx(x_pred_naive));

        approx_post_mean_sde = reshape(mean.(post_marginals_sde), length(r_pred), :);
        approx_post_mean_naive = reshape(mean.(post_marginals_naive), length(r_pred), :);
        @assert maximum(abs.(approx_post_mean_sde .- approx_post_mean_naive)) < 1e-7

        approx_post_std_sde = reshape(std.(post_marginals_sde), length(r_pred), :);
        approx_post_std_naive = reshape(std.(post_marginals_naive), length(r_pred), :);
        @assert maximum(abs.(approx_post_std_sde .- approx_post_std_naive)) < 1e-7

        # Generate plot. This should be roughly half a page width.
        width_inches = 2.8
        width_px = width_inches * Plots.DPI
        height_px = width_px

        plt_post = plot(
            xlabel="time",
            ylabel="space",
            legend=:bottomright,
            legendfontsize=9,
            background_color_legend=RGBA(1, 1, 1, 1),
            size=(width_px, height_px),
            grid = (0, 0),
            colorbar=nothing,
            xlims=(minimum(t), maximum(t)),
            ylim=(minimum(r_pred), maximum(r_pred)),
        );

        scatter!(plt_post, fill(0.0, length(z_r)), z_r;
            color=:black,
            label="Pseudo-Inputs",
        );
        scatter!(plt_post, last.(x_naive), first.(x_naive);
            markersize=0.1,
            color=:black,
            label="Observation Inputs",
        );
        savefig(plt_post, joinpath(first(save_dirs), "post_irregular.tikz"));

        # This image will have to be manually added onto the plot after it's been made.
        mean_min = minimum(approx_post_mean_sde);
        mean_max = maximum(approx_post_mean_sde);
        img = colorview(Gray, (approx_post_mean_sde .- mean_min) ./ (mean_max - mean_min));
        FileIO.save(joinpath(first(save_dirs), "post_irregular.jpg"), img);
    end
end


text:
using Test
using UnitTestDesign

res1 = UnitTestDesign.full_factorial([2, 3, 4, 2, 2], x -> false)
@test size(res1, 2) == prod([2, 3, 4, 2, 2])
@test UnitTestDesign.test_coverage(res1, [2, 3, 4, 2, 2], 5).finish == 0

omit = x -> x[3] == 4 && x[5] == 1
res2 = UnitTestDesign.full_factorial([2, 3, 4, 2, 2], omit)
@test size(res2, 2) == prod([2, 3, 4, 2, 2]) - prod([2, 3, 2])
@test UnitTestDesign.test_coverage(res2, [2, 3, 4, 2, 2], 5).finish == 12


text:
include("1.3.1 P(x)^p.jl")
include("1.3.2 P(x) Q(x)^p.jl")
include("1.3.3 Miscellaneous algebraic functions.jl")
include("1.3.4 Normalizing algebraic functions.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.3 Miscellaneous/.jl


text:
function foo()
    lines = readlines(open("data/7.txt", "r"))
    basework = 60
    workers = 5
    deps = Dict{Char, Array{Char, 1}}()
    total = Set{Char}()
    for line in lines
        a = line[6]
        b = line[37]
        if !(a in keys(deps))
            deps[a] = []
        end
        push!(total, a)
        push!(total, b)
        push!(deps[a], b)
    end
    working = Dict{Char, Int}()
    total_time = 0
    while !(isempty(total) && isempty(working))
        j = Set(vcat(values(deps)...))
        avail = sort(collect(setdiff(total, j)))
        for c in avail[1:min(end, workers-length(working))]
            working[c] = c - 'A' + 1 + basework
            delete!(total, c)
        end
        worktime = minimum(values(working))
        @assert worktime > 0
        total_time += worktime
        for (char, time) in working
            time -= worktime
            @assert time >= 0
            if time == 0
                delete!(deps, char)
                delete!(working, char)
            else
                working[char] = time
            end
        end
    end
    return total_time
end

println(foo())


text:
using ITensors

#
# Qubit site type
#

ITensors.space(::SiteType"Qubit") = 2

ITensors.state(::SiteType"Qubit", ::StateName"0") = 1

ITensors.state(::SiteType"Qubit", ::StateName"1") = 2

#
# 1-Qubit gates
#

op_matrix(s::String) = op_matrix(OpName(s))

op_matrix(::OpName"Id") =
  [1 0
   0 1]

op_matrix(::OpName"I") =
  op_matrix("Id")

op_matrix(::OpName"X") =
  [0 1
   1 0]

op_matrix(::OpName"σx") =
  op_matrix("X")

op_matrix(::OpName"σ1") =
  op_matrix("X")

op_matrix(::OpName"√NOT") =
  [(1+im)/2 (1-im)/2
   (1-im)/2 (1+im)/2]

op_matrix(::OpName"√X") =
  op_matrix("√NOT")

op_matrix(::OpName"Y") =
  [ 0 -im
   im   0]

op_matrix(::OpName"σy") =
  op_matrix("Y")

op_matrix(::OpName"σ2") =
  op_matrix("Y")

op_matrix(::OpName"iY") =
  [ 0 1
   -1 0]

op_matrix(::OpName"iσy") =
  op_matrix("iY")

op_matrix(::OpName"iσ2") =
  op_matrix("iY")

# Rϕ with ϕ = π
op_matrix(::OpName"Z") =
  [1  0
   0 -1]

op_matrix(::OpName"σz") =
  op_matrix("Z")

op_matrix(::OpName"σ3") =
  op_matrix("Z")

op_matrix(::OpName"H") =
  [1/sqrt(2) 1/sqrt(2)
   1/sqrt(2) -1/sqrt(2)]

# Rϕ with ϕ = π/2
op_matrix(::OpName"Phase") =
  [1  0
   0 im]

op_matrix(::OpName"P") =
  op_matrix("Phase")

op_matrix(::OpName"S") =
  op_matrix("Phase")

# Rϕ with ϕ = π/4
op_matrix(::OpName"π/8") =
  [1  0
   0  1/sqrt(2) + im/sqrt(2)]

op_matrix(::OpName"T") =
  op_matrix("π/8")

# Rotation around X-axis
op_matrix(::OpName"Rx"; θ::Number) =
  [    cos(θ/2)  -im*sin(θ/2)
   -im*sin(θ/2)      cos(θ/2)]

# Rotation around Y-axis
op_matrix(::OpName"Ry"; θ::Number) =
  [cos(θ/2) -sin(θ/2)
   sin(θ/2)  cos(θ/2)]

# Rotation around Z-axis
op_matrix(::OpName"Rz"; ϕ::Number) =
  [1         0
   0 exp(im*ϕ)]

#  [exp(-im*ϕ/2)           0
#   0            exp(im*ϕ/2)]

# Rotation around generic axis n̂
op_matrix(::OpName"Rn";
    θ::Real, ϕ::Real, λ::Real) =
  [          cos(θ/2)    -exp(im*λ)*sin(θ/2)
   exp(im*ϕ)*sin(θ/2) exp(im*(ϕ+λ))*cos(θ/2)]

op_matrix(::OpName"Rn̂"; kwargs...) =
  op_matrix("Rn"; kwargs...)

#
# 2-Qubit gates
#

op_matrix(::OpName"CNOT") =
  [1 0 0 0
   0 1 0 0
   0 0 0 1
   0 0 1 0]

op_matrix(::OpName"CX") =
  op_matrix("CNOT")

op_matrix(::OpName"CY") =
  [1 0  0   0
   0 1  0   0
   0 0  0 -im
   0 0 im   0]

op_matrix(::OpName"CZ") =
  [1 0 0  0
   0 1 0  0
   0 0 1  0
   0 0 0 -1]

op_matrix(::OpName"SWAP") =
  [1 0 0 0
   0 0 1 0
   0 1 0 0
   0 0 0 1]

op_matrix(::OpName"√SWAP") =
  [1        0        0 0
   0 (1+im)/2 (1-im)/2 0
   0 (1-im)/2 (1+im)/2 0
   0        0        0 1]

# Ising (XX) coupling gate
op_matrix(::OpName"XX"; ϕ::Number) =
  [    cos(ϕ)          0          0 -im*sin(ϕ)
            0     cos(ϕ) -im*sin(ϕ)          0
            0 -im*sin(ϕ)     cos(ϕ)          0
   -im*sin(ϕ)          0          0     cos(ϕ)]

# Ising (YY) coupling gate
op_matrix(::OpName"YY"; ϕ::Number) =
  [    cos(ϕ)          0          0  im*sin(ϕ)
            0     cos(ϕ) -im*sin(ϕ)          0
            0 -im*sin(ϕ)     cos(ϕ)          0
    im*sin(ϕ)          0          0     cos(ϕ)]

# Ising (ZZ) coupling gate
op_matrix(::OpName"ZZ"; ϕ::Number) =
  [    exp(-im*ϕ)         0         0          0
            0     exp(im*ϕ)         0          0
            0             0 exp(im*ϕ)          0
            0             0         0 exp(-im*ϕ)]

#
# 3-Qubit gates
#

op_matrix(::OpName"Toffoli") =
  [1 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0
   0 0 1 0 0 0 0 0
   0 0 0 1 0 0 0 0
   0 0 0 0 1 0 0 0
   0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 1
   0 0 0 0 0 0 1 0]

op_matrix(::OpName"CCNOT") =
  op_matrix("Toffoli")

op_matrix(::OpName"CCX") =
  op_matrix("Toffoli")

op_matrix(::OpName"TOFF") =
  op_matrix("Toffoli")

op_matrix(::OpName"Fredkin") =
  [1 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0
   0 0 1 0 0 0 0 0
   0 0 0 1 0 0 0 0
   0 0 0 0 1 0 0 0
   0 0 0 0 0 0 1 0
   0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 1]

op_matrix(::OpName"CSWAP") =
  op_matrix("Fredkin")

op_matrix(::OpName"CS") =
  op_matrix("Fredkin")

#
# 4-Qubit gates
#

op_matrix(::OpName"CCCNOT") =
  [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]

op_matrix(::OpName"randn"; dim) =
    randn(dim)

op_matrix(::OpName"noise"; dim) =
  randn(dim)

op_matrix(gn::OpName; kwargs...) = error("Gate $gn not implemented.")

function op_matrix(on::OpName, s::Index...; kwargs...)
  rs = reverse(s)
  return itensor(op_matrix(on; kwargs...), prime.(rs)..., dag.(rs)...)
end

op_matrix(gn::OpName"randn", s::Index...) =
  itensor(op_matrix(gn; dim = dim(s) * dim(s)), prime.(s)..., dag.(s)...)

op_matrix(gn::OpName"noise", s::Index...; krausind = Index(2, "kraus")) =
  itensor(op_matrix(gn; dim = dim(s) * dim(s) * dim(krausind)), prime.(s)..., dag.(s)..., krausind)

op_matrix(gn::String, s::Index...; kwargs...) =
  op_matrix(OpName(gn), s...; kwargs...)

op_matrix(gn::String, s::Vector{<:Index}, ns::Int...; kwargs...) =
  op_matrix(OpName(gn), s[[ns...]]...; kwargs...)

ITensors.op(gn::OpName, ::SiteType"Qubit", s::Index...; kwargs...) =
  op_matrix(gn, s...; kwargs...)


text:
import ..UncertainValues:
    AbstractScalarPopulation

import Base.rand
import StatsBase.sample

function resample(p::AbstractScalarPopulation)
    rand(p)
end

function resample(p::AbstractScalarPopulation, n::Int)
    rand(p, n)
end

function resample(p::AbstractScalarPopulation, constraint::SamplingConstraint)
    rand(constrain(p, constraint))
end

function resample(p::AbstractScalarPopulation, constraint::SamplingConstraint, n::Int)
    rand(constrain(p, constraint), n)
end

constraints = [
    :(NoConstraint),
    :(TruncateLowerQuantile),
    :(TruncateUpperQuantile),
    :(TruncateQuantiles),
    :(TruncateMaximum),
    :(TruncateMinimum),
    :(TruncateRange),
    :(TruncateStd)
]

for constraint in constraints
    funcs = quote
        function resample(p::AbstractScalarPopulation{T, PW}, constraint::$(constraint)) where {T, PW}
            rand(constrain(p, constraint))
        end

        function resample(p::AbstractScalarPopulation{T, PW}, constraint::$(constraint), n::Int) where {T, PW}
            rand(constrain(p, constraint), n)
        end
    end
    eval(funcs)
end


text:
# # Three Dimensions
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/solid_body_rotation.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/solid_body_rotation.ipynb)
#
# Simulate the trajectory of an individual point, first in a perfectly circular flow (a.k.a. solid body rotation). Then add a convergent term to obtain a spiraling trajectory, and a constant vertical velocity for the third dimension. These simple flow configurations can be thought of as idealized models e.g. ocean meso-scale eddies.
#
# For additional documentation e.g. see :
# [1](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),
# [2](https://JuliaClimate.github.io/MeshArrays.jl/dev/),
# [3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),
# [4](https://en.wikipedia.org/wiki/Displacement_(vector))
#
# Exercise examples:
# - make the sinking velocity decrease with time (hint: it increases in the original notebook)
# - change the number of times the particle goes around the origin
#
# ![solid body rotation](https://github.com/JuliaClimate/IndividualDisplacements.jl/raw/master/examples/figs/SolidBodyRotation.gif)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # 1 Problem Configuration
#
# Here we set up software, grid, flow fields, initial conditions.
#
# ### 1.1 Import Software

using IndividualDisplacements, DataFrames
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/flow_fields.jl"))

#nb # %% {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### 1.2  Flow Fields
#
# The `simple_flow_field` function (defined in `helper_functions.jl`) defines a simple
# three-dimensional flow field. Exercise: locate `simple_flow_field` and modify the
# flow field parameters (e.g. intensity and sign of the convergent term).

np,nz=16,4 #gridded domain size (horizontal and vertical)

u,v,w=solid_body_rotation(np,nz) #staggered velocity arrays

𝐹=FlowFields(u,u,v,v,0*w,1*w,[0,19.95*2*pi]); #FlowFields data structure

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### 1.3 Initialize Individuals
#
# Let's just set up one individual at [np*1/3,np*1/3,nz*1/3] in the three-dimensional
# space where the flow fields have been configured

(x,y,z)=(np*1/3,np*1/3,nz*1/3)

𝐼=Individuals(𝐹,x,y,z)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### 1.4 A Closer Look (optional)
#
# The above `Individuals` constructor wraps up 𝐹, the initial position, and other needed components
# within 𝐼. **At this point, you can either jump to section 2 or read through this section**
# to learn more about how the details as needed e.g. if you wanted to overide default options
# that were selected for you by the section 1.3 constructor.
#
# Initial position is

📌=[x,y,z]

# and the data structure ([DataFrame](http://juliadata.github.io/DataFrames.jl/stable/))
# to record properties along the individual's path accordingly.

🔴 = DataFrame(ID=Int[], x=Float64[], y=Float64[], z=Float64[], t=Float64[])

# It is the postprocessing function's responsibility to provide the record. It is thus
# important that this intermediary (`postproc`) be consistent with the solver setup (`sol`)
# and the expected record format (`🔴`).

function postproc(sol,𝐹::FlowFields;id=missing,𝑇=missing)
    df=postprocess_xy(sol,𝐹,id=id,𝑇=𝑇)
    #add third coordinate
    z=sol[3,:]
    df.z=z[:]
    return df
end

# The velocity function `🚄` relies only on flow fields obtained from
# `𝐹` (which is defined above) to interpolate velocity at the specified
# space-time position (e.g. those of individuals).

🚄 = dxdt!

# Now that every thing needed to carry out the computation is in place,
# we wrap up the problem configuration in a struct (`Individuals`) which
# links to the initial positions, flow fields, etc. all that will be
# necessary to compute trajectories over time (`∫!(𝐼,𝑇)`).

#assemble as a NamedTuple:
I=(position=📌,record=🔴,velocity=🚄,
postprocessing=postproc,parameters=𝐹)

#construct Individuals from NamedTuple:
𝐼=Individuals(I)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 2 Trajectory Simulations
#
# The `∫!` function call below returns the final positions & updates `𝐼.📌` accordingly. It also records properties observed along the trajectory in `𝐼.🔴`.
# Simple methods to visualize the individual trajectory (plot or movie) are provided at the end.

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### 2.1 Compute Trajectories

𝑇=(0.0,𝐼.𝑃.𝑇[2])
∫!(𝐼,𝑇)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### 2.2 Visualize Trajectories
#
# - define `myplot` convenience function
# - single plot example using `myplot`
# - (generate animation using `myplot`)

#md p=dirname(pathof(IndividualDisplacements))
#md include(joinpath(p,"../examples/recipes_plots.jl"));
#md nt=length(𝐼.🔴.x)

#md myplot(i)=plot(𝐼.🔴.x[1:i],𝐼.🔴.y[1:i],𝐼.🔴.z[1:i],linewidth=2,arrow = 2,
#md     title="Solid body rotation / Spiral example",leg=false,
#md     xaxis="x",yaxis="y",zaxis="z",xlims=(0,np),ylims=(0,np));

#nb # %% {"slideshow": {"slide_type": "subslide"}}
# Single plot example:

#md plt=myplot(nt)
#md scatter!(plt,[📌[1]],[📌[2]],[📌[3]])
#md #scatter!(plt,[𝐼.🔴.x[end]],[𝐼.🔴.y[end]],[𝐼.🔴.z[end]])
#md scatter!(plt,[𝐼.📌[1]],[𝐼.📌[2]],[𝐼.📌[3]])

#nb # %% {"slideshow": {"slide_type": "subslide"}}
# Animation example:
#
# ```
# p=Int(ceil(nt/100))
# anim = @animate for i ∈ 1:p:nt
#     myplot(i)
# end
#
# pth=tempdir()*"/"
# gif(anim, pth*"SolidBodyRotation.gif", fps = 15)
# ```


text:
"""
FaultTree
"""

export ftkofn, ftand, ftor, ftnot

import Base

function Base.:*(x::AbstractFaultTreeNode, y::AbstractFaultTreeNode)
    ftand(x, y)
end

function Base.:&(x::AbstractFaultTreeNode, y::AbstractFaultTreeNode)
    ftand(x, y)
end

function Base.:+(x::AbstractFaultTreeNode, y::AbstractFaultTreeNode)
    ftor(x, y)
end

function Base.:|(x::AbstractFaultTreeNode, y::AbstractFaultTreeNode)
    ftor(x, y)
end

function Base.:!(x::AbstractFaultTreeNode)
    ftnot(x)
end

function Base.:~(x::AbstractFaultTreeNode)
    ftnot(x)
end

function ftand(x::Vararg{AbstractFaultTreeNode})
    s = union([y.params for y = x]...)
    args = [y for y = x]
    FaultTreeOperation(s, :AND, args)
end

function ftor(x::Vararg{AbstractFaultTreeNode})
    s = union([y.params for y = x]...)
    args = [y for y = x]
    FaultTreeOperation(s, :OR, args)
end

function ftnot(x::AbstractFaultTreeNode)
    FaultTreeOperation(x.params, :NOT, [x])
end

function ftkofn(k::Int, x::Vararg{AbstractFaultTreeNode})
    s = union([y.params for y = x]...)
    args = [y for y = x]
    FaultTreeKoutofN(s, :KofN, k, args)
end


text:
using Compat: mapreduce

function _treepast(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    branches = BL[]
    nodes = NL[node]
    while hasinbound(tree, node)
        inbound = getinbound(tree, node)
        branches = push!(branches, inbound)
        node = src(tree, inbound)
        nodes = push!(nodes, node)
    end
    return branches, nodes
end

function _treefuture(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    branches = BL[]
    nodestoprocess = NL[node]
    nodesprocessed = NL[]
    while !isempty(nodestoprocess)
        nextnode = pop!(nodestoprocess)
        push!(nodesprocessed, nextnode)
        outbounds = getoutbounds(tree, nextnode)
        append!(branches, outbounds)
        children = map(branch -> dst(tree, branch), outbounds)
        append!(nodestoprocess, children)
    end
    return branches, nodesprocessed
end

"""
    branchhistory(tree::AbstractTree, node)

Find the branch route between a node on a tree and its root
"""
function branchhistory(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return _treepast(tree, node)[1]
end

"""
    nodehistory(tree::AbstractTree, node)

Find the node route between a node on a tree and its root
"""
function nodehistory(tree::T, node::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    return _treepast(tree, node)[2]
end

"""
    branchroute(tree::AbstractTree, node1, node2)

Find the branch route between two nodes on a tree
"""
function branchroute(tree::T, node1::NL, node2::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    branches1, nodes1 = _treepast(tree, node1)
    branches2, nodes2 = _treepast(tree, node2)
    nodes1[end] == nodes2[end] ||
        return error("No route between nodes")
    common = branches1 ∩ branches2
    return append!(filter(b -> b ∉ common, branches1),
                   filter(b -> b ∉ common, reverse(branches2)))
end

"""
    noderoute(tree::AbstractTree, node1, node2)

Find the node route between two nodes on a tree
"""
function noderoute(tree::T, node1::NL, node2::NL) where {NL, BL, T <: AbstractTree{NL, BL}}
    branches1, nodes1 = _treepast(tree, node1)
    branches2, nodes2 = _treepast(tree, node2)
    nodes1[end] == nodes2[end] ||
        return error("No route between nodes")
    common = nodes1[end]
    while min(length(nodes1), length(nodes2)) > 0 && nodes1[end] == nodes2[end]
        common = nodes1[end]
        pop!(nodes1)
        pop!(nodes2)
    end
    push!(nodes1, common)
    return append!(nodes1, reverse(nodes2))
end

"""
    distance(tree::AbstractTree, node1, node2)

Distance between two nodes on a tree
"""
function distance(tree::AbstractTree, node1, node2)
    branches = branchroute(tree, node1, node2)
    return mapreduce(branch -> getlength(tree, branch), +, branches;
    init = 0.0)
end

"""
    distances(tree::AbstractTree)

Pairwise distances between all leaf nodes on a tree
"""
function distances(tree::AbstractTree)
    leaves = nodenamefilter(isleaf, tree)
    return [distance(tree, i, j) for i in leaves, j in leaves]
end

"""
    height(tree::AbstractTree, node)

Height of a node of the tree above the root
"""
function heighttoroot(tree::AbstractTree, node)
    return mapreduce(branch -> getlength(tree, branch), +,
                     branchhistory(tree, node); init = 0.0)
end

"""
    heights(tree::AbstractTree)

Height of all of the leaves of the tree above the root
"""
function heightstoroot(tree::AbstractTree)
    return [heighttoroot(tree, i) for i in nodenamefilter(isleaf, tree)]
end


text:
using Test
using BetweenFlags.PerFlagFuncs

@testset "split_by_consecutives edge case" begin
  @test [1] == split_by_consecutives([1])[1]
end

@testset "merge_odd_even" begin
  a = [1, 3, 5]
  b = [2, 4, 6]
  c = merge_even_odd(a, b)
  @test all([x==y for (x, y) in zip(c, [1, 2, 3, 4, 5, 6])])
end

function get_alternating_consecutive_filter(A, B)
  dupes = [x for x in A for y in B if x==y]
  A = [x for x in A if !any([x == y for y in dupes])]
  B = [x for x in B if !any([x == y for y in dupes])]
  sort!(A)
  sort!(B)
  (C, D) = get_alternating_consecutive_vector(A, B)
  M = merge_even_odd(C, D)
  return M
end

@testset "get_alternating_consecutive_vector" begin
  # A few general cases:
  A = [11, 5, 6, 20, 2, 1, 4, 15, 12, 7, 18, 14]
  B = [5, 3, 16, 17, 13, 9, 20, 14, 1, 11]
  M_correct = [2, 3, 4, 9, 6, 9, 7, 13, 12, 13, 15, 16]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  A = [1, 3, 4, 5, 6, 7, 10, 11, 15, 16, 17]
  B = [2, 8, 9]
  M_correct = [1, 2, 3, 8, 4, 8, 5, 9, 6, 9]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  A = [5, 7, 14, 20]
  B = [1, 4, 9, 13, 17]
  M_correct = [5, 9, 7, 9, 14, 17]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  # Empty result case:
  A = [8, 10]
  B = [4]
  M_correct = Int64[]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  # Some edge cases:
  A = Int64[]
  B = [4]
  M_correct = Int64[]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  A = [1]
  B = Int64[]
  M_correct = Int64[]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])

  A = Int64[]
  B = Int64[]
  M_correct = Int64[]
  M = get_alternating_consecutive_filter(A, B)
  @test all([x==y for (x,y) in zip(M, M_correct)])
end

@testset "Test generator" begin
  # This function is a bit difficult to test automatically.
  # So, this is a test generator, and its test results are
  # trivially satisfied.

  # print_IO = false
  N_A = Base.rand(10:50)[1]
  N_B = Base.rand(10:50)[1]
  A = unique([Base.rand(1:20)[1] for x in 1:Base.rand(1:N_A)])
  B = unique([Base.rand(1:20)[1] for x in 1:Base.rand(1:N_B)])
  dupes = [x for x in A for y in B if x==y]
  A = [x for x in A if !any([x == y for y in dupes])]
  B = [x for x in B if !any([x == y for y in dupes])]
  sort!(A)
  sort!(B)
  # if print_IO
  #   print("\n------------------------- Original\n")
  #   print("A = ", A,"\n")
  #   print("B = ", B,"\n")
  # end
  (C, D) = get_alternating_consecutive_vector(A, B)
  M = merge_even_odd(C, D)
  @test all([x==y for (x,y) in zip(M, M)])
  # if print_IO
  #   print("M = ", M,"\n")
  #   print("\n-------------------------\n")
  # end
end


text:
export
    edlimb,
    edterm,
    ekacec,
    ekaced,
    ekacei,
    ekaclc,
    ekacld,
    ekacli,
    ekappr,
    ekbseg,
    ekccnt,
    ekcii,
    ekcls,
    ekdelr,
    ekffld,
    ekfind,
    ekgc,
    ekgd,
    ekgi,
    ekinsr,
    ekifld,
    eklef,
    eknelt,
    eknseg,
    ekntab,
    ekopn,
    ekopr,
    ekops,
    ekopw,
    ekpsel,
    ekrcec,
    ekrced,
    ekrcei,
    ekssum,
    ektnam,
    ekucec,
    ekuced,
    ekucei,
    ekuef,
    el2cgv,
    elemc,
    elemd,
    elemi,
    eqncpv,
    eqstr,
    esrchc,
    et2lst,
    et2utc,
    etcal,
    eul2m,
    eul2xf,
    expool

"""
    edlimb(a, b, c, viewpt)

Find the limb of a triaxial ellipsoid, viewed from a specified point.

### Arguments ###

- `a`: Length of ellipsoid semi-axis lying on the x-axis
- `b`: Length of ellipsoid semi-axis lying on the y-axis
- `c`: Length of ellipsoid semi-axis lying on the z-axis
- `viewpt`: Location of viewing point

### Output ###

Returns the limb of the ellipsoid as seen from the viewing point.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/edlimb_c.html)
"""
function edlimb(a, b, c, viewpt)
    @checkdims 3 viewpt
    limb = Ref{Ellipse}()
    ccall((:edlimb_c, libcspice), Cvoid,
           (SpiceDouble, SpiceDouble, SpiceDouble, Ref{SpiceDouble}, Ref{Ellipse}),
           a, b, c, viewpt, limb)
    handleerror()
    limb[]
end

"""
    edterm(trmtyp, source, target, et, fixref, abcorr, obsrvr, npts)

Compute a set of points on the umbral or penumbral terminator of a specified
target body, where the target shape is modeled as an ellipsoid.

### Arguments ###

- `trmtyp`: Terminator type
- `source`: Light source
- `target`: Target body
- `et`: Observation epoch
- `fixref`: Body-fixed frame associated with target
- `abcorr`: Aberration correction
- `obsrvr`: Observer
- `npts`: Number of points in terminator set

### Output ###

- `trgepc`: Epoch associated with target center
- `obspos`: Position of observer in body-fixed frame
- `trmpts`: Terminator point set

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/edterm_c.html)
"""
function edterm(trmtyp, source, target, et, fixref, abcorr, obsrvr, npts)
    trgepc = Ref{SpiceDouble}()
    obspos = Array{SpiceDouble}(undef, 3)
    trmpts = Array{SpiceDouble}(undef, 3, npts)
    ccall((:edterm_c, libcspice), Cvoid,
          (Cstring, Cstring, Cstring, SpiceDouble, Cstring, Cstring, Cstring, SpiceInt,
           Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          trmtyp, source, target, et, fixref, abcorr, obsrvr, npts, trgepc, obspos, trmpts)
    handleerror()
    trgepc[], obspos, trmpts
end

"""
    ekacec(handle, segno, recno, column, cvals, isnull)

Add data to a character column in a specified EK record.

### Arguments ###

- `handle`: EK file handle
- `segno`: Index of segment containing record
- `recno`: Record to which data is to be added
- `column`: Column name
- `cvals`: Character values to add to column
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekacec_c.html)
"""
function ekacec(handle, segno, recno, column, cvals, isnull)
    cvals_, nvals, vallen = chararray(cvals)
    ccall((:ekacec_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, vallen, cvals_, isnull)
    handleerror()
end

"""
    ekaced(handle, segno, recno, column, dvals, isnull)

Add data to an double precision column in a specified EK record.

### Arguments ###

- `handle`: EK file handle
- `segno`: Index of segment containing record
- `recno`: Record to which data is to be added
- `column`: Column name
- `dvals`: Double precision values to add to column
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekaced_c.html)
"""
function ekaced(handle, segno, recno, column, dvals, isnull)
    nvals = length(dvals)
    ccall((:ekaced_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt, Ref{SpiceDouble}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, dvals, isnull)
    handleerror()
end

"""
    ekacei(handle, segno, recno, column, ivals, isnull)

Add data to an integer column in a specified EK record.

### Arguments ###

- `handle`: EK file handle
- `segno`: Index of segment containing record
- `recno`: Record to which data is to be added
- `column`: Column name
- `ivals`: Integer values to add to column
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekacei_c.html)
"""
function ekacei(handle, segno, recno, column, ivals, isnull)
    nvals = length(ivals)
    ivals_ = SpiceInt.(ivals)
    ccall((:ekacei_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt, Ref{SpiceInt}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, ivals_, isnull)
    handleerror()
end

"""
    ekaclc(handle, segno, column, cvals, nlflgs, rcptrs)

Add an entire character column to an EK segment.

### Arguments ###

- `handle`: EK file handle.
- `segno`: Number of segment to add column to.
- `column`: Column name.
- `cvals`: Character values to add to column.
- `nlflgs`: Array of null flags for column entries.
- `rcptrs`: Record pointers for segment.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekaclc_c.html)
"""
function ekaclc(handle, segno, column, cvals, nlflgs, rcptrs)
    cvals_, nrows, vallen = chararray(cvals)
    nlflgs_ = SpiceBoolean.(nlflgs)
    @checkdims nrows nlflgs rcptrs
    entszs = SpiceInt.(length.(cvals))
    entszs[nlflgs] .= 0
    wkindx = Array{SpiceInt}(undef, nrows)
    ccall((:ekaclc_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Cstring, SpiceInt, Ref{SpiceChar}, Ref{SpiceInt}, Ref{SpiceBoolean},
           Ref{SpiceInt}, Ref{SpiceInt}),
          handle, segno - 1, column, vallen, cvals_, entszs, nlflgs_, rcptrs, wkindx)
    handleerror()
end

"""
    ekacld(handle, segno, column, dvals, nlflgs, rcptrs)

Add an entire double precision column to an EK segment.

### Arguments ###

- `handle`: EK file handle
- `segno`: Number of segment to add column to
- `column`: Column name
- `dvals`: Double precision values to add to column
- `nlflgs`: Array of null flags for column entries
- `rcptrs`: Record pointers for segment

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekacld_c.html)
"""
function ekacld(handle, segno, column, dvals, nlflgs, rcptrs)
    nrows = length(dvals)
    entszs = SpiceInt.(length.(dvals))
    entszs[nlflgs] .= 0
    dvals_ = array_to_cmatrix(dvals)
    nlflgs_ = SpiceBoolean.(nlflgs)
    @checkdims nrows nlflgs rcptrs
    wkindx = Array{SpiceInt}(undef, nrows)
    ccall((:ekacld_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Cstring, Ref{SpiceDouble}, Ref{SpiceInt}, Ref{SpiceBoolean},
           Ref{SpiceInt}, Ref{SpiceInt}),
          handle, segno - 1, column, dvals_, entszs, nlflgs_, rcptrs, wkindx)
    handleerror()
end

"""
    ekacli(handle, segno, column, ivals, nlflgs, rcptrs)

Add an entire integer column to an EK segment.

### Arguments ###

- `handle`: EK file handle
- `segno`: Number of segment to add column to
- `column`: Column name
- `ivals`: Integer values to add to column
- `nlflgs`: Array of null flags for column entries
- `rcptrs`: Record pointers for segment

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekacli_c.html)
"""
function ekacli(handle, segno, column, ivals, nlflgs, rcptrs)
    nrows = length(ivals)
    entszs = SpiceInt.(length.(ivals))
    entszs[nlflgs] .= 0
    ivals_ = array_to_cmatrix(ivals)
    nlflgs_ = SpiceBoolean.(nlflgs)
    @checkdims nrows nlflgs rcptrs
    wkindx = Array{SpiceInt}(undef, nrows)
    ccall((:ekacli_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Cstring, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceBoolean},
           Ref{SpiceInt}, Ref{SpiceInt}),
          handle, segno - 1, column, ivals_, entszs, nlflgs_, rcptrs, wkindx)
    handleerror()
end

"""
    ekappr(handle, segno)

Append a new, empty record at the end of a specified E-kernel segment.

### Arguments ###

- `handle`: File handle
- `segno`: Segment number

### Output ###

Returns the number of appended record.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekappr_c.html)
"""
function ekappr(handle, segno)
    recno = Ref{SpiceInt}()
    ccall((:ekappr_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{SpiceInt}),
          handle, segno - 1, recno)
    handleerror()
    recno[] + 1
end

"""
    ekbseg(handle, tabnam, cnames, decls)

Start a new segment in an E-kernel.

### Arguments ###

- `handle`: File handle
- `tabnam`: Table name
- `cnames`: Names of columns
- `decls `: Declarations of columns

### Output ###

Returns the segment number.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekbseg_c.html)
"""
function ekbseg(handle, tabnam, cnames, decls)
    cnames_, ncols, cnmlen = chararray(cnames)
    @checkdims ncols decls
    decls_, _, declen = chararray(decls)
    segno = Ref{SpiceInt}()
    ccall((:ekbseg_c, libcspice), Cvoid,
          (SpiceInt, Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, SpiceInt, Ref{SpiceChar},
           Ref{SpiceInt}),
          handle, tabnam, ncols, cnmlen, cnames_, declen, decls_, segno)
    handleerror()
    segno[] + 1
end

"""
    ekccnt(table)

Return the number of distinct columns in a specified, currently loaded table

### Arguments ###

- `table`: Name of table

### Output ###

Returns the count of distinct, currently loaded columns.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekccnt_c.html)
"""
function ekccnt(table)
    ccount = Ref{SpiceInt}()
    ccall((:ekccnt_c, libcspice), Cvoid, (Cstring, Ref{SpiceInt}), table, ccount)
    handleerror()
    Int(ccount[])
end

"""
    ekcii(table, cindex, lenout=256)

Return attribute information about a column belonging to a loaded EK table, specifying the column
by table and index.

### Arguments ###

- `table`: Name of table containing column
- `cindex`: Index of column whose attributes are to be found
- `lenout`: Maximum allowed length of column name (default: 256)

### Output ###

- `column`: Name of column
- `attdsc`: Column attribute descriptor

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekcii_c.html)
"""
function ekcii(table, cindex, lenout=256)
    column = Array{SpiceChar}(undef, lenout)
    attdsc = Ref{EKAttDsc}()
    ccall((:ekcii_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{EKAttDsc}),
          table, cindex - 1, lenout, column, attdsc)
    handleerror()
    chararray_to_string(column), attdsc[]
end

"""
    ekcls(handle)

Close an E-kernel.

### Arguments ###

- `handle`: EK file handle

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekcls_c.html)
"""
function ekcls(handle)
    ccall((:ekcls_c, libcspice), Cvoid, (SpiceInt,), handle)
    handleerror()
end

"""
    ekdelr(handle, segno, recno)

Delete a specified record from a specified E-kernel segment.

### Arguments ###

- `handle`: File handle
- `segno`: Segment number
- `recno`: Record number

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekdelr_c.html)
"""
function ekdelr(handle, segno, recno)
    ccall((:ekdelr_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt),
          handle, segno - 1, recno - 1)
    handleerror()
end

"""
    ekffld(handle, segno, rcptrs)

Complete a fast write operation on a new E-kernel segment.

### Arguments ###

- `handle`: File handle
- `segno`: Segment number
- `rcptrs`: Record pointers

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekffld_c.html)
"""
function ekffld(handle, segno, rcptrs)
    ccall((:ekffld_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{SpiceInt}),
          handle, segno - 1, rcptrs)
    handleerror()
end

"""
    ekfind(query, lenout=256)

Find E-kernel data that satisfy a set of constraints.

### Arguments ###

- `query`: Query specifying data to be found.
- `lenout`: Declared length of output error message string (default: 256)

### Output ###

Returns the number of matching rows.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekfind_c.html)
"""
function ekfind(query, lenout=256)
    nmrows = Ref{SpiceInt}()
    error = Ref{SpiceBoolean}()
    errmsg = Array{SpiceChar}(undef, lenout)
    ccall((:ekfind_c, libcspice), Cvoid,
          (Cstring, SpiceInt, Ref{SpiceInt}, Ref{SpiceBoolean}, Ref{SpiceChar}),
          query, lenout, nmrows, error, errmsg)
    handleerror()
    if Bool(error[])
        throw(SpiceError(chararray_to_string(errmsg)))
    end
    Int(nmrows[])
end

"""
    ekgc(selidx, row, elment, lenout=256)

Return an element of an entry in a column of character type in a specified row.

### Arguments ###

- `selidx`: Index of parent column in SELECT clause
- `row`: Row to fetch from
- `elment`: Index of element, within column entry, to fetch
- `lenout`: Maximum length of column element (default: 256)

### Output ###

Returns the character string element of column entry or `missing` if it was null
or `nothing` if the column was not found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekgc_c.html)
"""
function ekgc(selidx, row, elment, lenout=256)
    cdata = Array{SpiceChar}(undef, lenout)
    null = Ref{SpiceBoolean}()
    found = Ref{SpiceBoolean}()
    ccall((:ekgc_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceBoolean}, Ref{SpiceBoolean}),
          selidx - 1, row - 1, elment - 1, lenout, cdata, null, found)
    handleerror()
    Bool(found[]) || return nothing
    Bool(null[]) && return missing
    chararray_to_string(cdata)
end

"""
    ekgd(selidx, row, element)

Return an element of an entry in a column of double precision type in a specified row.

### Arguments ###

- `selidx`: Index of parent column in SELECT clause
- `row`: Row to fetch from
- `elment`: Index of element, within column entry, to fetch

### Output ###

Returns the double precision element of column entry or `missing` if it was null
or `nothing` if the column was not found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekgd_c.html)
"""
function ekgd(selidx, row, elment)
    null = Ref{SpiceBoolean}()
    found = Ref{SpiceBoolean}()
    ddata = Ref{SpiceDouble}()
    ccall((:ekgd_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}, Ref{SpiceBoolean}, Ref{SpiceBoolean}),
          selidx - 1, row - 1, elment - 1, ddata, null, found)
    handleerror()
    Bool(found[]) || return nothing
    Bool(null[]) && return missing
    ddata[]
end

"""
    ekgi(selidx, row, element)

Return an element of an entry in a column of integer type in a specified row.

### Arguments ###

- `selidx`: Index of parent column in SELECT clause
- `row`: Row to fetch from
- `elment`: Index of element, within column entry, to fetch

### Output ###

Returns the integer element of column entry or `missing` if it was null
or `nothing` if the column was not found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekgi_c.html)
"""
function ekgi(selidx, row, elment)
    null = Ref{SpiceBoolean}()
    found = Ref{SpiceBoolean}()
    idata = Ref{SpiceInt}()
    ccall((:ekgi_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceBoolean}, Ref{SpiceBoolean}),
          selidx - 1, row - 1, elment - 1, idata, null, found)
    handleerror()
    Bool(found[]) || return nothing
    Bool(null[]) && return missing
    Int(idata[])
end

"""
    ekinsr(handle, segno, recno)

Add a new, empty record to a specified E-kernel segment at a specified index.

### Arguments ###

- `handle`: File handle
- `segno`: Segment number
- `recno`: Record number

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekinsr_c.html)
"""
function ekinsr(handle, segno, recno)
    ccall((:ekinsr_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt),
          handle, segno - 1, recno - 1)
    handleerror()
end

"""
    ekifld(handle, tabnam, nrows, cnames, decls)

Initialize a new E-kernel segment to allow fast writing.

### Arguments ###

- `handle`: File handle
- `tabnam`: Table name
- `nrows`: Number of rows in the segment
- `cnames`: Names of columns
- `decls`: Declarations of columns

### Output ###

- `segno`: Segment number
- `rcptrs`: Array of record pointers

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekifld_c.html)
"""
function ekifld(handle, tabnam, nrows, cnames, decls)
    cnames_, ncols, cnmlen = chararray(cnames)
    @checkdims ncols decls
    decls_, _, declen = chararray(decls)
    segno = Ref{SpiceInt}()
    rcptrs = Array{SpiceInt}(undef, nrows)
    ccall((:ekifld_c, libcspice), Cvoid,
          (SpiceInt, Cstring, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceChar}, SpiceInt,
           Ref{SpiceChar}, Ref{SpiceInt}, Ref{SpiceInt}),
          handle, tabnam, ncols, nrows, cnmlen, cnames_, declen, decls_, segno, rcptrs)
    handleerror()
    segno[] + 1, rcptrs
end

"""
    eklef(fname)

Load an EK file, making it accessible to the EK readers.

### Arguments ###

- `fname`: Name of EK file to load

### Output ###

Returns the file handle of loaded EK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eklef_c.html)
"""
function eklef(fname)
    handle = Ref{SpiceInt}()
    ccall((:eklef_c, libcspice), Cvoid,
          (Cstring, Ref{SpiceInt}),
          fname, handle)
    handleerror()
    handle[]
end

"""
    eknelt(selidx, row)

Return the number of elements in a specified column entry in the current row.

### Arguments ###

- `selidx`: Index of parent column in SELECT clause
- `row`: Row containing element

### Output ###

Returns the number of elements in entry in current row.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eknelt_c.html)
"""
function eknelt(selidx, row)
    res = ccall((:eknelt_c, libcspice), SpiceInt,
                (SpiceInt, SpiceInt),
                selidx - 1, row - 1)
    handleerror()
    Int(res)
end

"""
    eknseg(handle)

Return the number of segments in a specified EK.

### Arguments ###

- `handle`: EK file handle

### Output ###

Returns the number of segments in the specified E-kernel.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eknseg_c.html)
"""
function eknseg(handle)
    res = ccall((:eknseg_c, libcspice), SpiceInt, (SpiceInt,), handle)
    handleerror()
    Int(res)
end

"""
    ekntab()

Return the number of loaded EK tables.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekntab_c.html)
"""
function ekntab()
    n = Ref{SpiceInt}()
    ccall((:ekntab_c, libcspice), Cvoid, (Ref{SpiceInt},), n)
    Int(n[])
end

"""
    ekopn(fname, ifname, ncomch)

Open a new E-kernel file and prepare the file for writing.

### Arguments ###

- `fname`: Name of EK file
- `ifname`: Internal file name
- `ncomch`: The number of characters to reserve for comments

### Output ###

Return the handle attached to the new EK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekopn_c.html)
"""
function ekopn(fname, ifname, ncomch)
    handle = Ref{SpiceInt}()
    ccall((:ekopn_c, libcspice), Cvoid,
          (Cstring, Cstring, SpiceInt, Ref{SpiceInt}),
          fname, ifname, ncomch, handle)
    handleerror()
    handle[]
end

"""
    ekopr(fname)

Open an existing E-kernel file for reading.

### Arguments ###

- `fname`: Name of EK file

### Output ###

Returns the handle attached to the EK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekopr_c.html)
"""
function ekopr(fname)
    handle = Ref{SpiceInt}()
    ccall((:ekopr_c, libcspice), Cvoid,
          (Cstring, Ref{SpiceInt}),
          fname, handle)
    handleerror()
    handle[]
end

"""
    ekops()

Open a scratch (temporary) E-kernel file and prepare the file for writing.

### Output ###

Returns the handle attached to the EK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekops_c.html)
"""
function ekops()
    handle = Ref{SpiceInt}()
    ccall((:ekops_c, libcspice), Cvoid, (Ref{SpiceInt},), handle)
    handleerror()
    handle[]
end

"""
    ekopw(fname)

Open an existing E-kernel file for writing.

### Arguments ###

- `fname`: Name of EK file

### Output ###

Returns the handle attached to the EK file.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekopw_c.html)
"""
function ekopw(fname)
    handle = Ref{SpiceInt}()
    ccall((:ekopw_c, libcspice), Cvoid,
          (Cstring, Ref{SpiceInt}),
          fname, handle)
    handleerror()
    handle[]
end

"""
    ekpsel(query, msglen=256, tablen=256, collen=256)

Parse the SELECT clause of an EK query, returning full particulars concerning each selected item.

### Arguments ###

- `query`: EK query
- `msglen`: Available space in the output error message string (default: 256)
- `tablen`: Length of strings in `tabs' output array (default: 256)
- `collen`: Length of strings in `cols' output array (default: 256)

### Output ###

- `xbegs`: Begin positions of expressions in SELECT clause
- `xends`: End positions of expressions in SELECT clause
- `xtypes`: Data types of expressions
- `xclass`: Classes of expressions
- `tabs`: Names of tables qualifying SELECT columns
- `cols`: Names of columns in SELECT clause of query

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekpsel_c.html)
"""
function ekpsel(query, msglen=256, tablen=256, collen=256)
    maxqsel = 100
    n = Ref{SpiceInt}()
    xbegs = Array{SpiceInt}(undef, maxqsel)
    xends = Array{SpiceInt}(undef, maxqsel)
    xtypes = Array{SpiceInt}(undef, maxqsel)
    xclass = Array{SpiceInt}(undef, maxqsel)
    tabs = Array{SpiceChar}(undef, tablen, maxqsel)
    cols = Array{SpiceChar}(undef, collen, maxqsel)
    error = Ref{SpiceBoolean}()
    errmsg = Array{SpiceChar}(undef, msglen)
    ccall((:ekpsel_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceInt},
           Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceChar}, Ref{SpiceBoolean},
           Ref{SpiceChar}),
          query, msglen, tablen, collen, n, xbegs, xends, xtypes, xclass, tabs, cols, error, errmsg)
    handleerror()
    if Bool(error[])
        throw(SpiceError(chararray_to_string(errmsg)))
    end
    Int.(xbegs[1:n[]]), Int.(xends[1:n[]]), Int.(xtypes[1:n[]]), Int.(xclass[1:n[]]),
    chararray_to_string(tabs)[1:n[]], chararray_to_string(cols)[1:n[]]
end

"""
    ekrcec(handle, segno, recno, column, lenout=256, nelts=100)

Read data from a character column in a specified EK record.

### Arguments ###

- `handle`: Handle attached to EK file
- `segno`: Index of segment containing record
- `recno`: Record from which data is to be read
- `column`: Column name
- `lenout`: Maximum length of output strings
- `nelts`: Maximum number of elements to return (default: 100)

### Output ###

Returns the character values in column entry or `missing` if they are null.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekrcec_c.html)
"""
function ekrcec(handle, segno, recno, column, lenout=256, nelts=100)
    nvals = Ref{SpiceInt}()
    cvals = Array{SpiceChar}(undef, lenout, nelts)
    isnull = Ref{SpiceBoolean}()
    ccall((:ekrcec_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt,
           Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceBoolean}),
          handle, segno - 1, recno - 1, column, lenout, nvals, cvals, isnull)
    handleerror()
    Bool(isnull[]) && return missing
    strip.(chararray_to_string(cvals)[1:nvals[]])
end

"""
    ekrced(handle, segno, recno, column, nelts=100)

Read data from a double precision column in a specified EK record.

### Arguments ###

- `handle`: Handle attached to EK file
- `segno`: Index of segment containing record
- `recno`: Record from which data is to be read
- `column`: Column name
- `nelts`: Maximum number of elements to return (default: 100)

### Output ###

Returns the values in column entry.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekrced_c.html)
"""
function ekrced(handle, segno, recno, column, nelts=100)
    nvals = Ref{SpiceInt}()
    dvals = Array{SpiceDouble}(undef, nelts)
    isnull = Ref{SpiceBoolean}()
    ccall((:ekrced_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring,
           Ref{SpiceInt}, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          handle, segno - 1, recno - 1, column, nvals, dvals, isnull)
    handleerror()
    Bool(isnull[]) && return missing
    dvals[1:nvals[]]
end

"""
    ekrcei(handle, segno, recno, column, nelts=100)

Read data from an integer column in a specified EK record.

### Arguments ###

- `handle`: Handle attached to EK file
- `segno`: Index of segment containing record
- `recno`: Record from which data is to be read
- `column`: Column name
- `nelts`: Maximum number of elements to return (default: 100)

### Output ###

Returns the values in column entry.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekrcei_c.html)
"""
function ekrcei(handle, segno, recno, column, nelts=100)
    nvals = Ref{SpiceInt}()
    ivals = Array{SpiceInt}(undef, nelts)
    isnull = Ref{SpiceBoolean}()
    ccall((:ekrcei_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring,
           Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceBoolean}),
          handle, segno - 1, recno - 1, column, nvals, ivals, isnull)
    handleerror()
    Bool(isnull[]) && return missing
    Int.(ivals[1:nvals[]])
end

"""
    ekssum(handle, segno)

Return summary information for a specified segment in a specified EK.

### Arguments ###

- `handle`: Handle of EK
- `segno`: Number of segment to be summarized

### Output ###

Returns the EK segment summary.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekssum_c.html)
"""
function ekssum(handle, segno)
    sum = Ref{EKSegSum}()
    ccall((:ekssum_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{EKSegSum}),
          handle, segno - 1, sum)
    handleerror()
    sum[]
end

"""
    ektnam(n, lenout=256)

Return the name of a specified, loaded table.

### Arguments ###

- `n`: Index of table
- `lenout`: Maximum table name length (default: 256)

### Output ###

Returns the name of table.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ektnam_c.html)
"""
function ektnam(n, lenout=256)
    table = Array{SpiceChar}(undef, lenout)
    ccall((:ektnam_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{SpiceChar}),
          n - 1, lenout, table)
    handleerror()
    chararray_to_string(table)
end

"""
    ekucec(handle, segno, recno, column, cvals, isnull)

Update a character column entry in a specified EK record.

### Arguments ###

- `handle`: EK file handle
- `segno`: Index of segment containing record
- `recno`: Record to which data is to be updated
- `column`: Column name
- `cvals`: Character values comprising new column entry
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekucec_c.html)
"""
function ekucec(handle, segno, recno, column, cvals, isnull)
    cvals_, nvals, vallen = chararray(cvals)
    ccall((:ekucec_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt, SpiceInt,
           Ref{SpiceChar}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, vallen, cvals_, isnull)
    handleerror()
end

"""
    ekuced(handle, segno, recno, column, dvals, isnull)

Update a double precision column entry in a specified EK record.

### Arguments ###

- `handle`: Handle attached to EK file
- `segno`: Index of segment containing record
- `recno`: Record in which entry is to be updated
- `column`: Column name
- `dvals`: Double precision values comprising new column entry
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekuced_c.html)
"""
function ekuced(handle, segno, recno, column, dvals, isnull)
    nvals = length(dvals)
    ccall((:ekuced_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt,
           Ref{SpiceDouble}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, dvals, isnull)
    handleerror()
end

"""
    ekucei(handle, segno, recno, column, dvals, isnull)

Update an integer column entry in a specified EK record.

### Arguments ###

- `handle`: Handle attached to EK file
- `segno`: Index of segment containing record
- `recno`: Record in which entry is to be updated
- `column`: Column name
- `ivals`: Integer values comprising new column entry
- `isnull`: Flag indicating whether column entry is null

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekucei_c.html)
"""
function ekucei(handle, segno, recno, column, ivals, isnull)
    nvals = length(ivals)
    ivals_ = SpiceInt.(ivals)
    ccall((:ekucei_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, Cstring, SpiceInt,
           Ref{SpiceInt}, SpiceBoolean),
          handle, segno - 1, recno - 1, column, nvals, ivals_, isnull)
    handleerror()
end

"""
    ekuef(handle)

Unload an EK file, making its contents inaccessible to the EK reader routines, and clearing space
in order to allow other EK files to be loaded.

### Arguments ###

- `handle`: Handle of EK file

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/ekuef_c.html)
"""
function ekuef(handle)
    ccall((:ekuef_c, libcspice), Cvoid, (SpiceInt,), handle)
    handleerror()
end

"""

Convert an ellipse to a center vector and two generating vectors.
The selected generating vectors are semi-axes of the ellipse.

### Arguments ###

- `ellipse`: An ellipse

### Output ###

Returns the center and semi-axes of ellipse.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/el2cgv_c.html)
"""
function el2cgv(ellipse)
    center = Array{SpiceDouble}(undef, 3)
    smajor = Array{SpiceDouble}(undef, 3)
    sminor = Array{SpiceDouble}(undef, 3)
    ccall((:el2cgv_c, libcspice), Cvoid,
          (Ref{Ellipse}, Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          ellipse, center, smajor, sminor)
    center, smajor, sminor
end

@deprecate elemc(item, cell) item in cell
@deprecate elemd(item, cell) item in cell
@deprecate elemi(item, cell) item in cell

"""
    elem[c/d/i](item, cell)

!!! warning "Deprecated"
    Use `item in cell` instead.
"""
elemc
elemd
elemi

"""
    eqstr(a, b)

Determine whether two strings are equivalent.

### Arguments ###

- `a`, `b`: Arbitrary character strings

### Output ###

Returns `true` if `a` and `b` are equivalent.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eqstr_c.html)
"""
function eqstr(a, b)
    res = ccall((:eqstr_c, libcspice), SpiceBoolean,
                (Cstring, Cstring), a, b)
    handleerror()
    Bool(res)
end

"""
    eqncpv(et, epoch, eqel, rapol, decpol)

Compute the state (position and velocity of an object whose trajectory is described via equinoctial
elements relative to some fixed plane (usually the equatorial plane of some planet).

### Arguments ###

- `et`: Epoch in seconds past J2000 to find state
- `epoch`: Epoch of elements in seconds past J2000
- `eqel`: Array of equinoctial elements
- `rapol`: Right Ascension of the pole of the reference plane
- `decpol`: Declination of the pole of the reference plane

### Output ###

Returns the state of the object described by `eqel`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eqncpv_c.html)
"""
function eqncpv(et, epoch, eqel, rapol, decpol)
    @checkdims 9 eqel
    state = Array{SpiceDouble}(undef, 6)
    ccall((:eqncpv_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, Ref{SpiceDouble}, SpiceDouble, SpiceDouble, Ref{SpiceDouble}),
          et, epoch, eqel, rapol, decpol, state)
    handleerror()
    state
end

"""
    esrchc(value, array)

Search for a given value within a character string array.

### Arguments ###

- `value`: Key value to be found in array
- `array`: Character string array to search

### Output ###

Returns the index of the first equivalent array entry, or -1 if no equivalent element is found.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/esrchc_c.html)
"""
function esrchc(value, array)
    array_, ndim, lenvals = chararray(array)
    res = Int(ccall((:esrchc_c, libcspice), SpiceInt,
                    (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
                    value, ndim, lenvals, array_))
    handleerror()
    res == -1 ? res : res + 1
end

"""
    et2lst(et, body, lon, typ, timlen=128, ampmlen=128)

Given an ephemeris epoch, compute the local solar time for an object on the surface of a body at a
specified longitude.

### Arguments ###

- `et`: Epoch in seconds past J2000 epoch
- `body`: ID-code of the body of interest
- `lon`: Longitude of surface point (radians)
- `typ`: Type of longitude "PLANETOCENTRIC", etc
- `timlen`: Available room in output time string (default: 128)
- `ampmlen`: Available room in output `ampm' string (default: 128)

### Output ###

- `hr`: Local hour on a "24 hour" clock
- `mn`: Minutes past the hour
- `sc`: Seconds past the minute
- `time`: String giving local time on 24 hour clock
- `ampm`: String giving time on A.M./ P.M. scale

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/et2lst_c.html)
"""
function et2lst(et, body, lon, typ, timlen=128, ampmlen=128)
    hr = Ref{SpiceInt}()
    mn = Ref{SpiceInt}()
    sc = Ref{SpiceInt}()
    time = Array{SpiceChar}(undef, timlen)
    ampm = Array{SpiceChar}(undef, ampmlen)
    ccall((:et2lst_c, libcspice), Cvoid,
          (SpiceDouble, SpiceInt, SpiceDouble, Cstring, SpiceInt, SpiceInt,
           Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceChar}),
          et, body, lon, typ, timlen, ampmlen, hr, mn, sc, time, ampm)
    handleerror()
    Int(hr[]), Int(mn[]), Int(sc[]), chararray_to_string(time), chararray_to_string(ampm)
end

"""
    et2utc(et, format, prec)

Convert an input time from ephemeris seconds past J2000
to Calendar, Day-of-Year, or Julian Date format, UTC.

### Arguments ###

- `et`: Input epoch, given in ephemeris seconds past J2000
- `format`: Format of output epoch. It may be any of the following:
    - `:C`: Calendar format, UTC
    - `:D`: Day-of-Year format, UTC
    - `:J`: Julian Date format, UTC
    - `:ISOC`: ISO Calendar format, UTC
    - `:ISOD`: ISO Day-of-Year format, UTC
- `prec`: Digits of precision in fractional seconds or days

### Output ###

Returns an output time string equivalent to the input epoch, in the specified format.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/et2utc_c.html)
"""
function et2utc(et, format, prec)
    lenout = 32
    utcstr = Array{SpiceChar}(undef, lenout)
    ccall((:et2utc_c, libcspice), Cvoid,
          (SpiceDouble, Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}),
          et, string(format), prec, lenout, utcstr)
    handleerror()
    chararray_to_string(utcstr)
end

"""
    etcal(et, lenout=128)

Convert from an ephemeris epoch measured in seconds past the epoch of J2000 to
a calendar string format using a formal calendar free of leapseconds.

### Arguments ###

- `et`: Ephemeris time measured in seconds past J2000
- `lenout`: Length of output string (default: 128)

### Output ###

Returns a standard calendar representation of `et`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/etcal_c.html)
"""
function etcal(et, lenout=128)
    string = Array{SpiceChar}(undef, lenout)
    ccall((:etcal_c, libcspice), Cvoid, (SpiceDouble, SpiceInt, Ref{SpiceChar}),
          et, lenout, string)
    chararray_to_string(string)
end

"""
    eul2m(angle3, angle2, angle1, axis3, axis2, axis1)

Construct a rotation matrix from a set of Euler angles.

### Arguments ###

- `angle3`, `angle2`, `angle1`: Rotation angles about third, second, and first rotation axes (radians)
- `axis3`, `axis2`, `axis1`: Axis numbers of third, second, and first rotation axes

### Output ###

A rotation matrix corresponding to the product of the 3 rotations.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eul2m_c.html)
"""
function eul2m(angle3, angle2, angle1, axis3, axis2, axis1)
    r = Matrix{SpiceDouble}(undef, 3, 3)
    ccall((:eul2m_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          angle3, angle2, angle1, axis3, axis2, axis1, r)
    handleerror()
    permutedims(r)
end

"""
    eul2xf(eulang, axisa, axisb, axisc)

Compute a state transformation from an Euler angle factorization of a rotation and the derivatives
of those Euler angles.

### Arguments ###

- `eulang`: An array of Euler angles and their derivatives
- `axisa`: Axis A of the Euler angle factorization
- `axisb`: Axis B of the Euler angle factorization
- `axisc`: Axis C of the Euler angle factorization

### Output ###

Returns a state transformation matrix.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/eul2xf_c.html)
"""
function eul2xf(eulang, axisa, axisb, axisc)
    @checkdims 6 eulang
    xform = Array{SpiceDouble}(undef, 6, 6)
    ccall((:eul2xf_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}),
          eulang, axisa, axisb, axisc, xform)
    handleerror()
    permutedims(xform)
end

"""
    expool(name)

Confirm the existence of a kernel variable in the kernel pool.

### Arguments ###

- `name`: Name of the variable whose value is to be returned

### Output ###

Returns `true` when the variable is in the pool.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/expool_c.html)
"""
function expool(name)
    found = Ref{SpiceBoolean}()
    ccall((:expool_c, libcspice), Cvoid,
          (Cstring, Ref{SpiceBoolean}),
          name, found)
    handleerror()
    Bool(found[])
end


text:
mutable struct Plan{T}
    n_districts::T
    assignment::Array{T}
    district_populations::Array{T}
    cut_edges::BitSet
    n_cut_edges::T
end

function Plan(graph::GraphIndex, assignment::Array{Int64})::Plan{Int64}
    n_districts = maximum(assignment)
    district_populations = zeros(Int64, n_districts)
    for (index, node_pop) in enumerate(graph.population)
        district_populations[assignment[index]] += node_pop
    end
    cut_edges = BitSet(Int64[])
    for index in 1:graph.n_edges
        left_assignment = assignment[graph.edges[1, index]]
        right_assignment = assignment[graph.edges[2, index]]
        if left_assignment != right_assignment
            push!(cut_edges, index)
        end
    end
    return Plan{Int64}(n_districts, assignment, district_populations,
                       cut_edges, length(cut_edges))
end

max_stars_repo_name: pjrule/gridwalk

max_issues_repo_name: pjrule/gridwalk

max_forks_repo_name: pjrule/gridwalk