text:
#### Misc. utils

mxdict(args...) = Dict{String, Any}(args...)
mxtuple(args...) = Any[args...]
mxempty() = zeros(Float64, 0, 0)

#### Temporary workspace

function reset_active_project(f)
    curr_proj = Base.active_project()
    try
        return f()
    finally
        Pkg.activate(curr_proj; io = devnull)
    end
end

function jlcall_workspace()
    workspace = ENV["MATDAEMON_WORKSPACE"]
    isfile(joinpath(workspace, "Project.toml")) && return workspace

    reset_active_project() do
        Pkg.activate(workspace; io = devnull)
        Pkg.develop(PackageSpec(path = realpath(joinpath(@__DIR__, ".."))); io = devnull)
        return workspace
    end
end

function jlcall_test_project()
    test_proj = joinpath(@__DIR__, "TestProject")
    isfile(joinpath(test_proj, "Manifest.toml")) && return test_proj

    reset_active_project() do
        Pkg.activate(test_proj; io = devnull)
        Pkg.instantiate(; io = devnull)
        return test_proj
    end
end

#### Recursive, typed equality testing

recurse_is_equal(eq) = (x, y) -> recurse_is_equal(eq, x, y)
recurse_is_equal(eq, x, y) = eq(x, y) #fallback
recurse_is_equal(eq, x::Tuple, y::Tuple) = eq(x, y) && all(recurse_is_equal(eq, x[i], y[i]) for i in 1:length(x))
recurse_is_equal(eq, x::KeyValueContainer, y::KeyValueContainer) = eq(x, y) && all(recurse_is_equal(eq, x[k], y[k]) for k in keys(x))
recurse_is_equal(eq, x::JLCallOptions, y::JLCallOptions) = all(recurse_is_equal(eq, getproperty(x, k), getproperty(y, k)) for k in fieldnames(JLCallOptions))

typed_is_equal(eq) = (x, y) -> typed_is_equal(eq, x, y)
typed_is_equal(eq, x, y) = eq(x, y) && typeof(x) == typeof(y)

verbose_typed_is_equal(eq) = (x, y) -> verbose_typed_is_equal(eq, x, y)
verbose_typed_is_equal(eq, x, y) = typed_is_equal(eq, x, y) || pprint_compare(x, y)

is_eq(x, y) = recurse_is_equal(verbose_typed_is_equal(isequal), x, y)
is_eqq(x, y) = recurse_is_equal(verbose_typed_is_equal(===), x, y)

#### Custom types

# Struct without custom matlabify; fields will not be recursively matlabified
struct A
    x
    y
end

recurse_is_equal(eq, a1::A, a2::A) = recurse_is_equal(eq, (a1.x, a1.y), (a2.x, a2.y))

# Struct with custom matlabify; fields will be recursively matlabified
struct B
    x
    y
end

recurse_is_equal(eq, b1::B, b2::B) = recurse_is_equal(eq, (b1.x, b1.y), (b2.x, b2.y))

MATDaemon.matlabify(b::B) = mxdict("x" => matlabify(b.x), "y" => matlabify(b.y))

#### Building deeply nested types for testing

function deeply_nested_pairs(; roundtrip = false)
    # List of pairs of primitive base cases
    ps = Any[
        nothing                => mxempty(),
        missing                => mxempty(),
        1                      => 1,
        2.0                    => 2.0,
        "string"               => "string",
        :symbol                => "symbol",
        [1.0, 2.0]             => [1.0, 2.0],
        [3.0 4.0]              => [3.0 4.0],
        ones(Float32, 1, 1, 2) => ones(Float32, 1, 1, 2),
        [3.0]                  => roundtrip ? 3.0 : [3.0],
        ones(1, 1, 1)          => roundtrip ? 1.0 : ones(1, 1, 1),
        zeros(1, 1, 2, 1, 1)   => roundtrip ? zeros(1, 1, 2) : zeros(1, 1, 2, 1, 1),
        trues(1, 2)            => roundtrip ? fill(true, 1, 2) : trues(1, 2),
    ]
    nprimitives = length(ps)

    # Iteratively build nested containers from primitives and push them to the list
    for _ in 1:nprimitives
        jl1, mx1 = ps[rand(1:end)]
        jl2, mx2 = ps[rand(1:end)]
        jl3, mx3 = ps[rand(1:end)]
        mxdict_12 = mxdict("x" => mx1, "y" => mx2)

        # tuple -> mxtuple
        push!(ps, (jl1, jl2, jl3) => mxtuple(mx1, mx2, mx3))

        # various associative containers => mxdict
        push!(ps, (; x = jl1, y = jl2) => mxdict_12) # named tuple
        push!(ps, pairs((; x = jl1, y = jl2)) => mxdict_12) # pairs iterator
        push!(ps, Dict{String, Any}("x" => jl1, "y" => jl2) => mxdict_12) # dict with string keys
        push!(ps, Dict{Symbol, Any}(:x => jl1, :y => jl2) => mxdict_12) # dict with symbol keys

        # custom types
        if !roundtrip
            push!(ps, A(jl1, jl2) => A(jl1, jl2)) # default `matlabify` falls back to identity
        end
        push!(ps, B(jl1, jl2) => mxdict_12) # custom `matlabify` recursively converts fields
    end

    return map(deepcopy, ps)
end

#### Pretty printing for deeply nested arguments

struct PrettyCompare
    arg1
    arg2
end

function pprint_compare(x, y)
    GarishPrint.pprint(PrettyCompare(x, y); compact = false, limit = false)
    println("")
    return false
end

#### Convenience method for calling and testing `MATDaemon.jlcall`

function wrap_jlcall(f, f_args, f_kwargs, f_output; kwargs...)
    opts = JLCallOptions(;
        f = f,
        workspace = jlcall_workspace(),
        debug = false,
        kwargs...,
    )
    optsfile = joinpath(opts.workspace, MATDaemon.JL_OPTIONS)

    try
        f_args_dict = Dict(:args => f_args, :kwargs => f_kwargs)
        MAT.matwrite(opts.infile, matlabify(f_args_dict))
        MAT.matwrite(optsfile, matlabify(opts))

        @test isfile(opts.infile)
        @test isfile(optsfile)
        @test is_eq(opts, MATDaemon.load_options(opts.workspace))

        Main.include(MATDaemon.jlcall_script())

        @test isfile(opts.outfile)
        @test is_eq(f_output, MAT.matread(opts.outfile)["output"])
    finally
        rm(optsfile; force = true)
        rm(opts.infile; force = true)
        rm(opts.outfile; force = true)
    end
end

#### Convenience method for calling `jlcall.m` via MATLAB.jl

function mx_wrap_jlcall(
        nargout::Int,
        f::String = "(args...; kwargs...) -> nothing",
        f_args::Tuple = (),
        f_kwargs::NamedTuple = NamedTuple();
        kwargs...
    )

    opts = JLCallOptions(;
        f = f,
        workspace = jlcall_workspace(),
        debug = false,
        gc = true,
        port = 5678,
        kwargs...
    )
    optsfile = joinpath(opts.workspace, MATDaemon.JL_OPTIONS)

    mxargs = Any[
        opts.f,
        matlabify(f_args),
        matlabify(f_kwargs),
    ]

    for k in fieldnames(JLCallOptions)
        k === :f && continue
        push!(mxargs, string(k))
        push!(mxargs, getproperty(opts, k))
    end

    f_output = mxcall(:jlcall, nargout, mxargs...)

    @test xor(isfile(optsfile), opts.gc)
    @test xor(isfile(opts.infile), opts.gc)
    @test xor(isfile(opts.outfile), opts.gc)

    return f_output
end


text:
using FunctionalCollections

include("./bench.jl")

println("\nPersistentVector")
println("================\n")

const rands = rand(1:1000000, 100000)

function appending(::Type{PersistentVector})
    function ()
        v = PersistentVector{Int}()
        for i=1:250000
            v = append(v, i)
        end
    end
end
function appending(::Type{Array})
    function ()
        a = Int[]
        for i=1:250000
            a = push!(a, i)
        end
    end
end

@bench "Appending" 20 [appending(PersistentVector),
                       appending(Array)]

vec(r::Range1) = PersistentVector([r])

function iterating(::Type{PersistentVector})
    pv = vec(1:500000)
    function ()
        sum = 0
        for el in pv
            sum += el
        end
    end
end
function iterating(::Type{Array})
    arr = [1:500000]
    function ()
        sum = 0
        for el in arr
            sum += el
        end
    end
end

@bench "Iterating" 20 [iterating(PersistentVector),
                       iterating(Array)]

function indexing(v::PersistentVector{T}) where T
    function ()
        for idx in rands
            v[idx]
        end
    end
end
function indexing(a::Array{T}) where T
    function ()
        for idx in rands
            a[idx]
        end
    end
end

@bench "Indexing" 20 [indexing(vec(1:1000000)), indexing(Array(Int, 1000000))]

function popping(::Type{PersistentVector})
    v = vec(1:100000)
    function ()
        v2 = v
        for _ in 1:length(v2)
            v2 = pop(v2)
        end
    end
end
function popping(::Type{Array})
    function ()
        a = Array(Int, 100000)
        for _ in 1:length(a)
            pop!(a)
        end
    end
end

@bench "Popping" 20 [popping(PersistentVector), popping(Array)]

function updating(v::PersistentVector{Int})
    function ()
        for idx in rands
            assoc(v, idx, 1)
        end
    end
end
function updating(a::Array{Int})
    function ()
        for idx in rands
            a[idx] = 1
        end
    end
end

@bench "Updating" 20 [updating(vec(1:1000000)), updating(Array(Int, 1000000))]

ext: jl


text:
if length(ARGS) ≠ 1
    error("""
        Only one argument is required:
            `\"movie_title\"`
    """)
end

title = ARGS[1]

i =


text:
using ParametricMachinesDemos
using BenchmarkTools, ChainRulesCore, Flux, Zygote, CUDA, Random

# Set seed on CPU and GPU

Random.seed!(1234)
CUDA.seed!(1234)

# Benchmarking

timeratio(b1, b2) = minimum(b1.times) / minimum(b2.times)

function benchmark((y₀, z₀, ȳ, z̄, machine), device)
    σ, W = machine.σ, machine.W
    forward, backward = ParametricMachinesDemos.filtrations(machine, y₀)

    function forward_pass(y₀, z₀)
        y, z = copy(y₀), copy(z₀)
        W̃ = ParametricMachinesDemos.clean!(copy(W), forward)
        ParametricMachinesDemos.solve!(y, z, W̃, σ, forward)
        return y, z
    end

    # Benchmark forward step
    benchmark_forward = if device == "CPU"
        @benchmark $forward_pass($y₀, $z₀)
    elseif device == "GPU"
        @benchmark CUDA.@sync $forward_pass($y₀, $z₀)
    end

    y, z = forward_pass(y₀, z₀)

    function backward_pass(ȳ, z̄)
        Dσ = ParametricMachinesDemos.derivative.(z .- z₀, σ, y)
        W̃′ = ParametricMachinesDemos.clean!(ParametricMachinesDemos.flip(W), backward)
        ȳ₀, z̄₀ = copy(ȳ), copy(z̄)
        ParametricMachinesDemos.solve!(z̄₀, ȳ₀, W̃′, Dσ, backward)
        return ȳ₀, z̄₀
    end

    # Benchmark backward step
    benchmark_backward = if device == "CPU"
        @benchmark $backward_pass($ȳ, $z̄)
    elseif device == "GPU"
        @benchmark CUDA.@sync $backward_pass($ȳ, $z̄)
    end
    @show timeratio(benchmark_backward, benchmark_forward)
    return (forward=benchmark_forward, backward=benchmark_backward)
end

df = @NamedTuple{device::String, problemsize::String, machine::String, forward::Float64, backward::Float64}[]
devices = (CPU=f64, GPU=gpu)

for problemsize in ["small", "medium"]
    dimensions = problemsize == "small" ? [2, 2, 2, 2, 2] : [32, 32, 32, 32, 32]
    minibatch = problemsize == "small" ? 2 : 32
    y₀ = rand(Float32, sum(dimensions), minibatch)
    z₀ = rand(Float32, sum(dimensions), minibatch)
    ȳ = rand(Float32, sum(dimensions), minibatch)
    z̄ = rand(Float32, sum(dimensions), minibatch)
    machine = DenseMachine(dimensions, tanh)

    for device in ["CPU", "GPU"]
        to_device = devices[Symbol(device)]
        benchmarks = benchmark(map(to_device, (y₀, z₀, ȳ, z̄, machine)), device)
        push!(
            df,
            (; device, problemsize, machine="dense",
            forward=minimum(benchmarks.forward.times),
            backward=minimum(benchmarks.backward.times))
        )
    end
end

# Convolutional machine

for problemsize in ["small", "medium"]
    dimensions = problemsize == "small" ? [2, 2, 2, 2, 2] : [32, 32, 32, 32, 32]
    minibatch = problemsize == "small" ? 2 : 32
    duration = 16
    pad = (3, 3)
    y₀ = rand(Float32, duration, sum(dimensions), minibatch)
    z₀ = rand(Float32, duration, sum(dimensions), minibatch)
    ȳ = rand(Float32, duration, sum(dimensions), minibatch)
    z̄ = rand(Float32, duration, sum(dimensions), minibatch)
    machine = ConvMachine(dimensions, tanh; pad)

    for device in ["CPU", "GPU"]
        to_device = devices[Symbol(device)]
        benchmarks = benchmark(map(to_device, (y₀, z₀, ȳ, z̄, machine)), device)
        push!(
            df,
            (; device, problemsize, machine="convolution",
            forward=minimum(benchmarks.forward.times),
            backward=minimum(benchmarks.backward.times))
        )
    end
end

# Recurrent machine

for problemsize in ["small", "medium"]
    dimensions = problemsize == "small" ? [2, 2, 2, 2, 2] : [32, 32, 32, 32, 32]
    minibatch = problemsize == "small" ? 2 : 32
    duration = 16
    pad = 6
    timeblock = 5
    y₀ = rand(Float32, duration, sum(dimensions), minibatch)
    z₀ = rand(Float32, duration, sum(dimensions), minibatch)
    ȳ = rand(Float32, duration, sum(dimensions), minibatch)
    z̄ = rand(Float32, duration, sum(dimensions), minibatch)
    machine = RecurMachine(dimensions, tanh; pad, timeblock)

    for device in ["CPU", "GPU"]
        to_device = devices[Symbol(device)]
        benchmarks = benchmark(map(to_device, (y₀, z₀, ȳ, z̄, machine)), device)
        push!(
            df,
            (; device, problemsize, machine="recurrent",
            forward=minimum(benchmarks.forward.times),
            backward=minimum(benchmarks.backward.times))
        )
    end
end

# Log to file

using CSV

CSV.write(joinpath(@__DIR__, "benchmarks.csv"), df)


text:
using Pkg
tempdir = mktempdir()
Pkg.activate(tempdir)
Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
Pkg.add(["BenchmarkTools", "PkgBenchmark", "Random", "Distances", "LightGraphs"])
Pkg.resolve()

using NearestNeighborDescent
using BenchmarkTools

const SUITE = BenchmarkGroup()

SUITE["knn graphs"] = include("bench_knn_graphs.jl")
SUITE["descent"] = include("bench_descent.jl")


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Distributions
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity
using OceanTurbulenceParameterEstimation

#####
##### Parameters
#####

Nz = 32
Lz = 128
Qᵇ = 1e-8
Qᵘ = -1e-5
Δt = 20.0
f₀ = 1e-4
N² = 1e-5
stop_time = 4hours
save_interval = 1hour
experiment_name = "convective_adjustment"
data_path = experiment_name * ".jld2"
ensemble_size = 50
generate_observations = false

# "True" parameters to be estimated by calibration
convective_κz = 1.0
convective_νz = 0.9
background_κz = 1e-4
background_νz = 1e-5

θ★ = [convective_κz, background_κz, convective_νz, background_νz]

#####
##### Generate synthetic observations
#####

if generate_observations || !(isfile(data_path))
    grid = RegularRectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))
    closure = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz, background_κz, convective_νz, background_νz)
    coriolis = FPlane(f=f₀)

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        boundary_conditions = (; u=u_bcs, b=b_bcs),
                                        coriolis = coriolis,
                                        closure = closure)

    set!(model, b = (x, y, z) -> N² * z)

    simulation = Simulation(model; Δt, stop_time)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(save_interval),
                                                          prefix = experiment_name,
                                                          array_type = Array{Float64},
                                                          field_slicer = nothing,
                                                          force = true)

    run!(simulation)
end

#####
##### Load truth data as observations
#####

data_path = experiment_name * ".jld2"

observations = OneDimensionalTimeSeries(data_path, field_names=(:u, :b), normalize=ZScore)

#####
##### Set up ensemble model
#####

column_ensemble_size = ColumnEnsembleSize(Nz=Nz, ensemble=(ensemble_size, 1), Hz=1)
ensemble_grid = RegularRectilinearGrid(size=column_ensemble_size, z = (-Lz, 0), topology = (Flat, Flat, Bounded))
closure_ensemble = [ConvectiveAdjustmentVerticalDiffusivity(; convective_κz, background_κz) for i = 1:ensemble_grid.Nx, j = 1:ensemble_grid.Ny]
coriolis_ensemble = [FPlane(f=f₀) for i = 1:ensemble_grid.Nx, j = 1:ensemble_grid.Ny]

ensemble_b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))
ensemble_u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))

ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
                                             tracers = :b,
                                             buoyancy = BuoyancyTracer(),
                                             boundary_conditions = (; u=ensemble_u_bcs, b=ensemble_b_bcs),
                                             coriolis = coriolis_ensemble,
                                             closure = closure_ensemble)

set!(ensemble_model, b = (x, y, z) -> N² * z)

ensemble_simulation = Simulation(ensemble_model; Δt, stop_time)

pop!(ensemble_simulation.diagnostics, :nan_checker)
#####
##### Build free parameters
#####

priors = (
    convective_κz = lognormal_with_mean_std(1.2, 0.3),
    background_κz = lognormal_with_mean_std(1e-4, 1e-4),
    convective_νz = lognormal_with_mean_std(1.2, 0.3),
    background_νz = lognormal_with_mean_std(1e-5, 1e-5),
)

# priors = (
#     convective_κz = LogNormal(0, 0.7),
#     background_κz = LogNormal(0, 0.7),
#     convective_νz = LogNormal(0, 0.7),
#     background_νz = LogNormal(0, 0.7)
# )

free_parameters = FreeParameters(priors)

#####
##### Build the Inverse Problem
#####

calibration = InverseProblem(observations, ensemble_simulation, free_parameters)

# forward_map(calibration, [θ★ for _ in 1:ensemble_size])
x = forward_map(calibration, [θ★ for _ in 1:ensemble_size])
y = observation_map(calibration)

# Assert that G(θ*) ≈ y
@show x[:, 1:1] == y

using Plots, LinearAlgebra
# p = plot(collect(1:length(x)), [x...], label="forward_map")
# plot!(collect(1:length(y)), [y...], label="observation_map")
# savefig(p, "obs_vs_pred.png")
# display(p)

# Assert that G(θ*) ≈ y
@show forward_map(calibration, [θ★ for _ in 1:ensemble_size]) == observation_map(calibration)

iterations = 10
eki = EnsembleKalmanInversion(calibration; noise_covariance=1e-2)
params, mean_vars, mean_us = iterate!(eki; iterations = iterations)

@show params
y = eki.mapped_observations
a = [norm(forward_map(calibration, [mean_us[i] for _ in 1:ensemble_size])[:,1] - y) for i in 1:iterations]
plot(collect(1:iterations), a)


text:
using Base: Symbol
import LinearAlgebra

# We wanna implement some special stuff which might be bad for Base.Expr
struct MyExpr
    e::Expr

    function MyExpr(e::Expr, should_cleanup::Bool = true)
        if should_cleanup
            new(cleanup(e))
        else
            new(flatten_once(e))
        end
    end
end
MyExpr(args...; cleanup = true) = MyExpr(Expr(args...), cleanup)
Base.getproperty(e::MyExpr, field::Symbol) = getproperty(getfield(e, :e), field)
Base.:(==)(a::MyExpr, b::MyExpr) = getfield(a, :e) == getfield(b, :e)
expr2myexpr(x::Expr) = MyExpr(x)
expr2myexpr(x) = x
function Base.show(io::IO, me::MyExpr)
    print(io, "MyExpr(")
    e = getfield(me, :e)
    show(io, e.head)
    for arg in e.args
        print(io, ", ")
        show(io, arg)
    end
    print(io, ")")
end

function cleanup(e::Expr)
    if e.head == :call
        head = e.args[1]
        args = e.args[2:end]
    else
        head = e.head
        args = e.args
    end

    new_args = Any[]

    for child in args
        if child isa Expr
            child = cleanup(child)
        end
        if (child isa Expr || child isa MyExpr) && child.head == e.head
            append!(new_args, child.args)
        elseif child isa LineNumberNode; nothing;
        else
            push!(new_args, child)
        end
    end
    Expr(head, new_args...)
end

# I.e. a c or c^†
struct Operator{T1, T2}
    name::Symbol
    daggered::Bool
    index::T1
    time::T2
    function Operator(n, d, i, t)
        _i = expr2myexpr(i); _t = expr2myexpr(t)
        new{typeof(_i), typeof(_t)}(n, d, _i, _t)
    end
end

# constructors
create(idx, time=0; name = :c) = Operator(name, true, idx, time)
annihilate(idx, time=0; name = :c) = Operator(name, false, idx, time)
c(idx, time=0; name = :c) = annihilate(idx, time, name=name)
# cd(idx, time=0; name = :c) = create(idx, time, name=name)
function Base.:(==)(a::Operator, b::Operator)
    a.name == b.name && a.daggered == b.daggered && a.index == b.index && a.time == b.time
end

ExpectationValue(x::MyExpr) = MyExpr(:ExpectationValue, x)

struct ArrayElement
    name::Symbol
    indices::Vector{Any}
    ArrayElement(name::Symbol, idxs::Vector) = new(name, expr2myexpr.(idxs))
end
ArrayElement(name::Symbol, idxs...) = ArrayElement(name, collect(idxs))
function Base.:(==)(a::ArrayElement, b::ArrayElement)
    a.name == b.name && a.indices == b.indices
end

struct Sym
    name::Symbol
end
Base.:(==)(a::Sym, b::Sym) = a.name == b.name

# base extensions
LinearAlgebra.adjoint(o::Operator) = Operator(o.name, !o.daggered, o.index, o.time)
const OpOrExpr = Union{ArrayElement, Operator, MyExpr, Sym}
Base.:(-)(a::OpOrExpr) = MyExpr(:*, -1, a)
for T1 in (OpOrExpr, Number), T2 in (OpOrExpr, Number)
    T1 == Number && T2 == Number && continue
    Base.:(*)(a::T1, b::T2) = MyExpr(:*, a, b)
    Base.:(+)(a::T1, b::T2) = MyExpr(:+, a, b)
    Base.:(-)(a::T1, b::T2) = MyExpr(:+, a, -b)
    Base.:(*)(a::T1, b::T2, c, ds...) = *(MyExpr(:*, a, b), c, ds...)
    Base.:(+)(a::T1, b::T2, c, ds...) = +(MyExpr(:+, a, b), c, ds...)
end

########################################
### Printing
########################################

function Base.show(io::IO, o::Operator)
    print(io, o.name)
    o.daggered && print(io, "^†")
    print(io, "_{")
    _print_expr(io, o.index)
    print(io, "}(", o.time, ")")
end

print_expr(e::MyExpr) = print_expr(stdout, e)
function print_expr(io::IO, e::MyExpr)
    _print_expr(io, e)
    println(io)
end
function _print_expr(io, e::MyExpr)
    e.head == :+ && print(io, "(")
    e.head == :ExpectationValue && print(io, "⟨")
    for (i, child) in enumerate(e.args)
        i > 1 && print(io, " ", e.head, " ")
        _print_expr(io, child)
    end
    e.head == :ExpectationValue && print(io, "⟩")
    e.head == :+ && print(io, ")")
    nothing
end
function _print_expr(io, e::ArrayElement)
    print(io, e.name, "_{")
    for (i, idx) in enumerate(e.indices)
        i > 1 && print(io, ", ")
        _print_expr(io, idx)
    end
    print(io, "}")
end
_print_expr(io, e::Sym) = print(io, e.name)
_print_expr(io, e) = print(io, e)

########################################
### Expand brackets & flatten tree
########################################

# expand a * (b + c) -> a*b + a*c
# flatten nested operations
function expand(e::MyExpr)
    for i in eachindex(e.args)
        e.args[i] = expand(e.args[i])
    end
    new_expr = _expand(e)
    while new_expr != e
        e = new_expr
        new_expr = _expand(e)
    end
    return e
end
expand(e) = e

function _expand(e::MyExpr)
    # expand *(+(a, b), x) -> +(*(a, x), *(b, x))
    if e.head == :*
        i = 1
        while i <= length(e.args)
            child = e.args[i]
            if child isa MyExpr && child.head == :+
                products = map(child.args) do x
                    MyExpr(:*, e.args[1:i-1]..., x, e.args[i+1:end]..., cleanup=false)
                end
                new_expr = MyExpr(:+, products..., cleanup=false)
                # deal with other possible +(c, d) terms
                for j in eachindex(products)
                    new_expr.args[j] = _expand(new_expr.args[j])
                end
                # possibly flatten additions
                return _expand(new_expr)
            elseif child isa MyExpr && child.head == :*
                deleteat!(e.args, i)
                for x in child.args
                    insert!(e.args, i, x)
                    i += 1
                end
            else
                i += 1
            end
        end
    elseif e.head == :+
        i = 1
        while i < length(e.args)
            child = e.args[i]
            if child isa MyExpr && child.head == :+
                deleteat!(e.args, i)
                append!(e.args, child.args)
            else
                i += 1
            end
        end
    elseif e.head == :ExpectationValue
        if e.args[1].head == :+
            terms = [MyExpr(:ExpectationValue, x, cleanup=false) for x in e.args[1].args]
            return MyExpr(:+, terms..., cleanup=false)
        end
    end
    return e
end
_expand(e) = e

# flatten one level of nested operators
function flatten_once(e::Expr)
    new_args = Any[]
    for child in e.args
        if (child isa Expr || child isa MyExpr) && child.head == e.head
            append!(new_args, child.args)
        else
            push!(new_args, child)
        end
    end
    Expr(e.head, new_args...)
end

########################################
### Apply Wicks theorem
########################################

# Wicks theorem
# ⟨abcd⟩ ->  ⟨ab⟩⟨cd⟩ - ⟨ac⟩⟨bd⟩ + ⟨ad⟩⟨bc⟩ etc
function wicks(e)
    if e.head == :+
        products = wicks.(e.args)
        MyExpr(:+, products...)
    elseif e.head == :ExpectationValue
        wicks(e.args[1])
    else
        _wicks(e)
    end
end

function _wicks(e)
    @assert e.head == :*
    skip = filter(i -> !(e.args[i] isa Operator), eachindex(e.args))
    groups = generate_pairings(e.args, copy(skip))
    products = map(enumerate(groups)) do (i, g)
        evs = map(p -> MyExpr(:ExpectationValue, MyExpr(:*, p...)), g)
        if isodd(i)
            # index 1 3 5 ... have even numbers pair-permutations
            MyExpr(:*, e.args[skip]..., evs...)
        else
            MyExpr(:*, -1, e.args[skip]..., evs...)
        end
    end
    MyExpr(:+, products...)
end

function generate_pairings(v::Vector{T}, skip=Int[]) where T
    length(v) - length(skip) == 0 && return Vector{Any}[]
    @assert (length(v) - length(skip)) % 2 == 0 "$v"
    combos = Vector{Any}[]
    first_idx = findfirst(i -> !(i in skip), eachindex(v))
    push!(skip, first_idx)

    for i in first_idx+1:length(v)
        if !(i in skip)
            if length(skip)+1 < length(v)
                rest = generate_pairings(v, [skip; i])
                append!(combos, [vcat(Pair(v[first_idx], v[i]), r) for r in rest])
            else
                push!(combos, [Pair(v[first_idx], v[i])])
            end
        end
    end

    return combos
end

########################################
### convert expectation values to Greens functions
########################################

function to_greens(e::MyExpr)
    if e.head == :+
        new_args = Any[]
        for i in eachindex(e.args)
            converted = _to_greens(e.args[i])
            converted == 0.0 || push!(new_args, converted)
        end
        return MyExpr(:+, new_args...)
    else
        return _to_greens(e)
    end
end

function _to_greens(e::MyExpr)
    @assert e.head == :*

    outputs = []
    for (i, x) in enumerate(e.args)
        if x isa MyExpr && x.head == :ExpectationValue
            y = x.args[1]
            y.head == :* || error(y)
            length(y.args) == 2 || error(y)
            a, b = y.args
            if a.daggered && !b.daggered
                # c^† c - like
                delta = if a.time == b.time
                    ifelse(a.index == b.index, 1, ArrayElement(:I, a.index, b.index))
                elseif a.time isa Number && b.time isa Number
                    0
                else
                    ArrayElement(:I, a.index, b.index) * ArrayElement(:I, a.time, b.time)
                end
                x = delta - ArrayElement(Symbol(:G, b.time, a.time), b.index, a.index)
                push!(outputs, x)
            elseif !a.daggered && b.daggered
                # c c^†
                # GreensElement(a.indices[1], b.indices[2])
                push!(outputs, ArrayElement(Symbol(:G, a.time, b.time), a.index, b.index))
            else
                return 0.0
            end
        else
            push!(outputs, x)
        end
    end

    return *(outputs...)
end

########################################
### Replacements
########################################

struct PLACEHOLDER{T} end
PLACEHOLDER() = PLACEHOLDER{Any}()
struct REPLACEMENT
    x::Int64
end

replace(e, replacements::Pair...) = _replace1(e, replacements)
function _replace1(e, replacements::Tuple)
    e = _replace2(e, replacements)
    for (matchable, replacement) in replacements
        is_match, captures = match(e, matchable)
        if is_match
            return apply_replace(replacement, captures)
        end
    end
    return e
end
function _replace2(e::MyExpr, replacements)
    MyExpr(
        _replace1(e.head, replacements),
        [_replace1(child, replacements) for child in e.args]...
    )
end
function _replace2(ae::ArrayElement, replacements)
    ArrayElement(
        _replace1(ae.name, replacements),
        [_replace1(idx, replacements) for idx in ae.indices]...
    )
end
_replace2(e, replacements) = e

match(e, ::PLACEHOLDER{Any}) = true, Any[e]
match(e::T, ::PLACEHOLDER{T}) where {T} = true, Any[e]
match(a, b) = a == b, Any[]
function match(e::ArrayElement, m::ArrayElement)
    passthrough = Any[]
    if e.name != m.name || length(e.indices) != length(m.indices)
        return false, passthrough
    end
    for (a, b) in zip(e.indices, m.indices)
        is_match, x = match(a, b)
        is_match || return false, passthrough
        passthrough = vcat(passthrough, x)
    end
    return true, passthrough
end
function match(e::MyExpr, m::MyExpr)
    passthrough = Any[]
    if e.head != m.head || length(e.args) != length(m.args)
        return false, passthrough
    end
    for (a, b) in zip(e.args, m.args)
        is_match, x = match(a, b)
        is_match || return false, passthrough
        passthrough = vcat(passthrough, x)
    end
    return true, passthrough
end

apply_replace(r::REPLACEMENT, captures::Vector) = captures[r.x]
apply_replace(x, captures::Vector) = x
function apply_replace(r::MyExpr, captures::Vector)
    MyExpr(
        apply_replace(r.head, captures),
        (apply_replace(x, captures) for x in r.args)...
    )
end
function apply_replace(r::ArrayElement, captures::Vector)
    ArrayElement(
        apply_replace(r.name, captures),
        (apply_replace(x, captures) for x in r.indices)...
    )
end

function cleanup(e::MyExpr)
    replace(
        e,
        MyExpr(:+, PLACEHOLDER(), 0) => REPLACEMENT(1),
        MyExpr(:+, 0, PLACEHOLDER()) => REPLACEMENT(1),
        MyExpr(:*, PLACEHOLDER(), 0) => 0,
        MyExpr(:*, 0, PLACEHOLDER()) => 0,
        MyExpr(:*, PLACEHOLDER(), 1) => REPLACEMENT(1),
        MyExpr(:*, 1, PLACEHOLDER()) => REPLACEMENT(1),
    )
end

########################################
### To code
########################################

function to_expr(e::MyExpr)
    if e.head in (:+, :*)
        return Expr(:call, e.head, to_expr.(e.args)...)
    elseif e.head == :ExpectationValue
        error("Expectation values must be transformed.")
    else
        error("What is a $(e.head)?")
    end
end
to_expr(ae::ArrayElement) = Expr(:ref, ae.name, to_expr.(ae.indices)...)
to_expr(s::Symbol) = s
to_expr(x::Number) = x

function generate_kernel_function(
        e::Expr; mc_type = Any, model_type = Any,
        index_unpacking = default_index_unpacking(e),
        greens_unpacking = default_greens_unpacking(e),
    )

    func_head = Expr(
        :tuple,
        Expr(Symbol("::"), :mc, Symbol(mc_type)),
        Expr(Symbol("::"), :model, Symbol(model_type)),
        :packed_indices, :packed_greens,
    )
    func_body = Expr(
        :block,
        :(N = length(lattice(model))),
        index_unpacking, greens_unpacking,
        e
    )
    Expr(Symbol("->"), func_head, func_body)
end

function default_index_unpacking(e)
    indices = filter(idx -> idx != :N, collect_indices(e))
    if length(indices) == 2
        @assert indices == Set([:i, :j])
        return :(i, j = packed_indices)
    elseif length(indices) == 4
        @assert indices == Set([:i, :j, :k, :l])
        return :(i, j, k, l = packed_indices)
    else
        error("We can only deal with 2 or 4 indices atm. $(length(indices))")
    end
end

function default_greens_unpacking(e)
    Gs = filter(idx -> idx != :N, collect_greens(e))
    if length(Gs) == 1
        @assert Gs == Set([:G00])
        return :(G00 = packed_greens)
    elseif length(Gs) == 4
        @assert Gs == Set([:G00, :G0l, :Gl0, :Gll])
        return :(G00, G0l, Gl0, Gll = packed_greens)
    else
        error("We can only deal with 2 or 4 indices atm. $(length(Gs))")
    end
end

function collect_indices(e::Expr, indices=Set{Symbol}(), is_index=false)
    if e.head in (:ref, :call)
        for arg in e.args[2:end]
            collect_indices(arg, indices, true)
        end
    elseif e.head == :block
        for arg in e.args
            collect_indices(arg, indices, is_index)
        end
    else
        error("Failed to parse $(e.head)")
    end
    return indices
end
function collect_indices(s::Symbol, indices, is_index)
    is_index && push!(indices, s)
end
collect_indices(_, _, _) = nothing

function collect_greens(e::Expr, names=Set{Symbol}())
    if e.head == :ref
        @assert length(e.args) > 1
        push!(names, e.args[1])
    elseif e.head == :call
        for arg in e.args[2:end]
            collect_greens(arg, names)
        end
    elseif e.head == :block
        for arg in e.args
            collect_greens(arg, names)
        end
    else
        error("Failed to parse $(e.head)")
    end
    return filter(name -> startswith(string(name), 'G'), names)
end
collect_greens(_, _) = nothing

macro expand2formula(code)
    code |> expand |> wicks |> to_greens |> to_expr
end

macro expand2kernel(code)
    code |> expand |> wicks |> to_greens |> to_expr |> generate_kernel_function
end


text:
module Pkg2

struct PkgError <: Exception
    msg::AbstractString
    ex::Nullable{Exception}
end
PkgError(msg::AbstractString) = PkgError(msg, Nullable{Exception}())
function Base.showerror(io::IO, pkgerr::PkgError)
    print(io, pkgerr.msg)
    if !isnull(pkgerr.ex)
        pkgex = get(pkgerr.ex)
        if isa(pkgex, CompositeException)
            for cex in pkgex
                print(io, "\n=> ")
                showerror(io, cex)
            end
        else
            print(io, "\n")
            showerror(io, pkgex)
        end
    end
end

# DIR
const DIR_NAME = ".julia"
_pkgroot() = abspath(get(ENV,"JULIA_PKGDIR",joinpath(homedir(),DIR_NAME)))
isversioned(p::AbstractString) = ((x,y) = (VERSION.major, VERSION.minor); basename(p) == "v$x.$y")

function dir()
    b = _pkgroot()
    x, y = VERSION.major, VERSION.minor
    d = joinpath(b,"v$x.$y")
    if isdir(d) || !isdir(b) || !isdir(joinpath(b, "METADATA"))
        return d
    end
    return b
end
dir(pkg::AbstractString...) = normpath(dir(),pkg...)

include("types.jl")
include("reqs.jl")
include("query.jl")
# include("read.jl")
include("resolve.jl")

end


text:
#
include("mesh.jl")
include("melem.jl")

function buildmatrices_k0(msh::mesh, conn::Connectivity)
    #=
    Builds the matrix in full format (because I wanted to use the exact EV
    solver provided by Arpack).
    =#
    e2dof = conn.e2dof
    ndof  = conn.ndof
    # 1st: preallocate
    A = zeros(Float64, ndof, ndof)
    B = zeros(Float64, ndof, ndof)
    # 2nd: main loop
    for e=1:length(msh.elt)
        # @inbounds removes the check on the bounds of the array
        @inbounds verts  = msh.vtx[:, msh.elt[e]] # get the vertices
        # Choose if you want to compute the elementary matrix by projecting on
        # ∇P1 or on ∇P2 (for example for comparison with the RT0 :) )
        AE, BE = melemk0P1(verts)
        # AE, BE = melemk0P2(verts)
        nloc = length(e2dof[e])
        for i = 1:nloc
            @inbounds tmpI = e2dof[e][i]
            absI, sgnI = abs(tmpI), sign(tmpI)
            for j=1:nloc
                @inbounds tmpJ = e2dof[e][j]
                absJ, sgnJ = abs(tmpJ), sign(tmpJ)
                @inbounds A[absI, absJ] += sgnI*sgnJ*(AE[i,j] + BE[i,j])
                @inbounds B[absI, absJ] += sgnI*sgnJ*BE[i,j]
            end
        end
    end
    return A, B
end


text:
using SparseArrays, ExtendableSparse
using BenchmarkTools
using Test
import StatsBase: sample

##

function add_elements1!(mat::SparseMatrixCSC,ijs::Vector{CartesianIndex{2}})
  for ij in ijs
    mat[ij] = 123.456
  end
  return nothing
end

function add_elements2!(mat::SparseMatrixCSC,ijs::Vector{CartesianIndex{2}})
  mate = ExtendableSparseMatrix(mat)
  for ij in ijs
    i,j = Tuple(ij)
    setindex!(mate,123.456,i,j)
  end
  flush!(mate)
  copy!(mat,mate.cscmatrix)
  return nothing
end

##
# test that they do the same

N = 6000
mat = sprand(N,N,0.2)
nchange = 500
tochange = sample(CartesianIndices(mat),nchange;replace=false)
filter!(ij-> mat[ij]==0.0,tochange)
mat1 = deepcopy(mat)
mat2 = deepcopy(mat)

add_elements1!(mat1,tochange)
add_elements2!(mat2,tochange)

@test all( mat1 .== mat2)

##
@benchmark add_elements1!(_mat,elements) setup=(_mat=copy(mat);elements=tochange)
##
@benchmark add_elements2!(_mat,elements) setup=(_mat=copy(mat);elements=tochange)

##

function check_element2(mat::SparseMatrixCSC,ij::CartesianIndex{2})
  return ! iszero(mat[ij])
end

function check_element1(mat::SparseMatrixCSC,ij::CartesianIndex{2})
  i,j = Tuple(ij)
  return ! iszero(ExtendableSparse.findindex(mat,i,j))
end

@benchmark check_element1($mat,CartesianIndex(333,444))
@benchmark check_element2($mat,CartesianIndex(333,444))


text:
## This script gives a sensitivity analysis of models that fit to Nairobi data

using StatsPlots,Dates,JLD2,Statistics,Optim,Parameters,Distributions,DataFrames,CSV
using DynamicHMC,OrdinaryDiffEq,DiffEqCallbacks
using NamedArrays,TransformVariables,LogDensityProblems,LinearAlgebra,BlackBoxOptim,Random
using DynamicHMC.Diagnostics,MCMCDiagnostics,MCMCChains,ForwardDiff,LogExpFunctions
using Plots.PlotMeasures,LogExpFunctions
using CountTimeSeries
import KenyaCoVSD
##Load data for all Kenya
@load("data/linelist_data_with_pos_neg_20feb_to_27apr_c_age.jld2")
@load("data/cleaned_linelist20210521_c.jld2")
@load("data/deaths_data20210521.jld2")
@load("data/serological_data_with_20feb_to_31dec_age.jld2")
@load("data/serological_data_with_20feb_to_10Mar2021_age.jld2")
@load("data/N_kenya.jld2")

##Create a `CovAreaModel` object for Kenya which we will used with alternative models explaining the data
#`basic_nai_model` contains the PCR, serological and death data for Nairobi.

basic_nai_model = KenyaCoVSD.CoVAreaModel("Nairobi",θ -> 0.,θ -> 0.0;
                                            case_data = linelist_data_with_pos_neg,
                                            sero_data = serological_data,
                                            death_data = deaths_data,
                                            pop_data = N_kenya)

##Create a table of model score estimates

modelscores = DataFrame(model = String[],DIC = Union{Float64,String}[],ΔDIC = Union{Float64,String}[])

## Model 1: paper model
@load("modelfits/Nairobi_model.jld2")
nai_two_group = deepcopy(model)
DIC_two_group = KenyaCoVSD.modeldic(nai_two_group)
push!(modelscores,["Two-group SES",DIC_two_group,DIC_two_group-DIC_two_group])

## Model 2: One group SEIR with fitted R(t)
##Load code dependencies -- These define the transmission model for the one-group (daily varying R(t))
#They also define a log-likelihood function, priors and code for the EM algorithm using a combination
#of HMC inference (E-Step) and ADAM optimisation for fitting contact rate (ct) (M-step).
# R(t) = R₀ * ct (t)

include("one_group_model.jl");
include("one_group_inference.jl");
include("one_group_fitting_ct.jl");

##
# Set the ADAM optimisation parameters
baseline_adam_params = (β₁ = 0.9,
                        β₂ = 0.9,
                        η = 0.05,
                        ϵ = 1e-8,
                        total_num_steps = 500,
                        averaging_num_steps=100,
                        λ = 1e6)

                        ## Get first draw from MCMC with baseline ct estimate
## NB: The EM algorithm approach is slow, there are saved results accessed in the next code
inferparameters!(nai_one_group,2000,trans_one_group,0.05,D,cts₀)
# and get the initial E[log L]
final_LL = 0
for k = 1:size(nai_one_group.MCMC_results.chain,1)
        θ = NamedTuple{Tuple(keys(nai_one_group.MCMC_results.chain))}([nai_one_group.MCMC_results.chain[k,n,1] for n = 1:size(nai_one_group.MCMC_results.chain,2)])
        final_LL += ll_onegroup_newvariant_infboost_ct(θ,nai_one_group,0.0,cts₀)
end
final_LL = (final_LL/size(nai_one_group.MCMC_results.chain,1)) - ct_penalty(baseline_adam_params.λ,cts₀)

EM_steps = [(cts₀,final_LL)]

# Run alternate MCMC (for parameters not contact rate ct) and ADAM (for ct)
for k = 1:5
        log_ct_fit,mean_log_LL = ADAM_optim(log.(EM_steps[end][1][25:end]),
                                                random_grad_log_ct!,
                                                baseline_adam_params,
                                                nai_one_group,
                                                baseline_ct,
                                                trans_one_group,
                                                0.0)

        # Add each iteration to EM_steps
        push!(EM_steps,(vcat(baseline_ct,exp.(log_ct_fit)),mean_log_LL))
        inferparameters!(nai_one_group,500,trans_one_group,0.05,D,EM_steps[end][1])
end

#Uncomment to save the sequence of fitted contact
# @save("sensitivity_analysis/nai_EM_fit.jld2",nai_one_group,EM_steps)
## Plot the one-group model fit

@load("sensitivity_analysis/nairobi_fits/nai_EM_fit.jld2")
@load("forecasts/condensed_county_forecasts.jld2")
include("../analysis_scripts/fitting_methods.jl");
include("../analysis_scripts/plotting_methods.jl")

ct_fitted = EM_steps[end][1]
N = sum(N_kenya[:,"Nairobi"])
ct_plt = plot_ct(ct_fitted)

##Calculate predictions

fit_mats = gather_uncertainty_one_group(nai_one_group,ct_fitted)

#Join the % PCR pos and number PCR pos predictions based on whether each day has neg PCR test counts or not
nai_tests = vcat(nai_one_group.PCR_cases[1:(end-14),2],fill(-17,14+size(fit_mats.prop_PCR_pos_mat,1) - size(nai_one_group.PCR_cases,1)))
pred_num_PCR_pos_mat = fit_mats.prop_PCR_pos_mat.*nai_tests.*(nai_tests .>= 0) .+ fit_mats.no_neg_PCR_pos_mat.*(nai_tests .< 0)
#Get the posterior mean and credible intervals
pred_prop_sero_pos = get_credible_intervals(fit_mats.prop_sero_pos_mat)
pred_num_PCR_pos = get_credible_intervals(pred_num_PCR_pos_mat)
pred_incidence = get_credible_intervals(fit_mats.infection_mat)
#Compare to two-group fits
nai_fit = condensed_county_forecasts[[fit.name == "Nairobi" for fit in condensed_county_forecasts]][1]
sero_plt_compare = plot_pop_exposure(nai_fit,serological_data,serology_data,N_kenya);
PCR_plt_compare = plot_PCR(nai_fit,linelist_data_with_pos_neg,linelist_data);

plot!(PCR_plt_compare, 4:(length(pred_num_PCR_pos.pred)-3),
                        weekly_mv_av(pred_num_PCR_pos.pred),
                        title = "Nairobi PCR test, model comparison",
                        color = :green,lw = 3,
                        ribbon = (weekly_mv_av(pred_num_PCR_pos.lb),weekly_mv_av(pred_num_PCR_pos.ub)),
                        lab = "One group model: fit and forecast")

plot!(sero_plt_compare,pred_prop_sero_pos.pred./nai_one_group.N,
        title = "Nairobi population exposure, model comparison",
        lab = "One group model fit: seroposivity",lw = 3,color = :green,ls = :dash,
        ribbon = (pred_prop_sero_pos.lb./nai_one_group.N,pred_prop_sero_pos.ub./nai_one_group.N))

plot!(sero_plt_compare,cumsum(pred_incidence.pred,dims = 1)./nai_one_group.N,
        lab = "One group model fit: Overall population exposure",lw = 3,color = :red,ls = :dash,
        ribbon = (cumsum(pred_incidence.lb,dims = 1)./nai_one_group.N,cumsum(pred_incidence.ub,dims = 1)./nai_one_group.N )
        )

savefig(sero_plt_compare,"plots/sensitivity_analysis_plots/comparison_onegrp_vs_twogrp_pop_exposure.png")
savefig(PCR_plt_compare,"plots/sensitivity_analysis_plots/comparison_onegrp_vs_twogrp_PCR_pos.png")

## Fill table for model 2
@load("sensitivity_analysis/nairobi_fits/nai_EM_fit.jld2") #<---- Saved EM algorithm fits
DIC_one_group_model = KenyaCoVSD.modeldic(nai_one_group)

push!(modelscores,["One-group SES",DIC_one_group_model,DIC_one_group_model-DIC_two_group])

## Model 3: lower 25%, middle 50% and upper 25% wealth percentiles

#Include code for three group model, log-likelihood, prior
include("three_group_model.jl");

#Run MCMC for three group model
KenyaCoVSD.inferparameters!(nai_three_group,2000,trans_three_groups,0.05,D,q₀;serowaningrate = 0.0)
#uncomment to save
# @save("nai_three_group.jld2",nai_three_group)

##
@load("sensitivity_analysis/nairobi_fits/nai_three_group.jld2")

DIC_three_group_model = KenyaCoVSD.modeldic(nai_three_group)
push!(modelscores,["Three-group SES",DIC_three_group_model,DIC_three_group_model-DIC_two_group])

## Read-out the model-score table
display(modelscores)


text:
"""環境いろいろ"""

include("rmp_tree.jl")

#using .RMPTree: State

export ObsParam_point
export ObsParam_cylinder
export set_obs

using Random
using Parameters

include("utils.jl")

sd = Random.seed!(123)

### 目標位置 ###

"""点目標"""
@with_kw struct GoalParam_point{T}
    x::T
    y::T
    z::T
end

function _set_goal(p::GoalParam_point{T}) where T
    State(
        [p.x, p.y, p.z], zeros(T, 3)
    )
end

### 障害物 ###
@with_kw struct ObsParam_field{T, U}
    x::T
    y::T
    z::T
    lx::T
    ly::T
    alpha::T
    beta::T
    gamma::T
    n::U
end

"""
障害物（点）のパラメータ
"""
@with_kw struct ObsParam_point{T}
    x::T
    y::T
    z::T
end

"""障害物（点）を設置"""
function _set_obs(p::ObsParam_point{T}) where T
    obs = Vector{State{T}}(undef, 1)
    obs[1] = State(
        [p.x, p.y, p.z], zeros(T, 3)
    )
    return obs
end

"""
障害物（面）のパラメータ
"""
@with_kw struct ObsParam_plane{T, U}
    x::T
    y::T
    z::T
    lx::T
    ly::T
    alpha::T
    beta::T
    gamma::T
    n::U
end

"""
障害物（面）を設置
"""
function _set_obs(p::ObsParam_plane{T, U}) where {T, U}
    obs = Vector{State{T}}(undef, p.n)
    R = rotate(p.alpha, p.beta, p.gamma)
    t = [p.x, p.y, p.z]

    for i in 1:p.n
        X = [
            (rand(sd, T) - 0.5) * p.lx
            (rand(sd, T) - 0.5) * p.ly
            0.0
        ]

        obs[i] = State(
            R * X + t,
            zeros(T, 3)
        )
    end
    return obs
end

"""
障害物（球）のパラメータ

r : 半径
n : 障害物点の数
"""
@with_kw struct ObsParam_sphere{T, U}
    x::T
    y::T
    z::T
    r::T
    n::U
end

"""
障害物（球）を設置
"""
function _set_obs(p::ObsParam_sphere{T, U}) where {T, U}
    obs = Vector{State{T}}(undef, p.n)
    for i in 1:p.n
        theta = 2*rand(T)-1 |> acos
        phi = 2π * rand(T)
        X = [
            p.r * sin(theta) * cos(phi) + p.x
            p.r * sin(theta) * sin(phi) + p.y
            p.r * cos(theta) + p.z
        ]
        obs[i] = State(
            X, zeros(T, 3)
        )
    end
    return obs
end

"""
障害物（円筒）のパラメータ

r :
"""
@with_kw struct ObsParam_cylinder{T, U}
    x::T
    y::T
    z::T
    r::T
    L::T
    alpha::T
    beta::T
    gamma::T
    n::U
end

"""障害物（円筒）を設置"""
function _set_obs(p::ObsParam_cylinder{T}) where {T}
    obs = Vector{State{T}}(undef, p.n)
    R = rotate(p.alpha, p.beta, p.gamma)
    t = [p.x, p.y, p.z]
    for i in 1:p.n
        theta = 2π * rand(sd, T)
        X = [
            p.r * cos(theta)
            p.r * sin(theta)
            (p.L)*(rand(sd, T) - 1/2)
        ]
        obs[i] = State(
            R * X + t,
            zeros(T, 3)
        )
    end
    return obs
end


text:
using NPZ
using ArgParse

s = ArgParseSettings()
include("./args.jl")
s = args_data(s)
s = args_pq(s)

config = parse_args(s)

pq_codes = NPZ.npzread(config["pq_codes"]);
pq_codes .+= 1

pq_codebook = NPZ.npzread(config["pq_codebook"]);

M=config["M"]
Ks = config["Ks"]
D = config["D"]
Ds = Int(D/M)

if length(size(pq_codebook))==1
    pq_codebook_ = zeros(Float32, M, Ks, Ds)
    for i in 1:M
        cc = pq_codebook[(i-1)*Ks*Ds + 1 : i*Ks*Ds]
        cc =reshape(cc, Ds, Ks)
        cc = transpose(cc)
        pq_codebook_[i,:,:] = cc
    end
end
pq_codebook = pq_codebook_

queries = NPZ.npzread(config["queries"]);
tr100 = NPZ.npzread(config["tr100"])
tr100 .+= 1;

include("./neighbors.jl")
sp = SN_pq.SearchNeighbors_PQ(config["M"], config["Ks"], config["D"], pq_codebook, pq_codes, config["metric"]);

time = @elapsed neighbors_MA = SN_pq.get_neighbors(sp, queries, 512);
println(time);

include("./utils.jl")
recall_atN(neighbors_MA, tr100)

show(to)


text:
_transform(t, cols::NamedTuple) = transform(t, pairs(cols))
_transform(t, cols::Union{Columns, IndexedTables.AbstractIndexedTable}) = _transform(t, columns(cols))
_transform(t, args...) = transform(t, args...)

function transform_vec_helper(args...)
    d = gensym()
    func = Expr(:(->), d, Expr(:call, :(JuliaDBMeta._transform), d, with_helper(d, args[end])))
    Expr(:call, :(JuliaDBMeta._pipe), func, replace_keywords(args[1:end-1])...)
end

"""
`@transform_vec(d, x)`

Replace all symbols in expression `x` with the respective column in `d`: the result has to be
 a `NamedTuple` of vectors or a table and is horizontally merged with `d`. In this context,
`_` refers to the whole table `d`. To use actual symbols, escape them with `^`, as in `^(:a)`.
Use `{}` syntax for automatically named `NamedTuples`.
Use `cols(c)` to refer to column `c` where `c` is a variable that evaluates to a symbol. `c` must be available in
the scope where the macro is called.

## Examples

```jldoctest transform_vec
julia> t = table((a = [1,2,3], b = ["x","y","z"]));

julia> @transform_vec t {:a .+ 1}
Table with 3 rows, 3 columns:
a  b    a .+ 1
──────────────
1  "x"  2
2  "y"  3
3  "z"  4
```
"""
macro transform_vec(args...)
    esc(transform_vec_helper(args...))
end

function transform_helper(args...)
    d = gensym()
    func = Expr(:(->), d, Expr(:call, :(JuliaDBMeta._transform), d, map_helper(d, args[end])))
    Expr(:call, :(JuliaDBMeta._pipe_chunks), func, args[1:end-1]...)
end

"""
`@transform(d, x)`

Apply the expression `x` row by row in `d`: collect the result as a table
(elements returned by `x` must be `NamedTuples`) and merge it horizontally with `d`. In this context,
`_` refers to the whole row. To use actual symbols, escape them with `^`, as in `^(:a)`.

Use `{}` syntax for automatically named `NamedTuples`.
Use `cols(c)` to refer to field `c` where `c` is a variable that evaluates to a symbol. `c` must be available in
the scope where the macro is called.

## Examples

```jldoctest transform
julia> t = table((a = [1,2,3], b = ["x","y","z"]));

julia> @transform t {:a + 1}
Table with 3 rows, 3 columns:
a  b    a + 1
──────────────
1  "x"  2
2  "y"  3
3  "z"  4
```
"""
macro transform(d, x)
    esc(transform_helper(d, x))
end

macro transform(x)
    i = gensym()
    esc(Expr(:(->), i, transform_helper(i, x)))
end


text:
function __init__()
    global DEFAULT_RESOURCE = Ref{AbstractResource}(CPU1())

    # for testing asynchronous training of learing networks:
    global MACHINE_CHANNEL =
        RemoteChannel(() -> Channel(100), myid())

    TRAIT_FUNCTION_GIVEN_NAME[:measure]      = is_measure
    TRAIT_FUNCTION_GIVEN_NAME[:measure_type] = is_measure_type

    MLJModelInterface.set_interface_mode(MLJModelInterface.FullInterface())
end


text:
module TransferFunctionGPKernel

    export symGKGᵀ, testme

    using HCubature, Cuba, SimpleLinearPiecewiseApproximation, TransferFunctions, Printf
    using Distributed

    using SpecialFunctions, Interpolations

    include("integrals.jl")

    include("GKGT.jl")
    include("convkernel.jl")
    include("convkernel_cuba.jl")

    include("test_GKGT.jl")
    include("GKGT.jl")

end


text:
@doc doc"""
    Euclidean{T<:Tuple} <: Manifold

Euclidean vector space $\mathbb R^n$.

# Constructor

    Euclidean(n)

generates the $n$-dimensional vector space $\mathbb R^n$.

   Euclidean(m, n)

generates the $mn$-dimensional vector space $\mathbb R^{m \times n}$, whose
elements are interpreted as $m \times n$ matrices.
"""
struct Euclidean{T<:Tuple} <: Manifold where {T} end

Euclidean(n::Int) = Euclidean{Tuple{n}}()
Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}()

@generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...)

struct EuclideanMetric <: RiemannianMetric end

@traitimpl HasMetric{Euclidean,EuclideanMetric}

function local_metric(::MetricManifold{<:Manifold,EuclideanMetric}, x)
    return Diagonal(ones(SVector{size(x, 1),eltype(x)}))
end

function inverse_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x)
    return local_metric(M, x)
end

det_local_metric(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = one(eltype(x))

log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = zero(eltype(x))

inner(::Euclidean, x, v, w) = dot(v, w)
inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)

norm(::Euclidean, x, v) = norm(v)
norm(::MetricManifold{<:Manifold,EuclideanMetric}, x, v) = norm(v)

exp!(M::Euclidean, y, x, v) = (y .= x + v)

log!(M::Euclidean, v, x, y) = (v .= y - x)


text:
using DFTK
import DFTK: mpi_mean!
using Test
using Random
using MPI
include("testcases.jl")

@testset "Forces on silicon" begin
    function energy_forces(positions)
        model = model_DFT(silicon.lattice, silicon.atoms, positions, [:lda_x, :lda_c_pw])
        basis = PlaneWaveBasis(model; Ecut=7, kgrid=[2, 2, 2], kshift=[0, 0, 0],
                               symmetries_respect_rgrid=true,
                               fft_size=(18, 18, 18))  # FFT chosen to match QE
        is_converged = DFTK.ScfConvergenceDensity(1e-11)
        scfres = self_consistent_field(basis; is_converged)
        scfres.energies.total, compute_forces(scfres), compute_forces_cart(scfres)
    end

    # symmetrical positions, forces should be 0
    _, F0, _ = energy_forces([(ones(3)) / 8, -ones(3) / 8])
    @test norm(F0) < 1e-4

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8]  # displace a bit from equilibrium
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-5
    pos2 = [pos1[1] + ε * disp, pos1[2]]
    pos3 = [pos1[1] - ε * disp, pos1[2]]

    # second-order finite differences for accurate comparison; TODO switch the other tests to this too
    E1, F1, Fc1 = energy_forces(pos1)
    E2,  _,  _  = energy_forces(pos2)
    E3,  _,  _  = energy_forces(pos3)

    diff_findiff = -(E2 - E3) / (2ε)
    diff_forces = dot(F1[1], disp)
    @test abs(diff_findiff - diff_forces) < 1e-7

    # Rough test against QE reference (using PZ functional)
    reference = [[-5.81108123960e-3, -4.60222477677e-3, -3.37528153911e-3],
                 [ 5.81108123960e-3,  4.60222477677e-3,  3.37528153911e-3]]
    @test maximum(v -> maximum(abs, v), reference - Fc1) < 1e-5
end

@testset "Forces on silicon with spin and temperature" begin
    function silicon_energy_forces(positions; smearing=Smearing.FermiDirac())
        model = model_DFT(silicon.lattice, silicon.atoms, positions, :lda_xc_teter93;
                          temperature=0.03, smearing, spin_polarization=:collinear)
        basis = PlaneWaveBasis(model; Ecut=4, kgrid=[4, 1, 2], kshift=[1/2, 0, 0])

        n_bands = 10
        is_converged = DFTK.ScfConvergenceDensity(5e-10)
        scfres = self_consistent_field(basis; n_bands, is_converged)
        scfres.energies.total, compute_forces(scfres)
    end

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8]  # displace a bit from equilibrium
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-6
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    for (tol, smearing) in [(0.003, Smearing.FermiDirac()), (5e-5, Smearing.Gaussian())]
        E1, F1 = silicon_energy_forces(pos1; smearing)
        E2, _  = silicon_energy_forces(pos2; smearing)

        diff_findiff = -(E2 - E1) / ε
        diff_forces  = dot(F1[1], disp)
        @test abs(diff_findiff - diff_forces) < tol
    end
end

@testset "Forces on spin-polarised case" begin
    function oxygen_energy_forces(positions)
        O = ElementPsp(o2molecule.atnum, psp=load_psp("hgh/pbe/O-q6.hgh"))
        magnetic_moments = [1.0, 1.0]
        model = model_PBE(diagm([7.0, 7.0, 7.0]), [O, O], positions;
                          temperature=0.02, smearing=Smearing.Gaussian(), magnetic_moments)
        basis = PlaneWaveBasis(model; Ecut=4, kgrid=[1, 1, 1])

        scfres = self_consistent_field(basis;
                                       is_converged=DFTK.ScfConvergenceDensity(1e-7),
                                       ρ=guess_density(basis, magnetic_moments),
                                       damping=0.7, mixing=SimpleMixing())
        scfres.energies.total, compute_forces(scfres)
    end

    pos1 = [[0, 0, 0.1155], [0.01, -2e-3, -0.2]]
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-6
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    E1, F1 = oxygen_energy_forces(pos1)
    E2, _  = oxygen_energy_forces(pos2)

    diff_findiff = -(E2 - E1) / ε
    diff_forces  = dot(F1[1], disp)

    @test abs(diff_findiff - diff_forces) < 1e-4
end


text:
# ---
# title: 995. Minimum Number of K Consecutive Bit Flips
# id: problem995
# author: Indigo
# date: 2021-02-18
# difficulty: Hard
# categories: Greedy, Sliding Window
# link: <https://leetcode.com/problems/minimum-number-of-k-consecutive-bit-flips/description/>
# hidden: true
# ---
#
# In an array `A` containing only 0s and 1s, a _`K`-bit flip _consists of
# choosing a (contiguous) subarray of length `K` and simultaneously changing
# every 0 in the subarray to 1, and every 1 in the subarray to 0.
#
# Return the minimum number of `K`-bit flips required so that there is no 0 in
# the array.  If it is not possible, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [0,1,0], K = 1
#     Output: 2
#     Explanation: Flip A[0], then flip A[2].
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,1,0], K = 2
#     Output: -1
#     Explanation:  No matter how we flip subarrays of size 2, we can't make the array become [1,1,1].
#
#
# **Example 3:**
#
#
#
#     Input: A = [0,0,0,1,0,1,1,0], K = 3
#     Output: 3
#     Explanation:
#     Flip A[0],A[1],A[2]: A becomes [1,1,1,1,0,1,1,0]
#     Flip A[4],A[5],A[6]: A becomes [1,1,1,1,1,0,0,0]
#     Flip A[5],A[6],A[7]: A becomes [1,1,1,1,1,1,1,1]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 30000`
#   2. `1 <= K <= A.length`
#
#
## @lc code=start
using LeetCode, DataStructures

function min_k_bit_flips(A::Vector{Int}, K::Int)
    len, res = length(A), 0
    q = Queue{Int}()
    for i in 1:len
        (!isempty(q) && first(q) <= i - K) && dequeue!(q)

        if (length(q) & 1) == A[i]
            (i + K - 1 > len) && return -1
            enqueue!(q, i)
            res += 1
        end
    end
    res
end
## @lc code=end


text:
module _PermEn2D
export PermEn2D
    """
        Perm2D = PermEn2D(Mat)

    Returns the bidimensional permutation entropy estimate (`Perm2D`) estimated for
    the data matrix (`Mat`) using the default parameters: time delay = 1,
    logarithm = natural, template matrix size = [floor(H/10) floor(W/10)]
    (where H and W represent the height (rows) and width (columns) of the data matrix `Mat`) \n
    ** The minimum dimension size of Mat must be > 10.**

        Perm2D = PermEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
                                        tau::Int=1, Norm::Bool=true, Logx::Real=exp(1), Lock::Bool=true)

    Returns the bidimensional permutation entropy (`Perm2D`) estimates for the data
    matrix (`Mat`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Template submatrix dimensions, an integer scalar (i.e. the same
              height and width) or a two-element vector of integers
              [height, width] with a minimum value > 1.
              (default: [floor(H/10) floor(W/10)]) \n
    `tau`   - Time Delay, a positive integer       (default: 1) \n
    `Norm`  - Normalization of permutation entropy estimate, a boolean  (default: true) \n
    `Logx`  - Logarithm base, a positive scalar    (default: natural)  \n
    `Lock`  - By default, PermEn2D only permits matrices with a maximum
              size of 128 x 128 to prevent memory errors when storing data
              on RAM. e.g. For Mat = [200 x 200], m = 3, and tau = 1,
              SampEn2D creates a vector of 753049836 elements.
              To enable matrices greater than [128 x 128] elements,
              set `Lock` to false.  (default: true)  \n
              `WARNING: unlocking the permitted matrix size may cause your Julia
              IDE to crash.`\n\n

    **NOTE** - The original bidimensional permutation entropy algorithms
              [1][2] do not account for equal-valued elements of the embedding
              matrices.
              To overcome this, `PermEn2D` uses the lowest common rank for
              such instances. For example, given an embedding matrix A where,
              A = [3.4  5.5  7.3]
                  |2.1  6    9.9|
                  [7.3  1.1  2.1]
              would normally be mapped to an ordinal pattern like so,
              [3.4  5.5  7.3  2.1  6  9.9  7.3  1.1  2.1] =>
              [ 8    4    9    1   2   5    3    7    6 ]
              However, indices 4 & 9, and 3 & 7 have the same values, 2.1
              and 7.3 respectively. Instead, PermEn2D uses the ordinal pattern
              [ 8    4    4    1   2   5    3    3    6 ]
              where the lowest rank (4 & 3) are used instead (of 9 & 7).
              Therefore, the number of possible permutations is no longer
              (mx*my)!, but (mx*my)^(mx*my). Here, the PermEn2D value is
              normalized by the maximum Shannon entropy (Smax = log((mx*my)!)
              ``assuming that no equal values are found in the permutation
              motif matrices``, as presented in [1].

    #  See also `SampEn2D`, `FuzzEn2D`, `DispEn2D`, `DistEn2D`

    # References:
        [1] Haroldo Ribeiro et al.,
                "Complexity-Entropy Causality Plane as a Complexity Measure
                for Two-Dimensional Patterns"
                PLoS ONE (2012), 7(8):e40689,

        [2] Luciano Zunino and Haroldo Ribeiro,
                "Discriminating image textures with the multiscale
                two-dimensional complexity-entropy causality plane"
                Chaos, Solitons and Fractals,  91:679-688 (2016)

        [3] Matthew Flood and Bernd Grimm,
                "EntropyHub: An Open-source Toolkit for Entropic Time Series Analysis"
                PLoS ONE (2021) 16(11): e0259448.
    """
    function PermEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
        tau::Int=1, Norm::Bool=true, Logx::Real=exp(1), Lock::Bool=true)

    NL, NW = size(Mat)
    ((NL > 128 || NW > 128) && Lock) ?
        error("To prevent memory errors, matrix width & length must have <= 128 elements.
            To estimate SampEn2D  for the current matrix ($NL,$NW) change Lock to 'false'.
            Caution: unlocking the safe matrix size may cause the Julia IDE to crash.") :
            nothing

    length(m)==1 ?  (mL = m; mW = m) : (mL = m[1];  mW = m[2])
    (NL > 10 && NW > 10) ? nothing :
        error("Number of rows and columns in Mat must be > 10")
    (minimum(m)>1) ? nothing :
        error("m:     must be an integer > 1, or a 2 element tuple of integer values > 1")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (Logx>0) ? nothing :  error("Logx:  must be a positive number > 0")

    NL = NL - (mL-1)*tau
    NW = NW - (mW-1)*tau

    Temp = Mat[1:tau:(mL-1)*tau+1,1:tau:(mW-1)*tau+1]
    Sord = sort(Temp[:])
    Dix = sortperm(Temp[:])

    if any(diff(Sord).==0)
        for x in findall(diff(Sord).==0).+1
            Dix[x] = Dix[x-1]
        end
    end
    Counter = [0]

    for k in 1:NL
        for n in 1:NW
            Temp = Mat[k:tau:(mL-1)*tau+k,n:tau:(mW-1)*tau+n]
            Sord = sort(Temp[:])
            Dx = sortperm(Temp[:])

            if any(diff(Sord).==0)
                for x in findall(diff(Sord).==0).+1
                    Dx[x] = Dx[x-1]
                end
            end

            if any(all((Dix .- Dx).==0,dims=1))
                Counter .+= all((Dix .- Dx).==0,dims=1).*1
            else
                Dix = [Dix Dx]
                Counter = [Counter 1]
            end
        end
    end

    sum(Counter) != NL*NW ? @warn("Potential error with permutation comparisons.") : nothing

    Pi = Counter/sum(Counter)
    Perm2D = -sum(Pi.*log.(Logx, Pi))
    Norm ? Perm2D /= log(Logx, factorial(big(mL*mW))) : nothing

    return Perm2D
    end

end
"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
mutable struct Review <: GitHubType
    pr::Union{PullRequest, Nothing}
    id::Union{Int, Nothing}
    user::Union{Owner, Nothing}
    body::Union{String, Nothing}
    state::Union{String, Nothing}
end

function Review(pr::PullRequest, data::Dict)
    rev = json2github(Review, data)
    rev.pr = pr
    rev
end

namefield(rev::Review) = rev.id

@api_default function reviews(api::GitHubAPI, repo, pr::PullRequest; options...)
    path = "/repos/$(name(repo))/pulls/$(name(pr))/reviews"
    results, page_data = gh_get_paged_json(api, path; options...)
    return map(x->Review(pr, x), results), page_data
end

@api_default function comments(api::GitHubAPI, repo, rev::Review; options...)
    path = "/repos/$(name(repo))/pulls/$(name(rev.pr))/reviews/$(name(rev))/comments"
    results, page_data = gh_get_paged_json(api, path; options...)
    return map(Comment, results), page_data
end

@api_default function reply_to(api::GitHubAPI, repo, r::Review, c::Comment, body; options...)
    create_comment(api, repo, r.pr, :review; params = Dict(
        :body => body,
        :in_reply_to => c.id
    ), options...)
end

@api_default function dismiss_review(api::GitHubAPI, repo::Repo, r::Review; options...)
    gh_put(api, "/repos/$(name(repo))/pulls/$(name(rev.pr))/reviews/$(name(rev))/dismissals")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "tree_sitter_c"
version = v"0.16.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/tree-sitter/tree-sitter.git", "d8c3f472d23ad79f519651d5cf715b56467d35d0")
    GitSource("https://github.com/tree-sitter/tree-sitter-c.git", "6002fcd5e86bb1e8670157bb008b97dbaf656d95")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
echo 'cmake_minimum_required(VERSION 3.13)
project(treesitter)
set(CMAKE_C_STANDARD 99)

include_directories(tree-sitter/lib/include tree-sitter/lib/src)

add_library(treesitter SHARED tree-sitter/lib/src/lib.c)

include_directories(tree-sitter-c/src)
add_library(treesitter_c SHARED tree-sitter-c/src/parser.c)

install(TARGETS treesitter_c DESTINATION lib CONFIGURATIONS Release)
install(TARGETS treesitter DESTINATION lib CONFIGURATIONS Release)' > CMakeLists.txt

BUILD_FLAGS=(-DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
mkdir build && cd build
cmake .. "${BUILD_FLAGS[@]}"
make -j${nproc}
make install
mkdir ${prefix}/include
cp ../tree-sitter/lib/include/tree_sitter/api.h ${prefix}/include/
install_license ../tree-sitter/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libtreesitter_c", :libtreesitter_c),
    LibraryProduct("libtreesitter", :libtreesitter)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
@testset "pointer" begin

@testset "unsafe_load & unsafe_store!" begin

@eval struct LoadableStruct
    a::Int64
    b::UInt8
end
Base.one(::Type{LoadableStruct}) = LoadableStruct(1,1)
Base.zero(::Type{LoadableStruct}) = LoadableStruct(0,0)

@testset for T in (Int8, UInt16, Int32, UInt32, Int64, UInt64, Int128,
                   Float32, Float64,
                   LoadableStruct),
             cached in (false, true)
    d_a = CuArray(ones(T))
    d_b = CuArray(zeros(T))

    ptr_a = CUDAnative.DevicePtr{T,AS.Global}(pointer(d_a))
    ptr_b = CUDAnative.DevicePtr{T,AS.Global}(pointer(d_b))
    @test Array(d_a) != Array(d_b)

    let ptr_a=ptr_a, ptr_b=ptr_b #JuliaLang/julia#15276
        if cached && capability(dev) >= v"3.2"
            @on_device unsafe_store!(ptr_b, unsafe_cached_load(ptr_a))
        else
            @on_device unsafe_store!(ptr_b, unsafe_load(ptr_a))
        end
    end
    @test Array(d_a) == Array(d_b)
end

@testset "indexing" begin
    function kernel(src, dst)
        unsafe_store!(dst, CUDAnative.unsafe_cached_load(src, 4))
        return
    end

    T = Complex{Int8}

    src = CuArray([T(1) T(9); T(3) T(4)])
    dst = CuArray([0])

    @cuda kernel(
        CUDAnative.DevicePtr{T,AS.Global}(pointer(src)),
        CUDAnative.DevicePtr{T,AS.Global}(pointer(dst))
    )

    @test Array(src)[4] == Array(dst)[1]
end

end

end


text:
using ConvertBNG
using Test

const PREC_EN = 0.002       #  2 mm for easting, northing absolute accuracy
const PREC_LL = 0.0000001   #  8 dp for lat, lon absolute accuracy
const PREC_FL = 0.000000001 # 10 dp for floating point comparison

const AIRY_1830_SEMI_MAJOR = 6377563.396
const AIRY_1830_SEMI_MINOR = 6356256.909

approx(x, y, atol) = all(((xx, yy),) -> isapprox(xx, yy, atol=atol), zip(x, y))
approx_ll(x, y) = approx(x, y, PREC_LL)
approx_en(x, y) = approx(x, y, PREC_EN)
approx_fl(x, y) = approx(x, y, PREC_FL)

@testset "All tests" begin

    @testset "caister water tower" begin
        # From Transformations and OSGM15 User Guide - Annexe A
        # Caister Water Tower example
        # ETRS89 longitude, latitude: (1.716073973, 52.658007833)
        # ETRS89 eastings, northings: (651307.003, 313255.686)
        # Shifts eastings, northings: (102.801, -78.236)
        # OSGB36 eastings, northings: (651409.804, 313177.450)
        @test approx_en(convert_etrs89_to_osgb36(651307.003, 313255.686), (651409.804, 313177.450))
        @test approx_en(convert_bng(1.716073973, 52.658007833), (651409.804, 313177.450))
        @test approx_en(convert_etrs89(1.716073973, 52.658007833), (651307.003, 313255.686))
        @test approx_ll(convert_etrs89_to_ll(651307.003, 313255.686), (1.716073973, 52.658007833))
        @test approx_ll(convert_lonlat(651409.804, 313177.450), (1.716073973, 52.658007833))
        # Check OSTN15 shift retrival
        @test approx_fl(get_ostn_ref(651, 313), (102.787, -78.242, 44.236))
        @test approx_fl(ostn15_shifts(651307.003, 313255.686), (102.800790098, -78.235984076, 44.228405219))

        # From Transformations and OSGM15 User Guide - Annexe C
        # Caister Water Tower example
        # E,N to lon,lat using Airy 1830 ellipsoid
        @test approx_ll(convert_to_ll(651409.903, 313177.270, AIRY_1830_SEMI_MAJOR, AIRY_1830_SEMI_MINOR), (1.717921583, 52.65757031))
    end

    @testset "min/max bounds tests" begin
        # failed ostn15 array retrieval
        @test get_ostn_ref(999, 999)[1] === NaN
        @test get_ostn_ref(999, 999)[2] === NaN
        @test get_ostn_ref(999, 999)[3] === NaN

        # Test min/max easting, northing, longitude, latitude
        # bad max easting
        @test convert_lonlat(700000.1, 1240000.0)[1] === NaN
        @test convert_lonlat(700000.1, 1240000.0)[2] === NaN
        # bad max northing
        @test convert_lonlat(600000.0, 1250000.1)[1] === NaN
        @test convert_lonlat(600000.0, 1250000.1)[2] === NaN
        # bad lon
        @test convert_bng(181., 51.44533267)[1] === NaN
        @test convert_bng(181., 51.44533267)[2] === NaN
        # bad lat
        @test convert_bng(-0.32824866, -90.01)[1] === NaN
        @test convert_bng(-0.32824866, -90.01)[2] === NaN

        # Test extent bounds of check function
        # below min lon
        max_lon = 1.768960
        min_lon = -6.379880
        @test check(-6.379881, (min_lon, max_lon)) === NaN
        # below min lat
        max_lat = 55.811741
        min_lat = 49.871159
        @test check(49.871158, (min_lat, max_lat)) === NaN
        # above max lon
        max_lon = 1.768960
        min_lon = -6.379880
        @test check(1.768961, (min_lon, max_lon)) === NaN
        # above max lat
        max_lat = 55.811741
        min_lat = 49.871159
        @test check(55.811742, (min_lat, max_lat)) === NaN
    end

    @testset "OS latlon to OSGB36 testset" begin
        # TP01 to TP40
        @test approx_en(convert_bng(-6.299777520, 49.92226394),  (  91492.146,   11318.804))
        @test approx_en(convert_bng(-5.203046100, 49.96006138),  ( 170370.718,   11572.405))
        @test approx_en(convert_bng(-4.108645636, 50.43885826),  ( 250359.811,   62016.569))
        @test approx_en(convert_bng(-1.297822772, 50.57563665),  ( 449816.371,   75335.861))
        @test approx_en(convert_bng(-1.450514337, 50.93127938),  ( 438710.920,  114792.250))
        @test approx_en(convert_bng(-3.551283492, 51.40078220),  ( 292184.870,  168003.465))
        @test approx_en(convert_bng( 1.444547304, 51.37447026),  ( 639821.835,  169565.858))
        @test approx_en(convert_bng(-2.544076183, 51.42754743),  ( 362269.991,  169978.690))
        @test approx_en(convert_bng(-0.119925572, 51.48936565),  ( 530624.974,  178388.464))
        @test approx_en(convert_bng(-4.308524770, 51.85890896),  ( 241124.584,  220332.641))
        @test approx_en(convert_bng( 0.897243270, 51.89436637),  ( 599445.590,  225722.826))
        @test approx_en(convert_bng(-2.154586144, 52.25529382),  ( 389544.190,  261912.153))
        @test approx_en(convert_bng(-0.912489570, 52.25160951),  ( 474335.969,  262047.755))
        @test approx_en(convert_bng( 0.401535471, 52.75136687),  ( 562180.547,  319784.995))
        @test approx_en(convert_bng(-1.197476559, 52.96219109),  ( 454002.834,  340834.943))
        @test approx_en(convert_bng(-2.640493208, 53.34480280),  ( 357455.843,  383290.436))
        @test approx_en(convert_bng(-4.289180698, 53.41628516),  ( 247958.971,  393492.909))
        @test approx_en(convert_bng(-4.289177929, 53.41630925),  ( 247959.241,  393495.583))
        @test approx_en(convert_bng(-3.040454907, 53.77911026),  ( 331534.564,  431920.794))
        @test approx_en(convert_bng(-1.663791682, 53.80021520),  ( 422242.186,  433818.701))
        @test approx_en(convert_bng(-4.634521682, 54.08666318),  ( 227778.330,  468847.388))
        @test approx_en(convert_bng(-0.077731332, 54.11685144),  ( 525745.670,  470703.214))
        @test approx_en(convert_bng(-4.388491181, 54.32919541),  ( 244780.636,  495254.887))
        @test approx_en(convert_bng(-2.938277411, 54.89542340),  ( 339921.145,  556034.761))
        @test approx_en(convert_bng(-1.616576852, 54.97912274),  ( 424639.355,  565012.703))
        @test approx_en(convert_bng(-4.296490163, 55.85399953),  ( 256340.925,  664697.269))
        @test approx_en(convert_bng(-3.294792193, 55.92478266),  ( 319188.434,  670947.534))
        @test approx_en(convert_bng(-5.828366919, 57.00606696),  ( 167634.202,  797067.144))
        @test approx_en(convert_bng(-2.048560307, 57.13902519),  ( 397160.491,  805349.736))
        @test approx_en(convert_bng(-4.219263986, 57.48625001),  ( 267056.768,  846176.972))
        @test approx_en(convert_bng(-8.578544561, 57.81351838),  (   9587.909,  899448.996))
        @test approx_en(convert_bng(-7.592555606, 58.21262247),  (  71713.132,  938516.404))
        @test approx_en(convert_bng(-6.260914555, 58.51560361),  ( 151968.652,  966483.780))
        @test approx_en(convert_bng(-3.726310221, 58.58120461),  ( 299721.891,  967202.992))
        @test approx_en(convert_bng(-3.214540011, 59.03743871),  ( 330398.323, 1017347.016))
        @test approx_en(convert_bng(-4.417576746, 59.09335035),  ( 261596.778, 1025447.602))
        @test approx_en(convert_bng(-5.827993398, 59.09671617),  ( 180862.461, 1029604.114))
        @test approx_en(convert_bng(-1.625169661, 59.53470794),  ( 421300.525, 1072147.239))
        @test approx_en(convert_bng(-1.274869104, 59.85409914),  ( 440725.073, 1107878.448))
        @test approx_en(convert_bng(-2.073828228, 60.13308092),  ( 395999.668, 1138728.951))
    end

    @testset "OS OSGB36 to latlon testset" begin
        # TP01 to TP40
        @test approx_ll(convert_lonlat(  91492.146,   11318.804), (-6.299777520, 49.92226394))
        @test approx_ll(convert_lonlat( 170370.718,   11572.405), (-5.203046100, 49.96006138))
        @test approx_ll(convert_lonlat( 250359.811,   62016.569), (-4.108645636, 50.43885826))
        @test approx_ll(convert_lonlat( 449816.371,   75335.861), (-1.297822772, 50.57563665))
        @test approx_ll(convert_lonlat( 438710.920,  114792.250), (-1.450514337, 50.93127938))
        @test approx_ll(convert_lonlat( 292184.870,  168003.465), (-3.551283492, 51.40078220))
        @test approx_ll(convert_lonlat( 639821.835,  169565.858), ( 1.444547304, 51.37447026))
        @test approx_ll(convert_lonlat( 362269.991,  169978.690), (-2.544076183, 51.42754743))
        @test approx_ll(convert_lonlat( 530624.974,  178388.464), (-0.119925572, 51.48936565))
        @test approx_ll(convert_lonlat( 241124.584,  220332.641), (-4.308524770, 51.85890896))
        @test approx_ll(convert_lonlat( 599445.590,  225722.826), ( 0.897243270, 51.89436637))
        @test approx_ll(convert_lonlat( 389544.190,  261912.153), (-2.154586144, 52.25529382))
        @test approx_ll(convert_lonlat( 474335.969,  262047.755), (-0.912489570, 52.25160951))
        @test approx_ll(convert_lonlat( 562180.547,  319784.995), ( 0.401535471, 52.75136687))
        @test approx_ll(convert_lonlat( 454002.834,  340834.943), (-1.197476559, 52.96219109))
        @test approx_ll(convert_lonlat( 357455.843,  383290.436), (-2.640493208, 53.34480280))
        @test approx_ll(convert_lonlat( 247958.971,  393492.909), (-4.289180698, 53.41628516))
        @test approx_ll(convert_lonlat( 247959.241,  393495.583), (-4.289177929, 53.41630925))
        @test approx_ll(convert_lonlat( 331534.564,  431920.794), (-3.040454907, 53.77911026))
        @test approx_ll(convert_lonlat( 422242.186,  433818.701), (-1.663791682, 53.80021520))
        @test approx_ll(convert_lonlat( 227778.330,  468847.388), (-4.634521682, 54.08666318))
        @test approx_ll(convert_lonlat( 525745.670,  470703.214), (-0.077731332, 54.11685144))
        @test approx_ll(convert_lonlat( 244780.636,  495254.887), (-4.388491181, 54.32919541))
        @test approx_ll(convert_lonlat( 339921.145,  556034.761), (-2.938277411, 54.89542340))
        @test approx_ll(convert_lonlat( 424639.355,  565012.703), (-1.616576852, 54.97912274))
        @test approx_ll(convert_lonlat( 256340.925,  664697.269), (-4.296490163, 55.85399953))
        @test approx_ll(convert_lonlat( 319188.434,  670947.534), (-3.294792193, 55.92478266))
        @test approx_ll(convert_lonlat( 167634.202,  797067.144), (-5.828366919, 57.00606696))
        @test approx_ll(convert_lonlat( 397160.491,  805349.736), (-2.048560307, 57.13902519))
        @test approx_ll(convert_lonlat( 267056.768,  846176.972), (-4.219263986, 57.48625001))
        @test approx_ll(convert_lonlat(   9587.909,  899448.996), (-8.578544561, 57.81351838))
        @test approx_ll(convert_lonlat(  71713.132,  938516.404), (-7.592555606, 58.21262247))
        @test approx_ll(convert_lonlat( 151968.652,  966483.780), (-6.260914555, 58.51560361))
        @test approx_ll(convert_lonlat( 299721.891,  967202.992), (-3.726310221, 58.58120461))
        @test approx_ll(convert_lonlat( 330398.323, 1017347.016), (-3.214540011, 59.03743871))
        @test approx_ll(convert_lonlat( 261596.778, 1025447.602), (-4.417576746, 59.09335035))
        @test approx_ll(convert_lonlat( 180862.461, 1029604.114), (-5.827993398, 59.09671617))
        @test approx_ll(convert_lonlat( 421300.525, 1072147.239), (-1.625169661, 59.53470794))
        @test approx_ll(convert_lonlat( 440725.073, 1107878.448), (-1.274869104, 59.85409914))
        @test approx_ll(convert_lonlat( 395999.668, 1138728.951), (-2.073828228, 60.1330809281))
    end

    @testset "Array input-output" begin
        # Check results from single- and multiple-point methods are the same
        let n = 100, lons = 5.0.*(rand(n) .- 0.5), lats = 55 .+ 5.0.*rand(n)
            en = vcat([hcat(convert_bng(lon, lat)...) for (lon, lat) in zip(lons, lats)]...)
            # Broadcast identity comparison to check NaNs are the same
            @test all(convert_bng(hcat(lons, lats)) .=== en)
            ll = vcat([hcat(convert_lonlat(east, north)...) for (east, north) in eachrow(en)]...)
            @test all(convert_lonlat(en) .=== ll)
        end
    end

end


text:
abstract type          Supervised <: Model end
abstract type        Unsupervised <: Model end
abstract type           Annotator <: Model end

abstract type Probabilistic <: Supervised end
abstract type Deterministic <: Supervised end
abstract type      Interval <: Supervised end

abstract type JointProbabilistic <: Probabilistic end

abstract type Static                <: Unsupervised end

abstract type SupervisedAnnotator   <: Annotator end
abstract type UnsupervisedAnnotator <: Annotator end

abstract type UnsupervisedDetector <: UnsupervisedAnnotator end
abstract type SupervisedDetector   <: SupervisedAnnotator end

abstract type ProbabilisticSupervisedDetector   <: SupervisedDetector end
abstract type ProbabilisticUnsupervisedDetector <: UnsupervisedDetector end

abstract type DeterministicSupervisedDetector   <: SupervisedDetector end
abstract type DeterministicUnsupervisedDetector <: UnsupervisedDetector end

const ABSTRACT_MODEL_SUBTYPES =
    [:Supervised,
     :Unsupervised,
     :Probabilistic,
     :Deterministic,
     :Interval,
     :JointProbabilistic,
     :Static,
     :Annotator,
     :SupervisedAnnotator,
     :UnsupervisedAnnotator,
     :SupervisedDetector,
     :UnsupervisedDetector,
     :ProbabilisticSupervisedDetector,
     :ProbabilisticUnsupervisedDetector,
     :DeterministicSupervisedDetector,
     :DeterministicUnsupervisedDetector]

MLJBase.Unsupervised -> Model


text:
using ModelingToolkit
using Unitful
using LinearAlgebra
using ForwardDiff

# TODO: Check which of these are still in use and/or necessary.

# Avoid Plots compatibility issue with Unitful
Base.occursin(s::Any, q::Unitful.Quantity) = occursin(s, repr(q))
Base.pointer(q::Unitful.Quantity) = pointer(repr(q))

# Ensure that ForwardDiff values can get sent to C calls
Base.unsafe_convert(T::Type{<:Any}, x::ForwardDiff.Dual) = T(ForwardDiff.value.(x))

# Compute norms for arrays of symbolic variables (as needed by EphemerisNBP)
LinearAlgebra.norm(a::AbstractArray{<:ModelingToolkit.Num}) = sum(a .^ 2)^(1/2)

# --------------------------------#
# Common Subexpression Evaluation #
# --------------------------------#
using SymbolicUtils.Code: Let, Func, MakeArray, SetArray, (←), LiteralExpr, AtIndex
using Symbolics: Equation

# Code below copied from Shashi Gowda's work in https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/200
# Under the MIT License as per conditions below:

# Copyright (c) 2020: Shashi Gowda, Yingbo Ma, Mason Protter, Julia Computing.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

### Code-related utilities

using SymbolicUtils

using SymbolicUtils: Sym, Term, istree

using SymbolicUtils.Rewriters

using DataStructures

### Common subexprssion evaluation

newsym() = Sym{Number}(gensym("cse"))

function _cse(expr, dict=OrderedDict())
    r = @rule ~x::istree => haskey(dict, ~x) ? dict[~x] : dict[~x] = newsym()
    final = Postwalk(Chain([r]))(expr)
end

function cse(expr)
    !istree(expr) && return expr
    dict=OrderedDict()
    final = _cse(expr, dict)
    Let([var ← ex for (ex, var) in pairs(dict)], final)
end

function _cse(exprs::AbstractArray)
    dict = OrderedDict()
    final = map(ex->_cse(ex, dict), exprs)
    ([var ← ex for (ex, var) in pairs(dict)], final)
end

function cse(x::MakeArray)
    assigns, expr = _cse(x.elems)
    Let(assigns, MakeArray(expr, x.similarto, x.output_eltype))
end

function cse(x::SetArray)
    assigns, expr = _cse(x.elems)
    Let(assigns, SetArray(x.inbounds, x.arr, expr))
end
# --------------------------------#
# End of copyrighted code         #
# --------------------------------#

function _cse(eq::Num, dict)
    final = _cse(Symbolics.value(eq), dict)
    Num(final)
end

function _cse(eq::AtIndex, dict)
    final = _cse(eq.elem, dict)
    AtIndex(eq.i, final)
end

function _cse(eq::Equation, dict)
    lhs = _cse(eq.lhs, dict)
    rhs = _cse(eq.rhs, dict)
    Equation(lhs, rhs)
end

function cse(x::LiteralExpr)
    # XXX: Assumes that the LiteralExpr contains a LineNumber as the first argument
    LiteralExpr(quote
        $(cse.(x.ex.args[2:end])...)
    end)
end

function cse(x::Func)
    body = cse(x.body)
    return Func(x.args, x.kwargs, body)
end


text:
function shellsort!(a::Vector, isless=isless)
    let h = 1
        # calculate the first gap that is too large, i.e.
        # get h(k) such that h(k)=(3^k - 1)/2 > ⌈N/3⌉).
        while (h < length(a)/3) || (h-1 < length(a)/3 <= h)
            h = 3*h + 1
        end
        while h > 1
            # reduce gap
            h = (h-1) ÷ 3
            # iterate over every sublist
            for s = 1:h
                # do a modified insertion sort
                for i = h+s:h:length(a)
                    for j = i:-h:h+s
                        isless(a[j], a[j-h]) || break
                        swap!(a, j, j-h)
                    end
                end
            end
        end
    end
    a
end


text:
Base.flipsign(x::SSigned, y::SFloat16{Y}) where Y = flipsign(x, SInt16{Base.bitcast(Int16, Y::Float16)}())
Base.flipsign(x::SSigned, y::SFloat32{Y}) where Y = flipsign(x, SInt32{Base.bitcast(Int32, Y::Float32)}())
Base.flipsign(x::SSigned, y::SFloat64{Y}) where Y = flipsign(x, SInt{Base.bitcast(Int64, Y::Float64)}())
#Base.flipsign(x::SSigned, y::Real) = flipsign(x, -oftype(x, signbit(y)))

Base.copysign(x::SSigned, y::SSigned)  = flipsign(x, x ⊻ y)
Base.copysign(x::SSigned, y::SFloat16) = copysign(x, SInt16{Base.bitcast(Int16, Y::Float16)}())
Base.copysign(x::SSigned, y::SFloat32) = copysign(x, SInt32{Base.bitcast(Int32, Y::Float32)}())
Base.copysign(x::SSigned, y::SFloat64) = copysign(x, SInt{Base.bitcast(Int64, Y::Float64)}())
#Base.copysign(x::SSigned, y::SReal)    = copysign(x, -oftype(x, signbit(y)))

Base.:(==)(x::SSigned,   y::SUnsigned) = (x >= 0) & (unsigned(x) == y)
Base.:(==)(x::SUnsigned, y::SSigned  ) = (y >= 0) & (x == unsigned(y))
Base.:(<)( x::SSigned,   y::SUnsigned) = (x <  0) | (unsigned(x) <  y)
Base.:(<)( x::SUnsigned, y::SSigned  ) = (y >= 0) & (x <  unsigned(y))
Base.:(<=)(x::SSigned,   y::SUnsigned) = (x <  0) | (unsigned(x) <= y)
Base.:(<=)(x::SUnsigned, y::SSigned ) = (y >= 0) & (x <= unsigned(y))
# unsigned shift counts always shift in the same direction

Base.rem2pi(::SInt32{X}, r::RoundingMode) where X = SFloat64{rem2pi(X::Int32, r)}()
Base.rem2pi(::SInt64{X}, r::RoundingMode) where X = SFloat64{rem2pi(X::Int64, r)}()

Base.invmod(::SInteger, ::SIntegerZeroType) = throw(DomainError(m, "`m` must not be 0."))

for (ST,BT) in SI2BI
    @eval begin
        Base.:(>>)(x::$ST{X}, y::SUnsigned{Y}) where {X,Y} = $ST{>>(X::$BT, Y)}()
        Base.:(<<)(x::$ST{X}, y::SUnsigned{Y}) where {X,Y} = $ST{<<(X::$BT, Y)}()
        Base.:(>>>)(x::$ST{X}, y::SUnsigned{Y}) where {X,Y} = $ST{>>>(X::$BT, Y)}()
        Base.:(^)(::$ST{X}, y::SUnsigned{Y}) where {X,Y} = $ST{^(X::$BT, Y::$BT)}()
        @generated function Base.lcm(a::$ST{A}, b::$ST{B}) where {A,B}
            out = $ST{lcm(A::$BT,B::$BT)}()
            quote
                $out
            end
        end
        (/)(x::$ST, y::$ST) where {T<:Integer} = float(x) / float(y)
        function Base.invmod(n::$ST, m::$ST)
            g, x, y = gcdx(n, m)
            g != 1 && throw(DomainError((n, m), "Greatest common divisor is $g."))
            # Note that m might be negative here.
            # For unsigned T, x might be close to typemax; add m to force a wrap-around.
            r = mod(x + m, m)
            # The postcondition is: mod(r * n, m) == mod(T(1), m) && div(r, m) == 0
            r
        end
    end
end
Base.:(>>)(x::SInteger, y::SInt) =
    ifelse(0 <= y, x >> unsigned(y), x << unsigned(-y))
Base.:(<<)(x::SInteger, y::SInt) =
    ifelse(0 <= y, x << unsigned(y), x >> unsigned(-y))
Base.:(>>>)(x::SInteger, y::SInt) =
    ifelse(0 <= y, x >>> unsigned(y), x << unsigned(-y))

Base.:(<<)(x::SInteger, y::BaseInteger) = <<(values(x), y)

#
Base.:(>>)(x::SInteger, y::Int) = >>(values(x), y)
Base.:(>>)(x::BaseInteger, y::SInt) = >>(x, values(y))
Base.:(<<)(x::SInteger, y::Int) = <<(values(x), y)
Base.:(<<)(x::BaseInteger, y::SInt) = <<(x, values(y))
Base.:(>>>)(x::SInteger, y::Int) = >>(values(x), y)
Base.:(>>>)(x::BaseInteger, y::SInt) = >>(x, values(y))

#SIntegerNegOneType = Union{SInt8{-Int8(1)},SInt16{-Int16(1)},SInt32{-Int32(1)},SInt64{-1},SInt128{-Int128(1)}}

Base.div(x::SSigned, y::SUnsigned) = flipsign(signed(div(unsigned(abs(x)), y)), x)
Base.div(x::SUnsigned, y::SSigned) = unsigned(flipsign(signed(div(x, unsigned(abs(y)))), y))

Base.rem(x::SSigned, y::SUnsigned) = flipsign(signed(rem(unsigned(abs(x)), y)), x)
Base.rem(x::SUnsigned, y::SSigned) = rem(x, unsigned(abs(y)))

Base.fld(x::SSigned, y::SUnsigned) = div(x, y) - (signbit(x) & (rem(x, y) != 0))
Base.fld(x::SUnsigned, y::SSigned) = div(x, y) - (signbit(y) & (rem(x, y) != 0))

Base.bswap(x::SInt8) = x
Base.bswap(x::SUInt8) = x
# TODO: bswap other stuff (probably unnecessary, but why not)

@inline _gcdxrem(x::Tuple{A,B}) where {A<:SInteger,B<:SInteger} = (x[2], rem(x[1], x[2]))
@inline _gcdxsub(q::SInteger, x::Tuple{A,B}) where {A<:SInteger,B<:SInteger} = (x[2], x[1] - q*x[2])

@inline function _gcdxwhile(
    ab::Tuple{A,B},s::Tuple{S0,S1},t::Tuple{T0,T1}
   ) where {A<:SInteger,B<:SInteger,S0<:SInteger,S1<:SInteger,T0<:SInteger,T1<:SInteger}
    _gcdxwhile(_gcdxrem(ab),
               _gcdxsub(div(ab[1],ab[2]), s),
               _gcdxsub(div(ab[1],ab[2]), t))
end

@inline function _gcdxwhile(ab::Tuple{A,B},
                            s::Tuple{S0,S1},
                            t::Tuple{T0,T1}
   ) where {A<:SInteger,B<:SIntegerZeroType,
            S0<:SInteger,S1<:SInteger,
            T0<:SInteger,T1<:SInteger}
    if ab[1] < 0
        return (-ab[1], -s[1], -t[1])
    else
        return (ab[1], s[1], t[1])
    end
end
function Base.gcdx(a::SInteger, b::SInteger)
    _gcdxwhile((a, b), (oneunit(a), zero(b)), (zero(a), oneunit(b)))
end

# only alternating `A` and `B` so that they eventually get promoted
#Tuple{SInt64{_1} where _1,Union{SInt64{-1}, SInt64{1}},SInt64{0}}

"""
    gcd(::Union{Integer,SInteger}, Union{Integer,SInteger})

# Examples
```jldoctest
using StaticValues, Test

julia> gcd(SInt(6), SInt(9)) === SInt(3)
true

julia> gcd(SInt(6),SInt(-9)) === SInt(3)
true

julia> gcd(SInt(6), -9) === 3
true

julia> gcd(6,SInt(-9)) === 3
true

julia> @inferred(gcd(SInt(6), SInt(-9)))
3
```
"""
function Base.gcd(a::SInteger, b::SInteger)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u, v = _gcdwhile(unsigned(abs(a >> za)), unsigned(abs(b >> zb)))
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(eltype(a)) && Base.__throw_gcd_overflow(a, b)
    r % typeof(a)
end

@inline _gcdwhile(u::U, v::U) where {U<:SInteger} = u, v
function _gcdwhile(u::SInteger, v::SInteger)
    if u > v
        _gcdwhile(v, u)
    end
    _gcdwhile(u, (v-u) >> trailing_zeros(v-u))
end

Base.gcd(a::SIntegerZeroType, b::SInteger) = abs(b)
Base.gcd(a::SIntegerZeroType, b::SIntegerZeroType) = abs(b)
Base.gcd(a::SInteger,         b::SIntegerZeroType) = abs(a)

"""
# StaticValues Examples
```jldoctest
julia> powermod(SInt(2), SInt(6), SInt(5)) === SInt(4)
true
```
"""
#=
Base.powermod(x::SInteger, p::SInteger, m::SInteger) =
    oftype(m, _powermod(one(m), p, m, prevpow(oftype(p, 2), p),oftype(m, mod(x,m))))

function _powermod(r::SInteger, p::SInteger, m::SInteger, t::SInteger, b::SInteger)
    if p >= t
        _powermod(mod(widemul(r,b), m), p - t, m, t, b)
    else
        _powermod2(r, p, m, t >>> SOne, b)
    end
end

Base.powermod(x::SInteger, p::SIntegerZeroType, m::SInteger) = mod(one(m), m)
Base.powermod(x::SInteger, p::SIntegerNegOneType, m::SInteger) = powermod(invmod(x, m), -p, m)
Base.powermod(x::SInteger, p::SInteger, m::SIntegerOneType) = zero(m)
Base.powermod(x::SInteger, p::SInteger, m::SIntegerNegOneType) = zero(m)

function _powermod2(r::SInteger, p::SInteger, m::SInteger, t::SInteger, b::SInteger)
    if t <= 0
        return r
    else
        return _powermod(mod(widemul(r, r), m), p, m, t, b)
    end
end

# there's not static bigint, so can't widen past 128
function _powermod2(r::SInt128, p::SInteger, m::SInteger, t::SInteger, b::SInteger)
    if t <= 0
        return r
    else
        return _powermod(mod(r*r, m), p, m, t, b)
    end
end

function _powermod(r::SInt128, p::SInteger, m::SInteger, t::SInteger, b::SInteger)
    if p >= t
        _powermod(mod(r * b, m), p - t, m, t, b)
    else
        _powermod2(r, p, m, t >>> SOne, b)
    end
end
=#
#_prevpow2(x::Unsigned) = one(x) << unsigned((sizeof(x)<<3)-leading_zeros(x)-SOne)
#_prevpow2(x::Integer) = reinterpret(typeof(x), x < 0 ? -_prevpow2(unsigned(-x)) : _prevpow2(unsigned(x)))


text:
# using PyPlot

###############################################################################
# This code section is translated directly from Joe's matlab script. It uses
# the bent-cylinder equation from Stanton 1998 (II).
###############################################################################

freq = 10.0.^range(3., stop=6., length=1000);
# define contants
Theta = 0;
g = 1.0357; # rho2/rho1
rho1 = 1.025; # density of water
rho2 = g*rho1; # density of shrimp
h = 1.0279; # c2 / c1
c1 = 1500; # spped of sound in  seawater
c2 = h*c1; # speed of sound in shrimp

a = 1.9E-3; # in m
rho_c_div_L = 2.2;
L = 29E-3; #in m
Ldiva = 16;
a = L./Ldiva;

R_12 = ((rho2*c2)/(rho1*c1) - 1) ./ (rho2*c2/(rho2*c1) + 1);

k1 = 2*pi.* freq ./ c1;
k2 =  2*pi.* freq ./ c2;

T_12 = 2*(rho2*c2./(rho1*c1)) / (1 + (rho2*c2./(rho1*c1)));
T_21 = 2*(rho1*c1./(rho2*c2)) / (1 + (rho1*c1./(rho2*c2)));
rho_c = rho_c_div_L * L;

alpha_B = 0.8;
mu_term = -pi/2*k1*a ./ (k1*a .+ 0.4);

I_o = 1 .- T_12 * T_21 .* exp.(im*4*k2*a) .* exp.(im*mu_term);

f_shrimp = 1/2*sqrt(rho_c*a) * R_12 .* exp.(-im*2*k1*a) .* I_o .*
	exp.(-alpha_B*(2*Theta*rho_c / L)^2);

TS_s = 10*log10.(abs.(f_shrimp).^2);
TS_shrimp = TS_s;

# plot(freq,TS_s);
# axis([350E3 800E3 -100 -40]);

###############################################################################
# recreating above using DWBA
###############################################################################
n_segments = 20

phi = 1 / rho_c_div_L
curvature_angle = range(-phi / 2, stop=phi / 2, length=n_segments)
origin = (rho_c - a) * [sin(phi / 2), cos(phi / 2)]
x = (rho_c - a) .* sin.(curvature_angle) .- origin[1]
z = (rho_c - a) .* cos.(curvature_angle) .- origin[2]
# plot(x, z)

bent_cylinder = Scatterer(
	[x'; zeros(1, n_segments); z'],
	a * ones(n_segments),
	h * ones(n_segments),
	g * ones(n_segments),
)

dwba = freq_spectrum(bent_cylinder, freq[1], freq[end], c1, length(k1))
# plot(freq, TS_s)
# plot(dwba["freq"], dwba["TS"])

@test dwba["TS"][1:10] == [-146.4384727146246,-134.3981674296699,-127.35600801470818,
	-122.36054597020149,-118.48682943594616,-115.3228602668385,-112.64886618846771,
	-110.33366282447088,-108.29263365106235,-106.46800319058687]