text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--epochs"
        help = "Epochs to train"
        arg_type = Int
        default = 100
    "--save-snapshot"
        help = "Whether to Save Snapshots"
        arg_type = Bool
        default = true
    "--snapshot-dir"
        help = "Directory to Store Snapshots"
        arg_type = ASCIIString
        default = "snapshot-vae"
    "--snapshot_n_iter"
        help = "How many iterations between snapshots"
        arg_type = Int
        default = 10
    "--stat_n_iter"
        help = "How many iterations between saving statistics"
        arg_type = Int
        default = 10
end
const args = parse_args(ARGS, s)

ENV["SEEP_NO_GPU"]=false
using MNIST, Seep, JLD
T = Float64
const pool = BuddyPool(Array{T}(1<<24))
const ftrain = map(T, traindata()[1]/255)
import Seep: @named

# Constants
const input_size = 28*28
const latent_size = 5
const output_size = input_size
const batch_size = 100

# Network
X = ANode("X", zeros(T, input_size, batch_size))
μ = ANode("μ", zeros(T, latent_size, batch_size))
lnσ = ANode("lnσ", zeros(T, latent_size, batch_size))

x = X

# Generate latent variables (encoding)
encode_params = Flow[]
let encode_size=[20, 20]

    fe = Linear(T, "e_encode", input_size, encode_size[1])
    x = tanh(fe(x))

    fh = Flow[]
    for i in 1:length(encode_size)
        j = i==1 ? i : i-1
        fi = Linear(T, "h_$(i)_encode", encode_size[j], encode_size[i])
        x = tanh(fi(x))
        push!(fh, fi)
    end
    fs = Linear(T, "s_encode", encode_size[end], latent_size)
    x = tanh(fs(x))

    fμ = Linear(T, "μ_encode", latent_size, latent_size)
    μ = fμ(x)

    fσ = Linear(T, "lnσ_encode", latent_size, latent_size)
    lnσ = fσ(x)

    push!(encode_params, fe, fh..., fs, fμ, fσ)
end

# Sample Latent Variables
@named ϵ = randn!(ANode(zeros(T, latent_size, batch_size)))
@named z = μ + exp(lnσ) .* ϵ

# Inference of image from latent variables (decode)
yhat = z
decode_params = Flow[]
let decode_size=[20, 20]

    fe = Linear(T, "e_decode", latent_size, decode_size[1])
    yhat = tanh(fe(yhat))

    fh = Flow[]
    for i in 1:length(decode_size)
        j = i==1 ? i : i-1
        fi = Linear(T, "h_$(i)_decode", decode_size[j], decode_size[i])
        yhat = tanh(fi(yhat))
        push!(fh, fi)
    end
    fs = Linear(T, "s_decode", decode_size[end], output_size)
    yhat = tanh(fs(yhat))

    fy = Linear(T, "y_decode", output_size, output_size)
    yhat = fy(yhat)

    push!(decode_params, fe, fh..., fs, fy)
end

# Loss
@named begin
  kl   = (0.5 / latent_size) * sum(μ.^2 + exp(lnσ) - lnσ - 1, 1)
  loss = 0.5*sum((X - yhat).^2, 1)
  obj  = sum(loss) + sum(kl)
end

grad = gradients(obj)

# Gradient Calculation for SGD
updates = map(x->Seep.adam(x, grad[x]), get_params(Flow[encode_params; decode_params]))
const VAE = instance(pool, [X; yhat; kl; loss; updates])
Seep.print_stats(STDOUT, pool)

# Load Latest Snapshot
n_epoch_start = 0

if args["save-snapshot"]
    if isdir(args["snapshot-dir"])
        f = filter(x->contains(x, "snapshot-"), readdir(args["snapshot-dir"]))
        if length(f) > 0
            epochs = map(x->parse(Int, split(split(x, ".")[1], "-")[end]), f)
            ind = indmax(epochs)
            n_epoch_start = epochs[ind]
            @show n_epoch_start
            load_snapshot(joinpath(args["snapshot-dir"], f[ind]), loss)
        end
    else
        mkdir(args["snapshot-dir"])
    end
end

# Train
function epoch()
    perm = randperm(60000)
    for i in 1:batch_size:60000
        ii = perm[i+(1:batch_size)-1]

        # Input Image
        VAE.X[:,:] = ftrain[:, ii]

        VAE()
    end
end

Erec = zeros(args["stat_n_iter"])
Ekl = zeros(args["stat_n_iter"])
for i = 1:args["epochs"]
    epoch()

    rec_loss = sqrt(sum(VAE.loss))/batch_size
    kl_loss = sum(VAE.kl)/batch_size
    println("Rec Loss: ", rec_loss, " KL Loss:", kl_loss)

    Erec[mod(i-1, args["stat_n_iter"])+1] = rec_loss
    Ekl[mod(i-1, args["stat_n_iter"])+1] = kl_loss
    if mod(i, args["stat_n_iter"])==0
        f = joinpath(args["snapshot-dir"], "statistics.jld")
        if isfile(f)
            Erec = vcat(load(f, "Erec"), Erec)
            Ekl = vcat(load(f, "Ekl"), Ekl)
        end
        save(f, "Erec", Erec, "Ekl", Ekl)
        Erec = zeros(args["stat_n_iter"])
        Ekl = zeros(args["stat_n_iter"])
    end

    if mod(i, args["snapshot_n_iter"])==0 && args["save-snapshot"]
        save_snapshot(joinpath(args["snapshot-dir"], "snapshot-$(i+n_epoch_start).jld"), loss)
    end
end


text:
# we generate code in this module, so precompile where possible
#__precompile__(true)

module Rotations

using StaticArrays

import Base: convert, eltype, size, length, getindex, inv, *, Tuple, eye

include("core_types.jl")
include("quaternion_types.jl")
include("angleaxis_types.jl")
include("euler_types.jl")
include("mean.jl")
include("derivatives.jl")

export
    Rotation, RotMatrix,
    Quat, SPQuat,
    AngleAxis, RodriguesVec,
    RotX, RotY, RotZ,
    RotXY, RotYX, RotZX, RotXZ, RotYZ, RotZY,
    RotXYX, RotYXY, RotZXZ, RotXZX, RotYZY, RotZYZ,
    RotXYZ, RotYXZ, RotZXY, RotXZY, RotYZX, RotZYX,

    # check validity of the rotation (is it close to unitary?)
    isrotation,

    # angle and axis introspection
    rotation_angle,
    rotation_axis

    # derivatives (names clash with ForwarDiff?)
    #jacobian, hessian

end # module

ext: jl


text:
using SuperResolution
using Test

@testset "SuperResolution" begin
end


text:
"""
   `loadDataset(filepath)`

Reads a dataset from a csv file located at `filepath`.
Each row of the dataset should be an
example, each column a feature. The last column are the targets, which
are assumed to be binary labels 0 and 1.
A column of ones is added to the beginning of the dataset and the set is
equally split into training and test set.

### Returns
 * train_x: the data for training
 * train_y: the labels of the training data
 * test_x: the data for testing
 * test_y: the labels of the test data
"""
function loadDataset(filepath)
    ## load dataset
    data = readdlm(filepath, ',')

    ## preprocess dataset
    print("preprocessing dataset...\t")
    tic()
    # split into features and labels:
    x = data[:,1:end-1]
    y = data[:,end]
    # replace 0 entries with -1
    y[y .== 0] = -1
    # add a column of ones to the beginning of x
    # because we encorporated b into w
    x0 = ones(size(x,1))
    x = hcat(x0,x)

    # seed for random number generator:
    srand(1234)
    # shuffle the training data:
    perm = randperm(size(x,1))
    x = x[perm,:]
    y = y[perm]

    # use every second datapoint for training
    # TODO: use more data for training
    train_x = x[1:2:end,:]
    train_y = y[1:2:end]
    test_x = x[2:2:end,:]
    test_y = y[2:2:end]

    @printf "finished after %.3f seconds\n" toq()

    return (train_x,train_y,test_x,test_y)
end

"""
   `testResult(w, test_x, test_y)`

Measures the classification performance of a hyperplane with normal
vector `w` on test data.

### Returns
 * `acc`: accuracy on the test data set
"""
function testResult(w, test_x, test_y)
    predicted_y = sign(test_x*w)
    acc = mean((predicted_y .== test_y))
    @printf "Accuracy: %.2f%%\n" acc*100
    return acc
end


text:
module QD

using Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("QD not installed properly, run Pkg.build(\"QD\"), restart Julia, and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

# TODO: Check whether this is actually true
const fma_is_fast = true

# Basic (internal) Functions

# Computes fl(a+b) and err(a+b).  Assumes |a| >= |b|.
function quick_two_sum(a::T, b::T) where {T <: AbstractFloat}
    s = a + b
    err = b - (s - a)
    s, err
end

# Computes fl(a-b) and err(a-b).  Assumes |a| >= |b|.
function quick_two_diff(a::T, b::T) where {T <: AbstractFloat}
    s = a - b
    err = (a - s) - b
    s, err
end

# Computes fl(a+b) and err(a+b).
function two_sum(a::T, b::T) where {T <: AbstractFloat}
    s = a + b
    bb = s - a
    err = (a - (s - bb)) + (b - bb)
    s, err
end

# Computes fl(a-b) and err(a-b).
function two_diff(a::T, b::T) where {T <: AbstractFloat}
    s = a - b
    bb = s - a
    err = (a - (s - bb)) - (b + bb)
    s, err
end

# Computes high word and lo word of a
function split(a::T) where {T <: AbstractFloat}
    _QD_SPLITTER = 134217729.0               # = 2^27 + 1
    _QD_SPLIT_THRESH = 6.69692879491417e+299 # = 2^996
    if a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH
        a *= 3.7252902984619140625e-09 # 2^-28
        temp = _QD_SPLITTER * a
        hi = temp - (temp - a)
        lo = a - hi
        hi *= 268435456.0       # 2^28
        lo *= 268435456.0       # 2^28
    else
        temp = _QD_SPLITTER * a
        hi = temp - (temp - a)
        lo = a - hi
    end
    hi, lo
end

# Computes fl(a*b) and err(a*b).
function two_prod(a::T, b::T) where {T <: AbstractFloat}
    if fma_is_fast
        p = a * b
        err = fma(a, b, -p)
    else
        p = a * b
        a_hi, a_lo = split(a)
        b_hi, b_lo = split(b)
        err = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo
    end
    p, err
end

# Computes fl(a*a) and err(a*a).  Faster than the above method.
function two_sqr(a::T) where {T <: AbstractFloat}
    if fma_is_fast
        p = a * a
        err = fma(a, a, - p)
        p, err
    else
        q = a * a
        hi, lo = split(a)
        err = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo
        q, err
    end
end

# Additions

function three_sum(a::T, b::T, c::T) where {T <: AbstractFloat}
    t1, t2 = two_sum(a, b)
    a , t3 = two_sum(c, t1)
    b , c  = two_sum(t2, t3)
    a, b, c
end

function three_sum2(a::T, b::T, c::T) where {T <: AbstractFloat}
    t1, t2 = two_sum(a, b)
    a , t3 = two_sum(c, t1)
    b = t2 + t3
    a, b, c
end

# Renormalization
function quick_renorm(c0::T, c1::T, c2::T, c3::T, c4::T) where
        {T <: AbstractFloat}
    s , t3 = quick_two_sum(c3, c4)
    s , t2 = quick_two_sum(c2, s )
    s , t1 = quick_two_sum(c1, s )
    c0, t0 = quick_two_sum(c0, s )

    s , t2 = quick_two_sum(t2, t3)
    s , t1 = quick_two_sum(t1, s )
    c1, t0 = quick_two_sum(t0, s )

    s , t1 = quick_two_sum(t1, t2)
    c2, t0 = quick_two_sum(t0, s )

    c3 = t0 + t1

    c0, c1, c2, c3, c4
end

function renorm(c0::T, c1::T, c2::T, c3::T) where {T <: AbstractFloat}
    s2 = 0.0
    s3 = 0.0

    isinf(c0) && return c0, c1, c2, c3

    s0, c3 = quick_two_sum(c2, c3)
    s0, c2 = quick_two_sum(c1, s0)
    c0, c1 = quick_two_sum(c0, s0)

    s0 = c0
    s1 = c1
    if s1 != 0.0
        s1, s2 = quick_two_sum(s1, c2)
        if s2 != 0.0
            s2, s3 = quick_two_sum(s2, c3)
        else
            s1, s2 = quick_two_sum(s1, c3)
        end
    else
        s0, s1 = quick_two_sum(s0, c2)
        if s1 != 0.0
            s1, s2 = quick_two_sum(s1, c3)
        else
            s0, s1 = quick_two_sum(s0, c3)
        end
    end

    c0 = s0
    c1 = s1
    c2 = s2
    c3 = s3

    c0, c1, c2, c3
end

function renorm(c0::T, c1::T, c2::T, c3::T, c4::T) where {T <: AbstractFloat}
    s2 = 0.0
    s3 = 0.0

    isinf(c0) && return c0, c1, c2, c3, c4

    s0, c4 = quick_two_sum(c3, c4)
    s0, c3 = quick_two_sum(c2, s0)
    s0, c2 = quick_two_sum(c1, s0)
    c0, c1 = quick_two_sum(c0, s0)

    s0 = c0
    s1 = c1

    if s1 != 0.0
        s1, s2 = quick_two_sum(s1, c2)
        if s2 != 0.0
            s2, s3 = quick_two_sum(s2, c3)
            if s3 != 0.0
                s3 += c4
            else
                s2, s3 = quick_two_sum(s2, c4)
            end
        else
            s1, s2 = quick_two_sum(s1, c3)
            if s2 != 0.0
                s2, s3 = quick_two_sum(s2, c4)
            else
                s1, s2 = quick_two_sum(s1, c4)
            end
        end
    else
        s0, s1 = quick_two_sum(s0, c2)
        if s1 != 0.0
            s1, s2 = quick_two_sum(s1, c3)
            if s2 != 0.0
                s2, s3 = quick_two_sum(s2, c4)
            else
                s1, s2 = quick_two_sum(s1, c4)
            end
        else
            s0, s1 = quick_two_sum(s0, c3)
            if s1 != 0.0
                s1, s2 = quick_two_sum(s1, c4)
            else
                s0, s1 = quick_two_sum(s0, c4)
            end
        end
    end

    c0 = s0
    c1 = s1
    c2 = s2
    c3 = s3

    c0, c1, c2, c3, c4
end

include("QD128.jl")
include("QD256.jl")

end


text:
"""Symbolic Omega"""
struct SymΩ{I} <: Space.ΩBase{I}
  # vals::Dict{I, V}
end

struct SymVal{K, V}
  val::NamedTuple{K, V}
end

Base.:*(x::SymVal, y::SymVal) =
  isindep(x, y) ? SymVal((mean = x.val.mean * y.val.mean, ig = x.val.ig ∪ y.val.ig)) : error("Cannot multiply expectatios of non-independent randvars")
Base.:+(x::SymVal, y::SymVal) = SymVal((mean = x.val.mean + y.val.mean, ig = x.val.ig ∪ y.val.ig))
Base.:+(x::SymVal, y) = SymVal((mean = x.val.mean + y, ig = x.val.ig))
Base.:+(x, y::SymVal) = SymVal((mean = x + y.val.mean, ig = y.val.ig))
Base.:-(x::SymVal, y) = SymVal((mean = x.val.mean - y, ig = x.val.ig))
Base.:-(x, y::SymVal) = SymVal((mean = x - y.val.mean, ig = y.val.ig))

S(::Type{T}) where T = Union{SymVal, T}

const SymΩProj = Space.ΩProj{SymΩ{ID}, ID} where ID

function Prim.normal(ωπ::SymΩProj, μ::Real, σ::Real)
  Space.randinc!(ωπ, SymVal((mean = μ, ig = Set([ωπ.id]))))
end

function Prim.normal(ωπ::SymΩProj, μ::SymVal, σ::SymVal)
  # @assert false
  Space.randinc!(ωπ, SymVal((mean = μ.val.mean, ig =  union(Set([ωπ.id]), μ.val.ig, σ.val.ig))))
end

function Prim.uniform(ωπ::SymΩProj, a::Real, b::Real)
  Space.randinc!(ωπ, SymVal((mean = (a + b)*0.5, ig = Set([ωπ.id]))))
end

"Independence Group"
struct IndepGroup
  ids::Set{ID}
end

IndepGroup(::ID) = IndepSet(id)

Base.union(x::IndepGroup, y::IndepGroup) = IndepGroup(union(x.ids, y.ids))
unionig(xs...) = union(xs...)
isintersect(x::Set, y::Set) = isempty(intersect(x, y))
isindep(x::Set, y::Set) = isintersect(x, y)
isindep(x::SymVal, y::SymVal) = isindep(x.val.ig, y.val.ig)

function analyticalmean(x)
  ω = SymΩ{ID}()
  x(ω).val.mean
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

# data input
(V, FV, EV) = ([1.01181 0.215639 0.91979 0.123616 1.02252 0.226347 0.930498 0.134324 0.0458309 -0.301827 0.348275 0.0006172 0.579367 0.23171 0.881811 0.534154 -0.0521776 0.627953 -0.190635 0.489496 -0.0233884 0.656742 -0.161846 0.518285 0.27652 -0.0875132 0.52527 0.161237 0.509324 0.145291 0.758074 0.394041 0.27631 0.564484 0.0058279 0.294002 1.01221 1.30039 0.741732 1.02991; 0.160033 0.0680099 0.956278 0.864255 0.160649 0.0686266 0.956895 0.864872 -0.200245 0.102199 0.417839 0.720283 -0.35354 -0.0510965 0.264543 0.566987 0.682359 0.543901 0.0592036 -0.0792537 0.956374 0.817917 0.333219 0.194762 -0.102028 0.146722 0.324834 0.573584 -0.16916 0.0795901 0.257702 0.506452 0.452356 0.181874 1.15396 0.883481 0.816151 0.545669 1.51776 1.24728; 0.196256 0.206963 0.196872 0.20758 0.997729 1.00844 0.998346 1.00905 0.0677451 0.601282 -0.0855504 0.447986 0.502301 1.03584 0.349005 0.882542 0.159301 0.18809 0.433316 0.462105 0.797002 0.825792 1.07102 1.09981 0.1446 0.377404 0.0774682 0.310272 0.580364 0.813168 0.513232 0.746036 0.403805 1.13971 0.767599 1.5035 0.249344 0.985248 0.613139 1.34904], Array{Int64,1}[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8],[9,10,11,12],[13,14,15,16],[9,10,13,14],[11,12,15,16],[9,11,13,15],[10,12,14,16],[17,18,19,20],[21,22,23,24],[17,18,21,22],[19,20,23,24],[17,19,21,23],[18,20,22,24],[25,26,27,28],[29,30,31,32],[25,26,29,30],[27,28,31,32],[25,27,29,31],[26,28,30,32],[33,34,35,36],[37,38,39,40],[33,34,37,38],[35,36,39,40],[33,35,37,39],[34,36,38,40]],Array{Int64,1}[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8],[9,10],[11,12],[13,14],[15,16],[9,11],[10,12],[13,15],[14,16],[9,13],[10,14],[11,15],[12,16],[17,18],[19,20],[21,22],[23,24],[17,19],[18,20],[21,23],[22,24],[17,21],[18,22],[19,23],[20,24],[25,26],[27,28],[29,30],[31,32],[25,27],[26,28],[29,31],[30,32],[25,29],[26,30],[27,31],[28,32],[33,34],[35,36],[37,38],[39,40],[33,35],[34,36],[37,39],[38,40],[33,37],[34,38],[35,39],[36,40]])

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.space_arrangement(
	W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)

cc = [copEV, copFE, copCF]
LarModelString = Lar.lar2obj(V::Lar.Points, cc::Lar.ChainComplex)
f = open("../test/out3d.obj", "w")
print(f, LarModelString); close(f)
V,EVs,FVs = Lar.obj2lar("../test/out3d.obj")

GL.VIEW([
      GL.GLPolyhedron(V, FVs[1])
      GL.GLFrame
]);

V,CVs,FVs,EVs = Lar.pols2tria(V, copEV, copFE, copCF); # whole assembly
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,CVs[2:end],5,5,5,99,0.2));
GL.VIEW(GL.GLExplode(V,[CVs[1]],1,1,1,99,0.5));


text:
# Autogenerated wrapper script for BWA_jll for x86_64-linux-gnu
export bwa, libbwa

using Zlib_jll
JLLWrappers.@generate_wrapper_header("BWA")
JLLWrappers.@declare_executable_product(bwa)
JLLWrappers.@declare_library_product(libbwa, "libbwa.so")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_executable_product(
        bwa,
        "bin/bwa",
    )

    JLLWrappers.@init_library_product(
        libbwa,
        "lib/libbwa.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BrainForth: @bf, @run

# --
@bf intro = ["The Guessing Game: Guess a number from 0-255!", println]

# -- name
@bf getname = [["Enter your name: ", prompt], whileempty]
# -- n
@bf getnum = [["Enter your guess: ", prompt], whileempty, parseint]
# --
@bf greeting = ["Hello, ", print, print, "!", println]

# Lacking an RNG, we come up with a number based on the user's name.
# s -- s n
@bf roll = [[1, +], map, dupv, prod]

# n --
@bf game = [[dup, getnum, tuck, !=],
            [dupd, <, ["Too high!"], ["Too low!"], iff, println],
            loop, drop, drop,
            "Correct!", println]

@bf main = [intro, getname, dupv, greeting,
            [1], [roll, game, "", println], loop]

@run [main]

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
# ---
# title: 310. Minimum Height Trees
# id: problem310
# author: Indigo
# date: 2021-03-20
# difficulty: Medium
# categories: Breadth-first Search, Graph
# link: <https://leetcode.com/problems/minimum-height-trees/description/>
# hidden: true
# ---
#
# A tree is an undirected graph in which any two vertices are connected by
# _exactly_  one path. In other words, any connected graph without simple cycles
# is a tree.
#
# Given a tree of `n` nodes labelled from `0` to `n - 1`, and an array of `n -
# 1` `edges` where `edges[i] = [ai, bi]` indicates that there is an undirected
# edge between the two nodes `ai` and `bi` in the tree, you can choose any node
# of the tree as the root. When you select a node `x` as the root, the result
# tree has height `h`. Among all possible rooted trees, those with minimum
# height (i.e. `min(h)`)  are called **minimum height trees** (MHTs).
#
# Return _a list of all **MHTs '** root labels_. You can return the answer in
# **any order**.
#
# The **height** of a rooted tree is the number of edges on the longest downward
# path between the root and a leaf.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/e1.jpg)
#
#
#
#     Input: n = 4, edges = [[1,0],[1,2],[1,3]]
#     Output: [1]
#     Explanation: As shown, the height of the tree is 1 when the root is the node with label 1 which is the only MHT.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/e2.jpg)
#
#
#
#     Input: n = 6, edges = [[3,0],[3,1],[3,2],[3,4],[5,4]]
#     Output: [3,4]
#
#
# **Example 3:**
#
#
#
#     Input: n = 1, edges = []
#     Output: [0]
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, edges = [[0,1]]
#     Output: [0,1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 2 * 104`
#   * `edges.length == n - 1`
#   * `0 <= ai, bi < n`
#   * `ai != bi`
#   * All the pairs `(ai, bi)` are distinct.
#   * The given input is **guaranteed** to be a tree and there will be **no repeated** edges.
#
#
## @lc code=start
using LeetCode

function find_min_height_trees(n::Int, edges::Vector{Vector{Int}})
    (n == 1) && return [1]
    neighbor = [Set{Int}() for i in 1:n]
    for (u, v) in edges
        push!(neighbor[u], v)
        push!(neighbor[v], u)
    end
    q = [i for i in 1:n if length(neighbor[i]) == 1]
    dep_q = ones(Int, size(q))
    m_dep = 1
    visited = fill(false, n)
    visited[q] .= true
    res = q[:]
    while !isempty(q)
        u, dep = popfirst!(q), popfirst!(dep_q)
        for nei in neighbor[u]
            visited[nei] && continue
            visited[nei] = true
            push!(q, nei)
            push!(dep_q, dep + 1)
        end
        if dep > m_dep
            m_dep = dep
            res = [u; q]
        end
    end
    return res
end
## @lc code=end


text:
module PairsDB

using MAT

# actual cycle duration of Daniel's data
const CYCLE_DURATION = 0.24997832635983264

const db_dir = joinpath(@__DIR__, "..", "data")

export DB, get_database, get_data, get_type, get_param,
    get_uniform_param, get_ids
# ============================================================================ #
struct DB
    db::Vector{Pair{Int16, Vector{String}}}
    typ::String
end

function Base.getindex(db::DB, k::Integer)
    k < 1 || k > length(db.db) && error("Index \"$(k)\" is out of range")
    return db.db[k]
end

function Base.getindex(db::DB; id::Integer=1)
    for k in eachindex(db.db)
        if db.db[k].first == id
            return db.db[k]
        end
    end
    return nothing
end

Base.IteratorSize(db::DB) = Base.HasLength()
Base.IteratorEltype(db::DB) = Base.HasEltype()
@inline Base.length(db::DB) = length(db.db)
@inline Base.eltype(db::DB) = eltype(db.db)

Base.iterate(db::DB) = iterate(db, 1)
function Base.iterate(db::DB, state::Integer)
    state > length(db) && return nothing
    return db.db[state], state + 1
end

@inline get_type(db::DB) = db.typ
@inline get_ids(db::DB) = first.(db)
# ============================================================================ #
function get_uniform_param(pair::Pair{Int16,Vector{String}}, name::String)
    out = get_param(pair, name)
    if !all(isequal(out[1]), out[2:end])
        val = NaN
        mx = -Inf
        for x in unique(out)
            n = sum(isequal(x), out)
            if n > mx
                val = x
                mx = n
            end
        end
        k = findfirst(isequal(val), out)
        @warn("Values of \"$(name)\" are inconsistenct across files\n   values: $(out)")
    else
        k = 1
    end
    return out[k]
end
# ============================================================================ #
function get_param(pair::Pair{Int16,Vector{String}}, name::String)
    p = get_param(pair.second[1], name)
    p == nothing && error("Parameter name \"$(name)\" is not valid")
    out = Vector{typeof(p)}(undef, length(pair.second))
    out[1] = p
    for k in 2:length(pair.second)
        out[k] = get_param(pair.second[k], name)
    end
    return out
end
# ============================================================================ #
function get_param(ifile::String, name::String)
    par = matopen(ifile, "r") do io
        return read(io, "parameters")
    end
    return haskey(par, name) ? par[name] : nothing
end
# ============================================================================ #
"""
## get\\_database(typ::String, condition::Function=id->true)
    * type: one of ["contrast", "msequence", "area", "spatial_frequency"]
    * condition(id::Int16) -> (true|false): whether or not to include a pair
    based on their id
"""
function get_database(typ::String, condition::Function=x->true)

    typ = lowercase(typ)

    files = find_files(db_dir, Regex("\\d{8}_\\d{3}_" * typ * "\\-\\d{3}\\.mat"))

    isempty(files) && error("No files match pattern \"$(typ)\"")

    ids = Vector{Int16}(undef, length(files))
    for k in eachindex(files)
        m = match(r"\d{8}_(?<id>\d{3})_.*", files[k])
        m == nothing && error("Failed to extract id from \"$(files[k])\"")
        ids[k] = parse(Int16, m[:id])
    end

    db = Vector{Pair{Int16, Vector{String}}}(undef, 0)
    processed = Vector{Int16}(undef, 0)
    for id in filter(condition, sort(ids))
        k = findall(isequal(id), ids)
        if !(id in processed)
            push!(db, id => files[k])
            push!(processed, id)
        end
    end
    return DB(db, typ)
end
# ============================================================================ #
function trial_indicies(ts::Vector{Float64}, evt::Real, dur::Real)
    k1 = findfirst(x-> x >= evt, ts)
    k1 == nothing && return 1:0

    k2 = findlast(x-> x <= (evt + dur), ts[k1:end])
    k2 == nothing && return 1:0 #k1:length(ts)

    return k1:(k2 + k1 -1)
end
# ============================================================================ #
function trial_timestamps(ret::Vector{Float64}, lgn::Vector{Float64},
    evt::Vector{Float64}, dur::Real, offset::Real=2.0, start::Real=0.0, t0::Real=0.0)

    ret_out = Vector{Float64}(undef, 0)
    lgn_out = Vector{Float64}(undef, 0)
    evt_out = Vector{Float64}(undef, length(evt))

    last_end = start
    for k in eachindex(evt)
        # NOTE: in Rathbun's data the grating phase was *NOT* reset between
        # trials, so we apply another offset to each spike (drift_offset) that
        # correct each spike time so that all are relative to the same grating
        # phase (whatever phase the grating started at) before adding the trial
        # offset and run offset, if t0 is not specified no correction is made
        drift_offset = t0 > 0.0 ? mod(evt[k] .- t0, CYCLE_DURATION) : 0.0

        evt_out[k] = last_end
        ridx = trial_indicies(ret, evt[k], dur)
        lidx = trial_indicies(lgn, evt[k], dur)
        if !isempty(ridx)
            append!(ret_out, (ret[ridx] .- evt[k]) .+ (last_end + drift_offset))
        end
        if !isempty(lidx)
            append!(lgn_out, (lgn[lidx] .- evt[k]) .+ (last_end + drift_offset))
        end
        lmx = isempty(lgn_out) ? 0.0 : lgn_out[end]
        rmx = isempty(ret_out) ? 0.0 : ret_out[end]
        last_end = max(rmx, lmx) + offset
    end
    return ret_out, lgn_out, evt_out
end
# ============================================================================ #
"""
## get\\_data(db::DB, idx::Integer=0; id::Integer=-1, cond::Function=x->true)
    * db::DB - database object
    * idx::Integer [0] - index of pair to get data for (trumps id option)
### Options:
    * id::Integer [-1] - id of pair to get data for
    * ffile::Function [x->true] - filter condition for selecting files based on their parameter dict
    * ftrial::Function [x->true] - filter condition for selecting trials based on their stimulus value

"""
function get_data(db::DB, idx::Integer=0; id::Integer=-1, ffile::Function=x->true, ftrial::Function=x->true)
    offset = 2.0

    idx < 1 && id < 0 && error("You *MUST* provide either an index or a pair id!")

    if idx < 1
        files = db[id=id].second
    else
        files = db[idx].second
        id = db[idx].first
    end

    ret = Vector{Float64}(undef, 0)
    lgn = Vector{Float64}(undef, 0)
    evt = Vector{Float64}(undef, 0)
    lab = Vector{Float64}(undef, 0)

    last_end = 0.0

    for file in files
        matopen(file, "r") do mf
            par = read(mf, "parameters")

            if ffile(par)
                values = read(mf, "values")

                if isempty(values)
                    append!(ret, read(mf, "retina") .+ last_end)
                    append!(lgn , read(mf, "lgn") .+ last_end)

                    if get_type(db) == "contrast"
                        append!(lab, par["contrast"])
                    elseif get_type(db) == "msequence"
                        append!(evt, read(mf, "stimulus") .+ last_end)
                        append!(lab, par["msequence"])
                    end

                    mx = max(ret[end], lgn[end])
                    mx = isempty(evt) ? mx : max(mx, evt[end])
                    last_end =  mx + offset

                else
                    b = ftrial.(values)

                    stim = read(mf, "stimulus")
                    evt_cur = stim[b]

                    if 100 < id < 200
                        # for Rathbun data, supply trial_timestamps with the time of the
                        # very first stimulus onset so it can calculate relative
                        # stimulus phase for the entire run
                        t0 = stim[1]
                    else
                        # not needed for Tucker's data (and does not apply to Marty's)
                        t0 = 0.0
                    end

                    rtmp, ltmp, etmp = trial_timestamps(
                        read(mf, "retina"),
                        read(mf, "lgn"),
                        evt_cur,
                        par["stimulus_duration"],
                        offset,
                        last_end,
                        t0
                    )

                    append!(ret, rtmp)
                    append!(lgn, ltmp)
                    append!(evt, etmp)
                    append!(lab, values[b])

                    re = isempty(ret) ? 0 : ret[end]
                    le = isempty(lgn) ? 0 : lgn[end]
                    ee = isempty(evt) ? 0 : evt[end]
                    last_end = maximum([re, le, ee]) + offset
                end
            end
        end
    end

    return ret, lgn, evt, lab
end
# ============================================================================ #
find_files(dir::AbstractString, re=r".*") = return do_match(dir, re, isfile)
# --------------------------------------------------------------------------- #
find_directories(dir::AbstractString, re=r".*") = return do_match(dir, re, isdir)
# =========================================================================== #
function do_match(dir::AbstractString, re::Regex, f::Function)
    if !isdir(dir)
        error("Input is not a vaild directory path")
    end
    files = [joinpath(dir, x) for x in readdir(dir)]
    return filter(x->occursin(re, x) && f(x), files)
end
# =========================================================================== #
end


text:
module test_octo_subquery

using Test # @test
using Octo.Adapters.SQL # Schema.model from to_sql SELECT SUM FROM WHERE ORDER BY

struct Album
end
Schema.model(Album, table_name="albums")

struct Track
end
Schema.model(Track, table_name="tracks")

albums = from(Album)
tracks = from(Track)

sub = from([SELECT SUM(:bytes) FROM tracks WHERE tracks.AlbumId == albums.AlbumId])
@test sub isa Octo.SubQuery
@test to_sql(sub) == "(SELECT SUM(bytes) FROM tracks WHERE AlbumId = AlbumId)"
q = [SELECT (:albumid, :title) FROM albums WHERE 10000000 > sub ORDER BY :title]
@test to_sql(q) == "SELECT albumid, title FROM albums WHERE 10000000 > (SELECT SUM(bytes) FROM tracks WHERE AlbumId = AlbumId) ORDER BY title"

sub = from([SELECT SUM(:bytes) FROM tracks WHERE tracks.AlbumId == albums.AlbumId], :sub)
@test to_sql(sub) == "(SELECT SUM(bytes) FROM tracks WHERE AlbumId = AlbumId) AS sub"
q = [SELECT (:albumid, :title) FROM albums WHERE 10000000 > sub ORDER BY :title]
@test to_sql(q) == "SELECT albumid, title FROM albums WHERE 10000000 > (SELECT SUM(bytes) FROM tracks WHERE AlbumId = AlbumId) AS sub ORDER BY title"

# https://youtu.be/baOAbOdcnxs?t=1948
struct SA
end
struct SB
end
struct SC
end
Schema.model(SA, table_name="A")
Schema.model(SB, table_name="B")
Schema.model(SC, table_name="C")
❔ = Octo.PlaceHolder
@test to_sql([1 + ❔]) == "1 + ?"
@test to_sql([❔ + 1]) == "? + 1"
A = from(SA, :A)
B = from(SB, :B)
C = from(SC, :C)
sub = from([SELECT (B.id, COUNT(*)) FROM B WHERE B.val == ❔ + 1 GROUP BY B.id], :B)
q = [SELECT * FROM (A, C, sub) WHERE A.val == 123 AND A.id == C.a_id AND B.id == C.b_id]
@test to_sql(q) == "SELECT * FROM A, C, (SELECT B.id, COUNT(*) FROM B WHERE B.val = ? + 1 GROUP BY B.id) AS B WHERE A.val = 123 AND A.id = C.a_id AND B.id = C.b_id"

end # module test_octo_subquery


text:
module Test02

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test02" begin
    re = Automa.RegExp

    a = re.rep('a')
    b = re.cat('b', re.rep('b'))
    ab = re.cat(a, b)

    a.actions[:enter] = [:enter_a]
    a.actions[:exit] = [:exit_a]
    a.actions[:final] = [:final_a]
    b.actions[:enter] = [:enter_b]
    b.actions[:exit] = [:exit_b]
    b.actions[:final] = [:final_b]
    ab.actions[:enter] = [:enter_re]
    ab.actions[:exit] = [:exit_re]
    ab.actions[:final] = [:final_re]

    machine = Automa.compile(ab)

    last, actions = Automa.execute(machine, "ab")
    @test last == 0
    @test actions == [:enter_re,:enter_a,:final_a,:exit_a,:enter_b,:final_b,:final_re,:exit_b,:exit_re]

    for generator in (:table, :inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        init_code = Automa.generate_init_code(ctx, machine)
        exec_code = Automa.generate_exec_code(ctx, machine, :debug)
        validate = @eval function (data)
            logger = Symbol[]
            $(init_code)
            p_end = p_eof = lastindex(data)
            $(exec_code)
            return cs == 0, logger
        end
        @test validate(b"b") == (true, [:enter_re,:enter_a,:exit_a,:enter_b,:final_b,:final_re,:exit_b,:exit_re])
        @test validate(b"a") == (false, [:enter_re,:enter_a,:final_a])
        @test validate(b"ab") == (true, [:enter_re,:enter_a,:final_a,:exit_a,:enter_b,:final_b,:final_re,:exit_b,:exit_re])
        @test validate(b"abb") == (true, [:enter_re,:enter_a,:final_a,:exit_a,:enter_b,:final_b,:final_re,:final_b,:final_re,:exit_b,:exit_re])
    end
end

end


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/abc_setup.jl
## (c) 2018 Eric B. Ford & Danley C. Hsu
# Collection of functions which specific ABC simulation parameters

module EvalSysSimModel
  export setup, get_param_vector, get_ss_obs #, evaluate_model
  export gen_data, calc_summary_stats, calc_distance, is_valid, normalize_dirch
  export make_proposal_dist_multidim_beta
  export setup, get_param_vector, get_ss_obs #, evaluate_model
  using ExoplanetsSysSim
  #using KahanSummation
  using SpecialFunctions
  if VERSION < v"0.7"
      using ABC
      using Compat.Statistics
  else
      using ApproximateBayesianComputing
      const ABC = ApproximateBayesianComputing
      using Statistics
      import ApproximateBayesianComputing.CompositeDistributions.CompositeDist
      import ApproximateBayesianComputing.TransformedBetaDistributions.LinearTransformedBeta
  end
  include(joinpath(dirname(pathof(ExoplanetsSysSim)),"..","examples","poisson_tests", "christiansen_func.jl"))

  sim_param_closure = SimParam()
  #sim_param_closure = setup_sim_param_christiansen()
  summary_stat_ref_closure =  CatalogSummaryStatistics()

    function is_valid(param_vector::Vector{Float64})
      global sim_param_closure
      update_sim_param_from_vector!(param_vector,sim_param_closure)
      rate_tab::Array{Float64,2} = get_any(sim_param_closure, "obs_par", Array{Float64,2})
      #lambda = sum_kbn(rate_tab)  # TODO: Restore kbn
      lambda = sum(rate_tab)
      #if lambda > 10. || any(x -> x < 0., rate_tab)
      if any(x -> x < 0., rate_tab[2:size(rate_tab,1),:]) || any(x -> x > 3., rate_tab[1,:])# || lambda > 10.
         return false
      end
      return true
    end

function normalize_dirch(param_vector::Vector{Float64})
    global sim_param_closure
    p_dim = length(get_any(sim_param_closure, "p_lim_arr", Array{Float64,1}))-1
    r_dim = length(get_any(sim_param_closure, "r_lim_arr", Array{Float64,1}))-1

    for i in 1:p_dim
        param_vector[((i-1)*(r_dim+1)+2):((i-1)*(r_dim+1)+(r_dim+1))] ./= sum(param_vector[((i-1)*(r_dim+1)+2):((i-1)*(r_dim+1)+(r_dim+1))])
    end

    update_sim_param_from_vector!(param_vector,sim_param_closure)
    return param_vector
end

# https://en.wikipedia.org/wiki/Trigamma_function
function trigamma_x_gr_4(x::T) where T<: Real
   1/x + 0.5/x^2 + 1/(6*x^3) - 1/(30*x^5) + 1/(42*x^7) - 1/(30*x^9) + 5/(66*x^11) - 691/(2730*x^13) + 7/(6*x^15)
end

function trigamma_x_lt_4(x::T) where T<: Real
  n = floor(Int64,5-x)
  z = x+n
  val = trigamma_x_gr_4(z)
  for i in 1:n
    z -= 1
    val += 1/z^2
  end
  val
end

function trigamma(x::T) where T<: Real
   x >= 4 ? trigamma_x_gr_4(x) : trigamma_x_lt_4(x)
end

function make_proposal_dist_multidim_beta(theta::AbstractArray{Float64,2}, weights::AbstractArray{Float64,1},  tau_factor::Float64; verbose::Bool = false)
    global sim_param_closure
    p_dim = length(get_any(sim_param_closure, "p_lim_arr", Array{Float64,1}))-1
    r_dim = length(get_any(sim_param_closure, "r_lim_arr", Array{Float64,1}))-1
    max_col_rate = 3.0

    function mom_alpha(x_bar::T, v_bar::T) where T<: Real
        x_bar * (((x_bar * (1 - x_bar)) / v_bar) - 1)
    end
    function mom_beta(x_bar::T, v_bar::T) where T<: Real
        (1 - x_bar) * (((x_bar * (1 - x_bar)) / v_bar) - 1)
    end
    # For algorithm, see https://scholarsarchive.byu.edu/cgi/viewcontent.cgi?article=2613&context=etd
    function fit_beta_mle(x::AbstractArray{T,1}; tol::T = 1e-6, max_it::Int64 = 10, init_guess::AbstractArray{T,1} = Array{T}(undef,0), w::AbstractArray{T,1} = Array{T}(undef,0), verbose::Bool = false ) where T<: Real
        lnxbar =   length(w)>1 ? Compat.Statistics.mean(log.(x),AnalyticWeights(w)) : Compat.Statistics.mean(log.(x))
        ln1mxbar = length(w)>1 ? Compat.Statistics.mean(log.(1.0.-x),AnalyticWeights(w)) : Compat.Statistics.mean(log.(1.0.-x))

        function itterate( mle_guess::Vector{T} ) where T<:Real
            (alpha, beta) = (mle_guess[1], mle_guess[2])
            dgab = digamma(alpha+beta)
            g1 = dgab - digamma(alpha) + lnxbar
            g2 = dgab - digamma(beta) + ln1mxbar
            tgab = trigamma(alpha+beta)
            G = [dgab-trigamma(alpha) tgab; tgab tgab-trigamma(beta)]
            mle_guess -= G \ [g1, g2]
        end

        local mle_new
        if length(init_guess) != 2
            xbar = length(w)>1 ? Compat.Statistics.mean(x,AnalyticWeights(w)) : Compat.Statistics.mean(x)
            vbar = length(w)>1 ? Compat.Statistics.varm(x,xbar,AnalyticWeights(w)) : Compat.Statistics.varm(x,xbar)
            mle_new = (vbar < xbar*(1.0-xbar)) ? [mom_alpha(xbar, vbar), mom_beta(xbar,vbar)] : ones(T,2)
        else
            mle_new = init_guess
        end
        if verbose
            println("it = 0: ", mle_new)
        end
        if any(mle_new.<=zero(T))
            println("# Warning: mean= ", xbar, " var= ",var," (alpha,beta)_init= ",mle_new," invalid, reinitializing to (1,1)")
            verbose = true
            mle_new = ones(T,2)
        end
        for i in 1:max_it
            mle_old = mle_new
            mle_new = itterate( mle_old )
            epsilon = max(abs.(mle_old.-mle_new))
            if verbose
                println("# it = ", i, ": ", mle_new, " max(Delta alpha, Delta beta)= ", epsilon)
            end
            if epsilon < tol
                break
            end
        end
        return mle_new
    end
    function make_beta(x::AbstractArray{T,1}, w::AbstractArray{T,1};
                       mean::T = Compat.Statistics.mean(x,AnalyticWeights(w)),
                       var::T = Compat.Statistics.varm(x,xbar,AnalyticWeights(w)) ) where T<:Real
        alpha_beta = (var < mean*(1.0-mean)) ? [mom_alpha(mean, var), mom_beta(mean,var)] : ones(T,2)
        if any(alpha_beta.<=zero(T))
            alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
        end
        if any(alpha_beta.<=zero(T))
            alpha_beta = ones(T,2)
        end
        Beta(alpha_beta[1], alpha_beta[2])
    end
    function make_beta_transformed(x::AbstractArray{T,1}, w::AbstractArray{T,1}; xmin::T=zero(T), xmax::T=one(T), mean::T = Compat.Statistics.mean(x,AnalyticWeights(w)), var::T = Compat.Statistics.varm(x,xbar,AnalyticWeights(w)) ) where T<:Real
        alpha_beta = (var < mean*(1.0-mean)) ? [mom_alpha(mean, var), mom_beta(mean,var)] : ones(T,2)
        if any(alpha_beta.<=zero(T))
            alpha_beta = fit_beta_mle(x, w=w, init_guess=alpha_beta, verbose=true)
        end
        if any(alpha_beta.<=zero(T))
            alpha_beta = ones(T,2)
        end
        LinearTransformedBeta(alpha_beta[1], alpha_beta[2], xmin=xmin, xmax=xmax)
    end

    theta_mean =  sum(theta.*weights',2) # weighted mean for parameters
    tau = tau_factor*ABC.var_weighted(theta'.-theta_mean',weights)  # scaled, weighted covar for parameters

    #=
    println("mean= ",theta_mean)
    println("var= ",tau)
    for i in 1:length(theta_mean)
        println("a= ",alpha(theta_mean[i],tau[i]), "  b= ",beta(theta_mean[i],tau[i]))
    end
    =#

    dist_arr = ContinuousDistribution[]
    for j in 1:p_dim
        col_startidx = (j-1)*(r_dim+1)+1
        dist_arr = vcat(dist_arr, make_beta_transformed(theta[col_startidx,:], weights, xmin=0.0, xmax=max_col_rate, mean=theta_mean[col_startidx]/max_col_rate, var=tau[col_startidx]/max_col_rate^2), ContinuousDistribution[ make_beta(theta[i,:], weights, mean=theta_mean[i], var=tau[i]) for i in (col_startidx+1):(col_startidx+r_dim)]   )
    end

    dist = CompositeDist(dist_arr)
end

function make_proposal_dist_multidim_beta(pop::ABC.abc_population_type, tau_factor::Float64; verbose::Bool = false)
    make_proposal_dist_multidim_beta(pop.theta, pop.weights, tau_factor, verbose=verbose)
end

function gen_data(param_vector::Vector{Float64})
    global sim_param_closure
    update_sim_param_from_vector!(param_vector,sim_param_closure)
    cat_phys = generate_kepler_physical_catalog(sim_param_closure)
    #cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys, sim_param_closure)
    #cat_obs = ExoplanetsSysSim.observe_kepler_targets_single_obs(cat_phys_cut, sim_param_closure)
    cat_obs = ExoplanetsSysSim.observe_kepler_targets_sky_avg(cat_phys, sim_param_closure)
    return cat_obs
end

# TODO OPT: Eventually, could adapt ABC.jl to use distance from first pass to decide if should compute additional summary statistics
function calc_summary_stats(cat::KeplerObsCatalog)
    global sim_param_closure
    sum_stat = calc_summary_stats_obs_binned_rates(cat, sim_param_closure, obs_skyavg = true)
    return sum_stat
end

function calc_distance(sum_stat_obs::CatalogSummaryStatistics,sum_stat_sim::CatalogSummaryStatistics, n::Integer = 0)
    global sim_param_closure
    dist1 = calc_distance_vector_binned(sum_stat_obs,sum_stat_sim, 1, sim_param_closure)
    num_available = length(dist1)
    num_to_use = n>0 ? min(n,num_available) : num_available
    return calc_scalar_distance(dist1[1:num_to_use])
end

function setup()
    global sim_param_closure = setup_sim_param_christiansen()
    sim_param_closure = set_test_param(sim_param_closure)

    ### Use simulated planet candidate catalog data
    #add_param_fixed(sim_param_closure,"num_kepler_targets",150000)  # For "observed" catalog
    #cat_obs = simulated_read_kepler_observations(sim_param_closure)
    ###

    WindowFunction.setup(sim_param_closure)
    println("# Finished reading in window function data")
    # add_param_fixed(sim_param_closure,"win_func_data",ExoplanetsSysSim.WindowFunction.get_window_function_data() )  # TODO OPT DETAIL: Does this serve any purpose?

    ### Use real planet candidate catalog data
    df_star = setup_star_table_christiansen(sim_param_closure)
    println("# Finished reading in stellar data")
    df_koi,usable_koi = read_koi_catalog(sim_param_closure)
    println("# Finished reading in KOI data")
    cat_obs = setup_actual_planet_candidate_catalog(df_star, df_koi, usable_koi, sim_param_closure)
    ###

    global summary_stat_ref_closure = calc_summary_stats_obs_binned_rates(cat_obs,sim_param_closure, trueobs_cat = true)
end

get_param_vector() = make_vector_of_sim_param(sim_param_closure)
get_ss_obs() = summary_stat_ref_closure

function set_simparam_ss(sim_param::ExoplanetsSysSim.SimParam, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics)
    global sim_param_closure = sim_param
    global summary_stat_ref_closure = ss_true
end

end  # module EvalSysSimModel

#include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))

module SysSimABC
  import ExoplanetsSysSim
  export setup_abc, run_abc, run_abc_largegen
  #import EvalSysSimModel
  using ..EvalSysSimModel
  import Distributions
  if VERSION < v"0.7"
      using ABC
      #using Compat
      using Compat.Statistics
  else
      using ApproximateBayesianComputing
      ABC = ApproximateBayesianComputing
      using Statistics
      using ApproximateBayesianComputing.CompositeDistributions
  end
  using Compat
  using Distributed

  include(joinpath(dirname(pathof(ExoplanetsSysSim)),"..","examples","poisson_tests", "christiansen_func.jl"))

  function setup_abc(num_dist::Integer = 0; max_generations::Int64=100 )
    EvalSysSimModel.setup()
    theta_true = EvalSysSimModel.get_param_vector()
    #param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 0.3) for x in 1:length(theta_true)] )
    limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    r_dim = length(get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1}))-1
    prior_arr = ContinuousDistribution[]
    for i in 1:(length(limitP)-1)
        max_in_col = floor(3*log(limitP[i+1]/limitP[i])/log(2))
        lambda_col = Uniform(0.0, max_in_col)
        dirch = Dirichlet(ones(r_dim))
        prior_arr = vcat(prior_arr, [lambda_col, dirch])
    end
    param_prior = CompositeDist(prior_arr)
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, make_proposal_dist=EvalSysSimModel.make_proposal_dist_multidim_beta, is_valid=EvalSysSimModel.is_valid, normalize=EvalSysSimModel.normalize_dirch, num_part=200, num_max_attempt=50, num_max_times=max_generations, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel, adaptive_quantiles = false, epsilon_reduction_factor=0.9, tau_factor=2.0);
end

function run_abc_largegen(pop::ABC.abc_population_type, ss_true::ExoplanetsSysSim.CatalogSummaryStatistics, epshist_targ::Float64, npart::Integer = 1000, num_dist::Integer = 0)
    sim_param_closure = setup_sim_param_christiansen()
    sim_param_closure = set_test_param(sim_param_closure)
    setup_star_table_christiansen(sim_param_closure)
    EvalSysSimModel.set_simparam_ss(sim_param_closure, ss_true)

    theta_true = EvalSysSimModel.get_param_vector()
    #param_prior = CompositeDist( Distributions.ContinuousDistribution[Distributions.Uniform(0., 0.3) for x in 1:length(theta_true)] )
    limitP::Array{Float64,1} = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
    r_dim = length(get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1}))-1
    prior_arr = ContinuousDistribution[]
    for i in 1:(length(limitP)-1)
        max_in_col = floor(3*log(limitP[i+1]/limitP[i])/log(2))
        lambda_col = Uniform(0.0, max_in_col)
        dirch = Dirichlet(ones(r_dim))
        prior_arr = vcat(prior_arr, [lambda_col, dirch])
    end
    param_prior = CompositeDist(prior_arr)
    in_parallel = nworkers() > 1 ? true : false

    calc_distance_ltd(sum_stat_obs::ExoplanetsSysSim.CatalogSummaryStatistics,sum_stat_sim::ExoplanetsSysSim.CatalogSummaryStatistics) = EvalSysSimModel.calc_distance(sum_stat_obs,sum_stat_sim,num_dist)

    global abc_plan = ABC.abc_pmc_plan_type(EvalSysSimModel.gen_data,EvalSysSimModel.calc_summary_stats,calc_distance_ltd, param_prior, make_proposal_dist=EvalSysSimModel.make_proposal_dist_multidim_beta, is_valid=EvalSysSimModel.is_valid, normalize=EvalSysSimModel.normalize_dirch, num_part=npart, num_max_attempt=50, num_max_times=1, epsilon_init=9.9e99, target_epsilon=1.0e-100, in_parallel=in_parallel);

    println("# run_abc_largegen: ",EvalSysSimModel.sim_param_closure)
    #if (std(pop.theta, 2)[1])/(mean(pop.theta, 2)[1]) > 0.3
    #    sampler_largegen = abc_plan.make_proposal_dist(pop, 1.2)
    #elseif (std(pop.theta, 2)[1])/(mean(pop.theta, 2)[1]) > 0.1
    #    sampler_largegen = abc_plan.make_proposal_dist(pop, 1.6)
    #else
    sampler_largegen = abc_plan.make_proposal_dist(pop, abc_plan.tau_factor)
    #end
    theta_largegen = Array{Float64}(size(pop.theta, 1), npart)
    weight_largegen = Array{Float64}(npart)
    for i in 1:npart
        theta_val, dist_largegen, attempts_largegen = ABC.generate_theta(abc_plan, sampler_largegen, ss_true, epshist_targ)
        theta_largegen[:,i] = theta_val
        prior_logpdf = Distributions.logpdf(abc_plan.prior,theta_val)
        sampler_logpdf = ABC.logpdf(sampler_largegen, theta_val)
        weight_largegen[i] = exp(prior_logpdf-sampler_logpdf)
    end
    return theta_largegen, weight_largegen
end

function run_abc(abc_plan::ABC.abc_pmc_plan_type)
    #global sim_param_closure
#    println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    #println("True catalog SS: ", ss_true)
    pop_out = ABC.run_abc(abc_plan,ss_true;verbose=true)
end

function run_abc(abc_plan::ABC.abc_pmc_plan_type, pop::ABC.abc_population_type)
    #global sim_param_closure
    dist_threshold = maximum(pop.dist)
    EvalSysSimModel.add_param_fixed(EvalSysSimModel.sim_param_closure,"minimum ABC dist skip pass 2",dist_threshold)
    #println("# run_abc: ",EvalSysSimModel.sim_param_closure)
    ss_true = EvalSysSimModel.get_ss_obs()
    pop_out = ABC.run_abc(abc_plan,ss_true,pop;verbose=true)
end

end # module SysSimABC


text:
"""
    gc2gd(n, xyz)

Transform geocentric coordinates to geodetic using the specified
reference ellipsoid.

### Given ###

- `n`: Ellipsoid identifier (Note 1)
- `xyz`: Geocentric vector (Note 2)

### Returned ###

- `elong`: Longitude (radians, east +ve, Note 3)
- `phi`: Latitude (geodetic, radians, Note 3)
- `height`: Height above ellipsoid (geodetic, Notes 2,3)

### Notes ###

1. The identifier n is a number that specifies the choice of
   reference ellipsoid.  The following are supported:

        - `WGS84`
        - `GRS80`
        - `WGS72`

2. The geocentric vector (xyz, given) and height (height, returned)
   are in meters.

3. An error status -1 means that the identifier n is illegal.  An
   error status -2 is theoretically impossible.  In all error cases,
   all three results are set to -1e9.

4. The inverse transformation is performed in the function eraGd2gc.

### Called ###

- `eraEform`: Earth reference ellipsoids
- `eraGc2gde`: geocentric to geodetic transformation, general

"""
function gc2gd(n, xyz)
    elong = Ref(0.0)
    phi = Ref(0.0)
    height = Ref(0.0)
    i = ccall((:eraGc2gd, liberfa), Cint,
              (Cint, Ptr{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
              n, xyz, elong, phi, height)
    if i == -1
        throw(ERFAException("illegal identifier"))
    elseif i == -2
        throw(ERFAException("internal error"))
    end
    elong[], phi[], height[]
end

"""
    gc2gde(a, f, xyz)

Transform geocentric coordinates to geodetic for a reference
ellipsoid of specified form.

### Given ###

- `a`: Equatorial radius (Notes 2,4)
- `f`: Flattening (Note 3)
- `xyz`: Geocentric vector (Note 4)

### Returned ###

- `elong`: Longitude (radians, east +ve)
- `phi`: Latitude (geodetic, radians)
- `height`: Height above ellipsoid (geodetic, Note 4)

### Notes ###

1. This function is based on the GCONV2H Fortran subroutine by
   Toshio Fukushima (see reference).

2. The equatorial radius, a, can be in any units, but meters is
   the conventional choice.

3. The flattening, f, is (for the Earth) a value around 0.00335,
   i.e. around 1/298.

4. The equatorial radius, a, and the geocentric vector, xyz,
   must be given in the same units, and determine the units of
   the returned height, height.

5. If an error occurs (status < 0), elong, phi and height are
   unchanged.

6. The inverse transformation is performed in the function
   eraGd2gce.

7. The transformation for a standard ellipsoid (such as ERFA_WGS84) can
   more conveniently be performed by calling eraGc2gd, which uses a
   numerical code to identify the required A and F values.

### Reference ###

- Fukushima, T., "Transformation from Cartesian to geodetic
    coordinates accelerated by Halley's method", J.Geodesy (2006)
    79: 689-693

"""
function gc2gde(a, f, xyz)
    elong = Ref(0.0)
    phi = Ref(0.0)
    height = Ref(0.0)
    i = ccall((:eraGc2gde, liberfa), Cint,
              (Cdouble, Cdouble, Ptr{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
              a, f, xyz, elong, phi, height)
    if i == -1
        throw(ERFAException("illegal f"))
    elseif i == -2
        throw(ERFAException("internal a"))
    end
    elong[], phi[], height[]
end

"""
    gd2gc(n, elong, phi, height)

Transform geodetic coordinates to geocentric using the specified
reference ellipsoid.

### Given ###

- `n`: Ellipsoid identifier (Note 1)
- `elong`: Longitude (radians, east +ve)
- `phi`: Latitude (geodetic, radians, Note 3)
- `height`: Height above ellipsoid (geodetic, Notes 2,3)

### Returned ###

- `xyz`: Geocentric vector (Note 2)

### Notes ###

1. The identifier n is a number that specifies the choice of
   reference ellipsoid.  The following are supported:

      n    ellipsoid

      1     ERFA_WGS84
      2     ERFA_GRS80
      3     ERFA_WGS72

   The n value has no significance outside the ERFA software.  For
   convenience, symbols ERFA_WGS84 etc. are defined in erfam.h.

2. The height (height, given) and the geocentric vector (xyz,
   returned) are in meters.

3. No validation is performed on the arguments elong, phi and
   height.  An error status -1 means that the identifier n is
   illegal.  An error status -2 protects against cases that would
   lead to arithmetic exceptions.  In all error cases, xyz is set
   to zeros.

4. The inverse transformation is performed in the function eraGc2gd.

### Called ###

- `eraEform`: Earth reference ellipsoids
- `eraGd2gce`: geodetic to geocentric transformation, general
- `eraZp`: zero p-vector

"""
function gd2gc(n, elong, phi, height)
    xyz = zeros(3)
    i = ccall((:eraGd2gc, liberfa), Cint,
              (Cint, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
              n, elong, phi, height, xyz)
    if i == -1
        throw(ERFAException("illegal identifier"))
    elseif i == -2
        throw(ERFAException("illegal case"))
    end
    xyz
end

"""
    gd2gce(a, f, elong, phi, height)

Transform geodetic coordinates to geocentric for a reference
ellipsoid of specified form.

### Given ###

- `a`: Equatorial radius (Notes 1,4)
- `f`: Flattening (Notes 2,4)
- `elong`: Longitude (radians, east +ve)
- `phi`: Latitude (geodetic, radians, Note 4)
- `height`: Height above ellipsoid (geodetic, Notes 3,4)

### Returned ###

- `xyz`: Geocentric vector (Note 3)

### Notes ###

1. The equatorial radius, a, can be in any units, but meters is
   the conventional choice.

2. The flattening, f, is (for the Earth) a value around 0.00335,
   i.e. around 1/298.

3. The equatorial radius, a, and the height, height, must be
   given in the same units, and determine the units of the
   returned geocentric vector, xyz.

4. No validation is performed on individual arguments.  The error
   status -1 protects against (unrealistic) cases that would lead
   to arithmetic exceptions.  If an error occurs, xyz is unchanged.

5. The inverse transformation is performed in the function
   eraGc2gde.

6. The transformation for a standard ellipsoid (such as ERFA_WGS84) can
   more conveniently be performed by calling eraGd2gc,  which uses a
   numerical code to identify the required a and f values.

### References ###

- Green, R.M., Spherical Astronomy, Cambridge University Press,
    (1985) Section 4.5, p96.

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992),
    Section 4.22, p202.

"""
function gd2gce(a, f, elong, phi, height)
    xyz = zeros(3)
    i = ccall((:eraGd2gce, liberfa), Cint,
              (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
              a, f, elong, phi, height, xyz)
    if i == -1
        throw(ERFAException("illegal case"))
    end
    xyz
end

"""
    gst06(uta, utb, tta, ttb, rnpb)

Greenwich apparent sidereal time, IAU 2006, given the NPB matrix.

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)
- `tta`, `ttb`: TT as a 2-part Julian Date (Notes 1,2)
- `rnpb`: Nutation x precession x bias matrix

### Returned ###

- Greenwich apparent sidereal time (radians)

### Notes ###

1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

          Part A        Part B

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   rotation angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession-nutation.  If UT1 is
   used for both purposes, errors of order 100 microarcseconds
   result.

3. Although the function uses the IAU 2006 series for s+XY/2, it is
   otherwise independent of the precession-nutation model and can in
   practice be used with any equinox-based NPB matrix.

4. The result is returned in the range 0 to 2pi.

### Called ###

- `eraBpn2xy`: extract CIP X,Y coordinates from NPB matrix
- `eraS06`: the CIO locator s, given X,Y, IAU 2006
- `eraAnp`: normalize angle into range 0 to 2pi
- `eraEra00`: Earth rotation angle, IAU 2000
- `eraEors`: equation of the origins, given NPB matrix and s

### Reference ###

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
function gst06(uta, utb, tta, ttb, rnpb)
    ccall((:eraGst06, liberfa), Cdouble,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
          uta, utb, tta, ttb, rnpb)
end

"""
    gmst82(dj1, dj2)

Universal Time to Greenwich mean sidereal time (IAU 1982 model).

### Given ###

- `dj1`, `dj2`: UT1 Julian Date (see note)

### Returned ###

- Greenwich mean sidereal time (radians)

### Notes ###

1. The UT1 date dj1+dj2 is a Julian Date, apportioned in any
   convenient way between the arguments dj1 and dj2.  For example,
   JD(UT1)=2450123.7 could be expressed in any of these ways,
   among others:

           dj1            dj2

       2450123.7          0          (JD method)
        2451545        -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5         0.2         (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  The date & time method is
   best matched to the algorithm used:  maximum accuracy (or, at
   least, minimum noise) is delivered when the dj1 argument is for
   0hrs UT1 on the day in question and the dj2 argument lies in the
   range 0 to 1, or vice versa.

2. The algorithm is based on the IAU 1982 expression.  This is
   always described as giving the GMST at 0 hours UT1.  In fact, it
   gives the difference between the GMST and the UT, the steady
   4-minutes-per-day drawing-ahead of ST with respect to UT.  When
   whole days are ignored, the expression happens to equal the GMST
   at 0 hours UT1 each day.

3. In this function, the entire UT1 (the sum of the two arguments
   dj1 and dj2) is used directly as the argument for the standard
   formula, the constant term of which is adjusted by 12 hours to
   take account of the noon phasing of Julian Date.  The UT1 is then
   added, but omitting whole days to conserve accuracy.

### Called ###

- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Transactions of the International Astronomical Union,
    XVIII B, 67 (1983).

- Aoki et al., Astron. Astrophys. 105, 359-361 (1982).

"""
gmst82

"""
    gst00b(dr, dd)

Greenwich apparent sidereal time (consistent with IAU 2000
resolutions but using the truncated nutation model IAU 2000B).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)

### Returned ###

- Greenwich apparent sidereal time (radians)

### Notes ###

1. The UT1 date uta+utb is a Julian Date, apportioned in any
   convenient way between the argument pair.  For example,
   JD=2450123.7 could be expressed in any of these ways, among
   others:

           uta            utb

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. The result is compatible with the IAU 2000 resolutions, except
   that accuracy has been compromised for the sake of speed and
   convenience in two respects:

   . UT is used instead of TDB (or TT) to compute the precession
     component of GMST and the equation of the equinoxes.  This
     results in errors of order 0.1 mas at present.

   . The IAU 2000B abridged nutation model (McCarthy & Luzum, 2001)
     is used, introducing errors of up to 1 mas.

3. This GAST is compatible with the IAU 2000 resolutions and must be
   used only in conjunction with other IAU 2000 compatible
   components such as precession-nutation.

4. The result is returned in the range 0 to 2pi.

5. The algorithm is from Capitaine et al. (2003) and IERS
   Conventions 2003.

### Called ###

- `eraGmst00`: Greenwich mean sidereal time, IAU 2000
- `eraEe00b`: equation of the equinoxes, IAU 2000B
- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- McCarthy, D.D. & Luzum, B.J., "An abridged model of the
    precession-nutation of the celestial pole", Celestial Mechanics &
    Dynamical Astronomy, 85, 37-49 (2003)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
gst00b

"""
    gst94(dr, dd)

Greenwich apparent sidereal time (consistent with IAU 1982/94
resolutions).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)

### Returned ###

- Greenwich apparent sidereal time (radians)

### Notes ###

1. The UT1 date uta+utb is a Julian Date, apportioned in any
   convenient way between the argument pair.  For example,
   JD=2450123.7 could be expressed in any of these ways, among
   others:

           uta            utb

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. The result is compatible with the IAU 1982 and 1994 resolutions,
   except that accuracy has been compromised for the sake of
   convenience in that UT is used instead of TDB (or TT) to compute
   the equation of the equinoxes.

3. This GAST must be used only in conjunction with contemporaneous
   IAU standards such as 1976 precession, 1980 obliquity and 1982
   nutation.  It is not compatible with the IAU 2000 resolutions.

4. The result is returned in the range 0 to 2pi.

### Called ###

- `eraGmst82`: Greenwich mean sidereal time, IAU 1982
- `eraEqeq94`: equation of the equinoxes, IAU 1994
- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

- IAU Resolution C7, Recommendation 3 (1994)

"""
gst94

for name in ("gmst82",
             "gst00b",
             "gst94")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d1, d2) = ccall(($fc, liberfa), Cdouble, (Cdouble, Cdouble), d1, d2)
end

"""
    gmst00(uta, utb, tta, ttb)

Greenwich mean sidereal time (model consistent with IAU 2000
resolutions).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)
- `tta`, `ttb`: TT as a 2-part Julian Date (Notes 1,2)

### Returned ###

- Greenwich mean sidereal time (radians)

### Notes ###

1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

          Part A         Part B

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession.  If UT1 is used for
   both purposes, errors of order 100 microarcseconds result.

3. This GMST is compatible with the IAU 2000 resolutions and must be
   used only in conjunction with other IAU 2000 compatible
   components such as precession-nutation and equation of the
   equinoxes.

4. The result is returned in the range 0 to 2pi.

5. The algorithm is from Capitaine et al. (2003) and IERS
   Conventions 2003.

### Called ###

- `eraEra00`: Earth rotation angle, IAU 2000
- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
gmst00

"""
    gmst06(uta, utb, tta, ttb)

Greenwich mean sidereal time (consistent with IAU 2006 precession).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)
- `tta`, `ttb`: TT as a 2-part Julian Date (Notes 1,2)

### Returned ###

- Greenwich mean sidereal time (radians)

### Notes ###

1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

          Part A        Part B

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   rotation angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession.  If UT1 is used for
   both purposes, errors of order 100 microarcseconds result.

3. This GMST is compatible with the IAU 2006 precession and must not
   be used with other precession models.

4. The result is returned in the range 0 to 2pi.

### Called ###

- `eraEra00`: Earth rotation angle, IAU 2000
- `eraAnp`: normalize angle into range 0 to 2pi

### Reference ###

- Capitaine, N., Wallace, P.T. & Chapront, J., 2005,
    Astron.Astrophys. 432, 355

"""
gmst06

"""
    gst00a(uta, utb, tta, ttb)

Greenwich apparent sidereal time (consistent with IAU 2000
resolutions).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)
- `tta`, `ttb`: TT as a 2-part Julian Date (Notes 1,2)

### Returned  ###

- Greenwich apparent sidereal time (radians)

### Notes ###

1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

          Part A        Part B

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   Rotation Angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession-nutation.  If UT1 is
   used for both purposes, errors of order 100 microarcseconds
   result.

3. This GAST is compatible with the IAU 2000 resolutions and must be
   used only in conjunction with other IAU 2000 compatible
   components such as precession-nutation.

4. The result is returned in the range 0 to 2pi.

5. The algorithm is from Capitaine et al. (2003) and IERS
   Conventions 2003.

### Called ###

- `eraGmst00`: Greenwich mean sidereal time, IAU 2000
- `eraEe00a`: equation of the equinoxes, IAU 2000A
- `eraAnp`: normalize angle into range 0 to 2pi

### References ###

- Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
    implement the IAU 2000 definition of UT1", Astronomy &
    Astrophysics, 406, 1135-1149 (2003)

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
gst00a

"""
    gst06a(uta, utb, tta, ttb)

Greenwich apparent sidereal time (consistent with IAU 2000 and 2006
resolutions).

### Given ###

- `uta`, `utb`: UT1 as a 2-part Julian Date (Notes 1,2)
- `tta`, `ttb`: TT as a 2-part Julian Date (Notes 1,2)

### Returned ###

- Greenwich apparent sidereal time (radians)

### Notes ###

1. The UT1 and TT dates uta+utb and tta+ttb respectively, are both
   Julian Dates, apportioned in any convenient way between the
   argument pairs.  For example, JD=2450123.7 could be expressed in
   any of these ways, among others:

          Part A        Part B

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable (in the case of UT;  the TT is not at all critical
   in this respect).  The J2000 and MJD methods are good compromises
   between resolution and convenience.  For UT, the date & time
   method is best matched to the algorithm that is used by the Earth
   rotation angle function, called internally:  maximum precision is
   delivered when the uta argument is for 0hrs UT1 on the day in
   question and the utb argument lies in the range 0 to 1, or vice
   versa.

2. Both UT1 and TT are required, UT1 to predict the Earth rotation
   and TT to predict the effects of precession-nutation.  If UT1 is
   used for both purposes, errors of order 100 microarcseconds
   result.

3. This GAST is compatible with the IAU 2000/2006 resolutions and
   must be used only in conjunction with IAU 2006 precession and
   IAU 2000A nutation.

4. The result is returned in the range 0 to 2pi.

### Called ###

- `eraPnm06a`: classical NPB matrix, IAU 2006/2000A
- `eraGst06`: Greenwich apparent ST, IAU 2006, given NPB matrix

### Reference ###

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
gst06a

for name in ("gmst00",
             "gmst06",
             "gst00a",
             "gst06a")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d1, d2, t1, t2) = ccall(($fc, liberfa), Cdouble, (Cdouble, Cdouble, Cdouble, Cdouble), d1, d2, t1, t2)
end

"""
    g2icrs(dl, db)

Transformation from Galactic Coordinates to ICRS.

### Given ###

- `dl`: Galactic longitude (radians)
- `db`: Galactic latitude (radians)

### Returned ###

- `dr`: ICRS right ascension (radians)
- `dd`: ICRS declination (radians)

### Notes ###

1. The IAU 1958 system of Galactic coordinates was defined with
   respect to the now obsolete reference system FK4 B1950.0.  When
   interpreting the system in a modern context, several factors have
   to be taken into account:

   . The inclusion in FK4 positions of the E-terms of aberration.

   . The distortion of the FK4 proper motion system by differential
     Galactic rotation.

   . The use of the B1950.0 equinox rather than the now-standard
     J2000.0.

   . The frame bias between ICRS and the J2000.0 mean place system.

   The Hipparcos Catalogue (Perryman & ESA 1997) provides a rotation
   matrix that transforms directly between ICRS and Galactic
   coordinates with the above factors taken into account.  The
   matrix is derived from three angles, namely the ICRS coordinates
   of the Galactic pole and the longitude of the ascending node of
   the galactic equator on the ICRS equator.  They are given in
   degrees to five decimal places and for canonical purposes are
   regarded as exact.  In the Hipparcos Catalogue the matrix
   elements are given to 10 decimal places (about 20 microarcsec).
   In the present ERFA function the matrix elements have been
   recomputed from the canonical three angles and are given to 30
   decimal places.

2. The inverse transformation is performed by the function eraIcrs2g.

### Called ###

- `eraAnp`: normalize angle into range 0 to 2pi
- `eraAnpm`: normalize angle into range +/- pi
- `eraS2c`: spherical coordinates to unit vector
- `eraTrxp`: product of transpose of r-matrix and p-vector
- `eraC2s`: p-vector to spherical

### Reference ###

- Perryman M.A.C. & ESA, 1997, ESA SP-1200, The Hipparcos and Tycho
    catalogues.  Astrometric and photometric star catalogues
    derived from the ESA Hipparcos Space Astrometry Mission.  ESA
    Publications Division, Noordwijk, Netherlands.

"""
function g2icrs(a, b)
    r1 = Ref(0.0)
    r2 = Ref(0.0)
    ccall((:eraG2icrs, liberfa), Cvoid,
            (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
            a, b, r1, r2)
    r1[], r2[]
end


text:
module Gcd

export my_gcd

"""
   my_gcd(a, b)

Compute the greatest common divisor of a and b

# Arguments
- `a::Integer`: first argument, not required in real life documentation
- `b::Integer`: second argument, not required in real life documentation

Both arguments should be strictly positive, otherwise a `DomaiunErr` exception will be thrown.

# Examples

```julia-repl
julia> my_gcd(12, 15)
3
julia> my_gcd(13, 11)
1
```
"""
function my_gcd(a::Integer, b::Integer)
    a < 1 && throw(DomainError(a, "argument should be strictly positive"))
    b < 1 && throw(DomainError(b, "argument should be strictly positive"))
    while a != b
        if a > b
            a -= b
        else
            b -= a
        end
    end
    return a
end

end # module


text:
# Polynomial kernel function

"""
    Poly <: Kernel

Polynomial kernel (covariance)
```math
k(x,x') = σ²(xᵀx' + c)ᵈ
```
with signal standard deviation ``σ``, additive constant ``c``, and degree ``d``.
"""
mutable struct Poly{T<:Real} <: Kernel
    "Constant"
    c::T
    "Signal variance"
    σ2::T
    "Degree of polynomial"
    deg::Int
    "Priors for kernel parameters"
    priors::Array
end

"""
Polynomial kernel function

    Poly(lc::Real, lσ::Real, deg::Int)

# Arguments
  - `lc::Real`: additive constant (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
  - `deg::Int`: degree of polynomial
"""
Poly(lc::T, lσ::T, deg::Int) where T = Poly{T}(exp(lc), exp(2 * lσ), deg, [])

function KernelData(k::Poly, X1::AbstractMatrix, X2::AbstractMatrix)
    XtX=X1'*X2
	if X1==X2
    	LinearAlgebra.copytri!(XtX, 'U')
	end
    LinIsoData(XtX)
end
kernel_data_key(k::Poly, X1::AbstractMatrix, X2::AbstractMatrix) = "LinIsoData"

_cov(poly::Poly, xTy) = poly.σ2*(poly.c.+xTy).^poly.deg
function cov(poly::Poly, x::AbstractVector, y::AbstractVector)
    K = _cov(poly, dot(x,y))
end
function cov!(cK::AbstractMatrix, poly::Poly, X::AbstractMatrix, data::LinIsoData)
    cK .= _cov(poly, data.XtX)
end
@inline @inbounds function cov_ij(poly::Poly, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, dim::Int)
    return _cov(poly, dotij(X1, X2, i, j, dim))
end
@inline @inbounds function cov_ij(poly::Poly, X1::AbstractMatrix, X2::AbstractMatrix, data::LinIsoData, i::Int, j::Int, dim::Int)
    return _cov(poly, data.XtX[i, j])
end
cov(poly::Poly, X1::AbstractMatrix, X2::AbstractMatrix, data::LinIsoData) = _cov(poly, data.XtX)

get_params(poly::Poly) = [log(poly.c), log(poly.σ2) / 2]
get_param_names(poly::Poly) = [:lc, :lσ]
num_params(poly::Poly) = 2

function set_params!(poly::Poly, hyp::AbstractVector)
    length(hyp) == 2 || throw(ArgumentError("Polynomial function has two parameters"))
    poly.c = exp(hyp[1])
    poly.σ2 = exp(2 * hyp[2])
end

@inline dk_dlc(poly::Poly, xTy::Real) = poly.c*poly.deg*poly.σ2*(poly.c+xTy).^(poly.deg-1)
@inline dk_dlσ(poly::Poly, xTy::Real) = 2 * _cov(poly,xTy)
@inline function dKij_dθp(poly::Poly, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dlc(poly, dotij(X1,X2,i,j,dim))
    else
        return dk_dlσ(poly, dotij(X1, X2,i,j,dim))
    end
end
@inline function dKij_dθp(poly::Poly, X1::AbstractMatrix, X2::AbstractMatrix, data::LinIsoData, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dlc(poly, data.XtX[i,j])
    else
        return dk_dlσ(poly, data.XtX[i,j])
    end
end


text:
export DiagOp

"""
`DiagOp(domainType::Type, dim_in::Tuple, d::AbstractArray)`

`DiagOp(d::AbstractArray)`

Creates a `LinearOperator` which, when multiplied with an array `x`, returns the elementwise product `d.*x`.

```julia
julia> D = DiagOp(Float64, (2, 2,), [1. 2.; 3. 4.])
╲  ℝ^(2, 2) -> ℝ^(2, 2)

julia> D*ones(2,2)
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

```

"""
struct DiagOp{N, D, C, T <: Union{AbstractArray{C,N}, Number}} <: LinearOperator
	dim_in::NTuple{N,Int}
	d::T
end

# Constructors

###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int})
function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: AbstractArray}
	size(d) != DomainDim && error("dimension of d must coincide with DomainDim")
    C = eltype(d) <: Complex ? complex(DomainType) : DomainType
    DiagOp{N, DomainType, C, T}(DomainDim, d)
end

###standard constructor with Scalar
function DiagOp(DomainType::Type, DomainDim::NTuple{N,Int}, d::T) where {N, T <: Number}
    C = eltype(d) <: Complex ? Complex{DomainType} : DomainType
    DiagOp{N, DomainType, C, T}(DomainDim, d)
end

# other constructors
DiagOp(d::A) where {A <: AbstractArray} = DiagOp(eltype(d), size(d), d)
DiagOp(DomainDim::NTuple{N,Int}, d::A) where {N, A <: Number} = DiagOp(Float64, DomainDim, d)

# Mappings

function mul!(y::AbstractArray{C,N}, L::DiagOp{N, D, C, T}, b::AbstractArray{D,N}) where {N, D, C, T}
	y .= L.d.*b
end

function mul!(y::AbstractArray{D,N}, L::AdjointOperator{DiagOp{N, D, C, T}}, b::AbstractArray{C,N}) where {N, D, C, T}
	y .= conj.(L.A.d).*b
end

function mul!(y::AbstractArray{D,N}, L::AdjointOperator{DiagOp{N, D, C, T}}, b::AbstractArray{C,N}) where {N, D <: Real, C <: Complex, T}
    y .= real.(conj.(L.A.d).*b)
end

# Transformations (we'll see about this)
# inv(L::DiagOp) = DiagOp(L.domainType, L.dim_in, (L.d).^(-1))

# Properties

diag(L::DiagOp) = L.d
diag_AAc(L::DiagOp) = L.d.*conj.(L.d)
diag_AcA(L::DiagOp) = conj.(L.d).*L.d

domainType(L::DiagOp{N, D, C, T}) where {N, D, C, T} = D
codomainType(L::DiagOp{N, D, C, T}) where {N, D, C, T} = C

size(L::DiagOp) = (L.dim_in, L.dim_in)

fun_name(L::DiagOp) = "╲"

is_diagonal(L::DiagOp) = true

# TODO: probably the following allows for too-close-to-singular matrices
is_invertible(L::DiagOp) = all(L.d .!= 0.)
is_full_row_rank(L::DiagOp) = is_invertible(L)
is_full_column_rank(L::DiagOp) = is_invertible(L)


text:
using Documenter, DemoCards
using ImageFiltering
using OffsetArrays
# loading Plots before generating demos (might be helpful to avoid segfaults)
using Plots
ENV["GKSwstype"] = "nul" # needed for the GR backend on headless servers
gr()

demos, demos_cb, demos_assets = makedemos("demos")

assets = []
isnothing(demos_assets) || (push!(assets, demos_assets))
format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = assets)

makedocs(
    modules  = [ImageFiltering, OffsetArrays, Kernel, KernelFactors, ImageFiltering.MapWindow],
    format   = format,
    sitename = "ImageFiltering",
    pages    = [
        "index.md",
        demos,
        "Function reference" => "function_reference.md"
    ],
    # Note(johnnychen94): doctest is moved as part of unit test for two reasons:
    # - we want to run doctest for various Julia versions on various platforms.
    # - `doctest=true` for `makedocs` will trigger not only doctest of ImageFiltering but also
    #   that of other packages listed in `modules`. It usually throws warnings that don't
    #   belong to the scope of this package.
    doctest = false
)

demos_cb()

deploydocs(
    repo   = "github.com/JuliaImages/ImageFiltering.jl.git",
    # NOTE: remove the comment if the PR requires a push preview
    # push_preview = true
)


text:
hassubset(N::Type{<:Name}, type::Symbol) = hassubset(N, Val(type))
function hassubset(::Type{<:Name}, type::Val)
    throw(ArgumentError("$(type) subset not supported. Possible choices: :train, :valid, :test"))
end
hassubset(N::Type{<:Name}, ::Val{:train}) = ninstances(N)[1] > 0
hassubset(N::Type{<:Name}, ::Val{:valid}) = ninstances(N)[2] > 0
hassubset(N::Type{<:Name}, ::Val{:test}) = ninstances(N)[3] > 0

"""
    datapath(N::Type{<:Name}, type::Symbol)

Returns the path to the file in which the `type` subset of dataset `N` is stored and download the dataset if it not installed yet. Valid subset types are: `:train`, `:valid`, `:test`.
"""
function datapath(N::Type{<:Name}, type::Symbol)
    if !hassubset(N, type)
        throw(ArgumentError("$(type) subset not provided for $(nameof(N)) dataset"))
    end
    ext = saveformat(format(N))
    return joinpath(datadir(N), "data_$(type).$(ext)")
end

datadir(N::Type{<:Name}) = @datadep_str(datadepname(N))

function saveraw(N::Type{<:Name}, path, type, data; clearpath = true)
    saveraw(format(N), datapath(N, type), data)
    clearpath && rm(path)
    return
end

loadraw(N::Type{<:Name}, type) = loadraw(format(N), datapath(N, type))

hasmeta(N::Type{<:Name}) = isfile(joinpath(datadir(N), "meta.bson"))
savemeta(N::Type{<:Name}, meta) = BSON.bson(joinpath(datadir(N), "meta.bson"), meta)
loadmeta(N::Type{<:Name}) = BSON.load(joinpath(datadir(N), "meta.bson"))

# MLDatasets
loadraw_mldatasets(modul_name, type::Symbol) = loadraw_mldatasets(modul_name, Val(type))
function loadraw_mldatasets(::Type{<:Name}, type::Val)
    throw(ArgumentError("$(type) subset not supported. Possible choices: :train, :test"))
end
loadraw_mldatasets(dataset, ::Val{:train}, T = Float32) = dataset.traindata(T)
loadraw_mldatasets(dataset, ::Val{:test}, T = Float32) = dataset.testdata(T)

# unzip
function unzip(path; clearpath = true, rename = identity, skip = (f) -> false)
    r = ZipFile.Reader(path)
    for f in filter(!skip, r.files)
        write(joinpath(dirname(path), rename(f.name)), read(f, String))
    end
    clearpath && rm(path)
    return
end

"""
    remove(N::Type{<:Name})
Removes given dataset `N`, if installed.
"""
function remove(N::Type{<:Name})
    path = DataDeps.try_determine_load_path(datadepname(N), @__DIR__)
    if !isnothing(path)
        rm(path; recursive = true)
        printstyled(" ✖ $(nameof(N)) dataset removed \n"; color = :red)
    else
        printstyled(" ✖ $(nameof(N)) dataset not installed \n"; color = :red)
    end
    return
end

"""
    removeall()
Removes all downloaded datasets.
"""
function removeall()
    if DataDeps.input_bool("Do you want to remove all downloaded datasets?")
        remove.(name_subtypes())
    end
    return
end

"""
    listdatasets()
Prints the names of all available datasets and highlights datasets that have already been downloaded. Datasets are sorted into groups based on their problem type and format type.
"""
function listdatasets()
    Ns = name_subtypes()
    Ts = task_subtypes()
    Fs = format_subtypes()

    for T in Ts
        printstyled("$(nameof(T)): \n"; color = :blue, bold = true)
        for F in Fs
            datasets = [N for N in Ns if task(N) <: T && format(N) <: F]
            isempty(datasets) && continue
            printstyled("  $(nameof(F)): \n"; color = :yellow, bold = true)
            for N in datasets
                path = DataDeps.try_determine_load_path(datadepname(N), @__DIR__)
                if isnothing(path)
                    printstyled("    ✖ $(nameof(N)) \n"; color = :red)
                else
                    printstyled("    ✔ $(nameof(N)) \n"; color = :green, bold = true)
                end
            end
        end
    end
    return
end

# HDF5 array utilities
struct HDF5Array{T,N} <: AbstractArray{T,N}
    data
    inds

    function HDF5Array(data, inds = 1:size(data, ndims(data)))
        return new{eltype(data),ndims(data)}(data, inds)
    end
end

Base.size(A::HDF5Array) = (size(A.data)[1:end-1]..., length(A.inds))

Base.summary(A::HDF5Array{T,1}) where {T} = "$(length(A))-element HDF5Vector{$T}"
Base.summary(A::HDF5Array{T,2}) where {T} = "$(join(size(A), "x")) HDF5Matrix{$T}"
Base.summary(A::HDF5Array{T,N}) where {T,N} = "$(join(size(A), "x")) HDF5Array{$T}"

Base.show(io::IO, A::HDF5Array) = print(io, summary(A))
Base.show(io::IO, ::MIME"text/plain", A::HDF5Array) = show(io, A)

function Base.getindex(A::HDF5Array, I...)
    indices = Base.to_indices(A, I)
    return getindex(A.data, indices[1:end-1]..., A.inds[indices[end]])
end

function Base.setindex!(A::HDF5Array, v, I...)
    indices = Base.to_indices(A, I)
    return setindex!(A.data, v, indices[1:end-1]..., A.inds[indices[end]])
end

function Base.selectdim(A::HDF5Array{T,N}, d::Int, inds) where {T,N}
    if d != N
        throw(ArgumentError("only last dimension ($N) can be selected"))
    end
    return HDF5Array(A.data, A.inds[inds])
end


text:
## index_set.jl : unit testing for index_set.jl
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

@testset "IndexSet                     " begin

for (T, n, i) in Base.Iterators.zip(["SL" "ML"], [1 6], [5, 0])
    @eval begin
        t = $(Symbol(T))()
        @test isa(t, $(Symbol(T)))
        @test isa(t, $(Symbol(T)){1})
        idxset = collect(get_index_set(t, 5))
        @test length(idxset) == $n
        @test idxset[1] == Level($i)
        @test idxset[end] == Level(5)
        @test_throws MethodError $(Symbol(T))(2)
        @sprintf "%s" t
    end
end

for (T, n3, n2) in Base.Iterators.zip(["FT" "TD" "HC" "ZC"], [512 120 38 98], [28 17 11 15])
    @eval begin
        #
        # 3d index set
        #
        t3 = $(Symbol(T))(3)
        @test isa(t3, $(Symbol(T)){3})
        idxset3 = collect(get_index_set(t3, 7))
        @test length(idxset3) == $n3
        @test idxset3[1] == Index(0, 0, 0)
        @test Index(1, 1, 1) ∈ idxset3
        @test Index(8, 0, 0) ∉ idxset3
        @test Index(0, 8, 0) ∉ idxset3
        @test Index(0, 0, 8) ∉ idxset3
        #
        # weighted d2 index set
        #
        t2 = $(Symbol(T))(1/1, 3/5)
        @test isa(t2, $(Symbol(T)){2})
        idxset2 = collect(get_index_set(t2, 6))
        @test length(idxset2) == $n2
        @test idxset2[1] == Index(0, 0)
        @test Index(1, 1) ∈ idxset2
        @test Index(0, 4) ∉ idxset2
        @test Index(7, 0) ∉ idxset2
        @test_throws ArgumentError $(Symbol(T))(1)
        @test_throws ArgumentError $(Symbol(T))(1, 0)
        @sprintf "%s" t2
    end
end

ad = AD(2)
@test ad isa AD
@test ad isa AD{2}
@test_throws ArgumentError get_index_set(ad, 4)
@test_throws ArgumentError AD(1)
@test_throws MethodError AD()
@sprintf "%s" ad

u = U(1)
@test u isa U
@test u isa U{1}
@test_throws ArgumentError get_index_set(u, 1)
@test_throws MethodError U()
@test_throws MethodError U(1, 2)
@sprintf "%s" u

m = U(2)
@test m isa U{2}
@test_throws ArgumentError get_index_set(m, 4)
@sprintf "%s" m

@test length(collect(get_index_set(TD(2), 0))) == 1
@test isempty(get_index_set(TD(2), -1))

end


text:
module TableScraper

using Gumbo: parsehtml
using Cascadia: @sel_str, nodeText
import HTTP;

export scrape_tables

include("Tables.jl")

"""
    scrape_tables(url)
    scrape_tables(url, cell_transform[=nodeText])
    scrape_tables(url, cell_transform[=nodeText], header_transform[=nodeText])

This function will scrape `url` for any WELL-FORMED tables wrapped in `<table>` tags and return
them in a vector.

# Arguments

    - `url`: The URL to look for tables
    - `cell_transform`: By default, each of the table cells wrapped in `<td>` is transformed by
        the callable (i.e. `Function` or type definition) `cell_transform`. The default
        `cell_transform` is `Cascadia.nodeText` which extracts the node's text. You may wish to use
        `identity` to extract just the cell as a `Gumbo.HTMLNode` type for more advanced processing,
        e.g. `scrape_tables(url, identity)`
    - `header_transform`: By default, each of the table header wrapped in `<th>` is transformed by
        the callable (i.e. `Function` or type definition) `header_transform`. The default
        `header_transform` is `Cascadia.nodeText` which extracts the node's text. You may wish to
        use `identity` to extract just the cell as a `Gumbo.HTMLNode` type for more advanced
        processing, e.g. `scrape_tables(url, identity, identity)`

# Return

    `Vector{TableScraper.Table}`

The `TableScraper.Table` is a Tables.jl-compatible row-accessible type. So you can convert it to
another Tables.jl compatible type if you wish e.g. `DataFrame.(scrape_tables(url))` will return a
vector of `DataFrame`s
"""
function scrape_tables(url, cell_transform=nodeText, header_transform=nodeText)::Vector{Table}
    result_tables = []

    response::HTTP.Messages.Response =
        try
            response = HTTP.get(url)
        catch e
            println("Error when attempting to get $url. Make you are connected to the internet and the URL is accessible")
            raise(e)
        end

    # the body is the html content
    parsed_html = parsehtml(String(response.body))

    # look for tables in the parsed_html
    tables_elems = eachmatch(sel"table", parsed_html.root)

    n_tables = length(tables_elems)

    headers = [[] for _ in 1:n_tables]

    for (header, table_elem) in zip(headers, tables_elems)
        for header1 in eachmatch(sel"tr th", table_elem)
            # check the header span
            # you are on your won if you don't use nodeText
            if (nodeText == header_transform) & haskey(header1.attributes, "colspan")
                colspan = parse(Int, header1.attributes["colspan"])
                for i in 1:colspan
                    push!(header, nodeText(header1)*"$i")
                end
            else
                push!(header, header_transform(header1))
            end
        end
    end

    result_tables = [[] for _ in 1:n_tables]

    # parse the tables
    # Fill dataframe with rows from the table
    for (results_arr, table_elem) in zip(result_tables, tables_elems)
        for row in eachmatch(sel"tbody tr", table_elem)
            tds = cell_transform.(eachmatch(sel"td", row))
            if length(tds) > 0 # likely to be header column so don't add
                push!(results_arr, tds)
            end
        end
    end

    [Table(rows, header) for (rows, header) in zip(result_tables, headers)]
end

end


text:
# ---
# title: 1405. Longest Happy String
# id: problem1405
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/longest-happy-string/description/>
# hidden: true
# ---
#
# A string is called _happy_ if it does not have any of the strings `'aaa'`,
# `'bbb'` or `'ccc'` as a substring.
#
# Given three integers `a`, `b` and `c`, return **any** string `s`, which
# satisfies following conditions:
#
#   * `s` is _happy  _and longest possible.
#   * `s` contains **at most** `a` occurrences of the letter `'a'`, **at most** `b` occurrences of the letter `'b'` and **at most** `c` occurrences of the letter `'c'`.
#   * `s `will only contain `'a'`, `'b'` and `'c'` letters.
#
# If there is no such string `s` return the empty string `""`.
#
#
#
# **Example 1:**
#
#
#
#     Input: a = 1, b = 1, c = 7
#     Output: "ccaccbcc"
#     Explanation: "ccbccacc" would also be a correct answer.
#
#
# **Example 2:**
#
#
#
#     Input: a = 2, b = 2, c = 1
#     Output: "aabbc"
#
#
# **Example 3:**
#
#
#
#     Input: a = 7, b = 1, c = 0
#     Output: "aabaa"
#     Explanation: It's the only correct answer in this case.
#
#
#
#
# **Constraints:**
#
#   * `0 <= a, b, c <= 100`
#   * `a + b + c > 0`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end