text:
include("../aimajulia.jl");

using Test;
using Random;

using Main.aimajulia;

#The following probability tests are from the aima-python doctests

cpt = variable_node(aimajulia.burglary_network, "Alarm");
event = Dict([Pair("Burglary", true), Pair("Earthquake", true)]);

@test probability(cpt, true, event) == 0.95;

event["Burglary"] = false;

@test probability(cpt, false, event) == 0.71;

s = Dict([Pair("A", true),
        Pair("B", false),
        Pair("C", true),
        Pair("D", false)]);

@test consistent_with(s, Dict());

@test consistent_with(s, s);

@test !consistent_with(s, Dict(Pair("A", false)));

@test !consistent_with(s, Dict(Pair("D", true)));

p = ProbabilityDistribution(variable_name="Flip");
p["H"], p["T"] = 0.25, 0.75;

@test p["H"] == 0.25;

p = ProbabilityDistribution(variable_name="X", frequencies=Dict([Pair("lo", 125),
                                                                Pair("med", 375),
                                                                Pair("hi", 500)]));

@test ((p["lo"], p["med"], p["hi"]) == (0.125, 0.375, 0.5));

p = JointProbabilityDistribution(["X", "Y"]);
p[(1,1)] = 0.25;

@test p[(1,1)] == 0.25;

p[Dict([Pair("X", 0), Pair("Y", 1)])] = 0.5;

@test p[Dict([Pair("X", 0), Pair("Y", 1)])] == 0.5;

@test (event_values(Dict([Pair("A", 10), Pair("B", 9), Pair("C", 8)]), ["C", "A"]) == (8, 10));

@test (event_values((1, 2), ["C", "A"]) == (1, 2));

p = JointProbabilityDistribution(["X", "Y"]);
p[(0, 0)], p[(0, 1)], p[(1, 1)], p[(2, 1)] = 0.25, 0.5, 0.125, 0.125;

@test enumerate_joint(["Y"], Dict([Pair("X", 0)]), p) == 0.75;

@test enumerate_joint(["X"], Dict([Pair("Y", 2)]), p) == 0;

@test enumerate_joint(["X"], Dict([Pair("Y", 1)]), p) == 0.75;

@test show_approximation(enumerate_joint_ask("X", Dict([Pair("Y", 1)]), p)) == "0: 0.6667, 1: 0.1667, 2: 0.1667";

bn = BayesianNetworkNode("X", "Burglary", Dict([Pair(true, 0.2), Pair(false, 0.625)]));

@test probability(bn, false, Dict([Pair("Burglary", false), Pair("Earthquake", true)])) == 0.375;

@test probability(BayesianNetworkNode("W", "", 0.75), false, Dict([Pair("Random", true)])) == 0.25;

X = BayesianNetworkNode("X", "Burglary", Dict([Pair(true, 0.2), Pair(false, 0.625)]));

@test (sample(X, Dict([Pair("Burglary", false), Pair("Earthquake", true)])) in (true, false));

Z = BayesianNetworkNode("Z", "P Q", Dict([Pair((true, true), 0.2),
                                            Pair((true, false), 0.3),
                                            Pair((false, true), 0.5),
                                            Pair((false, false), 0.7)]));

@test (sample(Z, Dict([Pair("P", true), Pair("Q", false)])) in (true, false));

@test show_approximation(enumeration_ask("Burglary",
                                        Dict([Pair("JohnCalls", true), Pair("MaryCalls", true)]),
                                        aimajulia.burglary_network)) == "false: 0.7158, true: 0.2842";

@test show_approximation(elimination_ask("Burglary",
                                        Dict([Pair("JohnCalls", true), Pair("MaryCalls", true)]),
                                        aimajulia.burglary_network)) == "false: 0.7158, true: 0.2842";

# RandomDevice() does not allow seeding.

mt_rng = MersenneTwister(21);

p = rejection_sampling("Earthquake", Dict(), aimajulia.burglary_network, 1000, mt_rng);

@test ((p[true], p[false]) == (0.002, 0.998));

mt_rng = Random.seed!(mt_rng, 71);

p = likelihood_weighting("Earthquake", Dict(), aimajulia.burglary_network, 1000, mt_rng);

@test ((p[true], p[false]) == (0.0, 1.0));

mt_rng = Random.seed!(mt_rng, 1017);

@test (show_approximation(likelihood_weighting("Burglary",
                            Dict([Pair("JohnCalls", true), Pair("MaryCalls", true)]),
                            aimajulia.burglary_network,
                            10000,
                            mt_rng)) == "false: 0.718, true: 0.282");

umbrella_prior = [0.5, 0.5];
umbrella_transition = [[0.7, 0.3], [0.3, 0.7]];
umbrella_sensor =  [[0.9, 0.2], [0.1, 0.8]];
umbrella_hmm = HiddenMarkovModel(umbrella_transition, umbrella_sensor);
umbrella_evidence = [true, true, false, true, true];    # Umbrella observation sequence (Fig. 15.5b)

@test (repr(forward_backward(umbrella_hmm, umbrella_evidence, umbrella_prior)) ==
        "Array{Float64,1}[[0.646936, 0.353064], [0.867339, 0.132661], [0.820419, 0.179581], [0.307484, 0.692516], [0.820419, 0.179581], [0.867339, 0.132661]]");

umbrella_evidence = [true, false, true, false, true];

@test (repr(forward_backward(umbrella_hmm, umbrella_evidence, umbrella_prior)) ==
        "Array{Float64,1}[[0.587074, 0.412926], [0.717684, 0.282316], [0.2324, 0.7676], [0.607195, 0.392805], [0.2324, 0.7676], [0.717684, 0.282316]]");

umbrella_prior = [0.5, 0.5];
umbrella_transition = [[0.7, 0.3], [0.3, 0.7]];
umbrella_sensor =  [[0.9, 0.2], [0.1, 0.8]];
umbrella_hmm = HiddenMarkovModel(umbrella_transition, umbrella_sensor);
umbrella_evidence = [true, false, true, false, true];
e_t = false;
t = 4;
d = 2;

@test (repr(fixed_lag_smoothing(e_t, umbrella_hmm, d, umbrella_evidence; t=t)) == "[0.111111, 0.888889]");

d = 5;

@test (fixed_lag_smoothing(e_t, umbrella_hmm, d, umbrella_evidence; t=t) == nothing);

umbrella_evidence = [true, true, false, true, true];
e_t = true;
d = 1;

@test (repr(fixed_lag_smoothing(e_t, umbrella_hmm, d, umbrella_evidence; t=t)) == "[0.993865, 0.00613497]");

N = 10;
umbrella_evidence = true;
umbrella_transition = [[0.7, 0.3], [0.3, 0.7]];
umbrella_sensor = [[0.9, 0.2], [0.1, 0.8]];
umbrella_hmm = HiddenMarkovModel(umbrella_transition, umbrella_sensor);
s = particle_filtering(umbrella_evidence, N, umbrella_hmm);

@test length(s) == N;

@test all(state in ("A", "B") for state in s);

# Probability Distribution Example (p.493)
p = ProbabilityDistribution(variable_name="Weather");
p["sunny"] = 0.6;
p["rain"] = 0.1;
p["cloudy"] = 0.29;
p["snow"] = 0.01;

@test p["rain"] == 0.1;

# Joint Probability Distribution Example (Fig. 13.3)
p = JointProbabilityDistribution(["Toothache", "Cavity", "Catch"]);
p[(true, true, true)] = 0.108;
p[(true, true, false)] = 0.012;
p[(false, true, true)] = 0.072;
p[(false, true, false)] = 0.008;
p[(true, false, true)] = 0.016;
p[(true, false, false)] = 0.064;
p[(false, false, true)] = 0.144;
p[(false, false, false)] = 0.576;

@test p[(true, true, true)] == 0.108;

# P(Cavity | Toothache) example from page 500
probability_cavity = enumerate_joint_ask("Cavity", Dict([Pair("Toothache", true)]), p);

@test show_approximation(probability_cavity) == "false: 0.4, true: 0.6";

@test (0.6 - 0.001 < probability_cavity[true] < 0.6 + 0.001);

@test (0.4 - 0.001 < probability_cavity[false] < 0.4 + 0.001);

# Seed the RNG for Monte Carlo localization Base.Tests
mt_rng = MersenneTwister(sum(Vector{UInt8}("aima-julia")));

m = MonteCarloLocalizationMap([0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 1 0;
                                0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0;
                                1 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1 0;
                                0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0;
                                0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0;
                                0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0;
                                0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0;
                                0 0 1 1 1 1 1 0 0 0 1 1 1 0 1 1 0;
                                0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 0;
                                0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0;
                                0 0 1 1 1 1 1 0 0 0 1 1 1 0 0 1 0],
                                rng=mt_rng);

"""
    P_motion_sample(kinematic_state::Tuple, v::Tuple, w::Int64)

Return a sample from the possible kinematic states (using a single element
probability distribution).
"""
function P_motion_sample(kinematic_state::Tuple, v::Tuple, w::Int64)
    local position::Tuple = kinematic_state[1:2];
    local orientation::Int64 = kinematic_state[3];

    # Rotate the robot.
    orientation = (orientation + w) % 4;

    for i in 1:orientation
        v = (v[2], -v[1]);
    end

    position = (position[1] + v[1], position[2] + v[2]);

    return (position..., orientation);
end

"""
    P_sensor(x::Int64, y::Int64)

Return the conditional probability for the range sensor noise reading.
"""
function P_sensor(x::Int64, y::Int64)
    if (x == y)
        return 0.8;
    elseif (abs(x - y) <= 2)
        return 0.05;
    else
        return 0.0;
    end
end

a = Dict([Pair("v", (0, 0)), Pair("w", 0.0)]);
z = (2, 4, 1, 6);
S = monte_carlo_localization(a, z, 1000, P_motion_sample, P_sensor, m);
grid_1 = fill(0, 11, 17);

for (x, y, v) in S
    if ((0 < x <= 11) && (0 < y <= 17))
        grid_1[x, y] = grid_1[x, y] + 1;
    end
end

println("GRID 1:");
for x in 1:size(grid_1)[1]
    for y in 1:size(grid_1)[2]
        print(grid_1[x, y], " ");
    end
    println();
end
println();

a = Dict([Pair("v", (0, 1)), Pair("w", 0.0)]);
z = (2, 3, 5, 7);
S = monte_carlo_localization(a, z, 1000, P_motion_sample, P_sensor, m, S);
grid_2 = fill(0, 11, 17);

for (x, y, v) in S
    if ((0 < x <= 11) && (0 < y <= 17))
        grid_2[x, y] = grid_2[x, y] + 1;
    end
end

println("GRID 2:");
for x in 1:size(grid_2)[1]
    for y in 1:size(grid_2)[2]
        print(grid_2[x, y], " ");
    end
    println();
end
println();

@test (grid_2[7, 8] > 700);

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
function nominal_diameter(method,P,Hn,Vdotn)
        if method == "Bier" # default
            D = 0.176*(P/Hn)^0.466
        elseif method == "Sarkaria"
            D = 0.71*P^0.43/Hn^0.65
        elseif method == "Warnick"
            D = 0.72*P^0.43/Hn^0.63
        elseif method == "Moffat"
            D = 0.52*P^0.43/Hn^0.6
        # methods not depending on efficiency
        elseif method == "USBR"
            D = 1.517*Vdotn^0.5/Hn^0.25
        elseif method == "Fahlbusch"
            D = 1.12*Vdotn^0.45/Hn^0.12
        end
    return D*1.5
end
# for darcys friction factor Sami and Jain
function f_Darcy(NRe, pipe_roughness,D)
    f = 1/(1.8*log10(6.9/NRe+(pipe_roughness/D/3.7)^1.1))^2
    return f
end
# for Reynolds number
function renoldsnum(rho,v,D,mu)
    return rho*v*D/mu
end
#

ext: jl


text:
include(joinpath("..", "common", "utils.jl"))
using LinearAlgebra: Diagonal, diag

⊗ = kron

abstract type AbstractProblem end
struct Isothermal{T} <: AbstractProblem
    T₀::T
    R::T
    p₀::T
    g::T
    ρ₀::T
    H::T
    Cᵥ::T
    γ::T
    function Isothermal{T}(;
        T₀ = 300,
        R = 287,
        p₀ = 10^5,
        g = T(98 // 100),
    ) where {T}
        H = R * T₀ / g
        ρ₀ = p₀ / (R * T₀)
        Cᵥ = 5R / 2
        γ = R / Cᵥ + 1
        new{T}(T₀, R, p₀, g, ρ₀, H, Cᵥ, γ)
    end
end

p̄_(s::AbstractProblem, z) = s.p₀ * exp(-z / s.H)
ρ̄_(s::AbstractProblem, z) = s.ρ₀ * exp(-z / s.H)
Ē_(s::AbstractProblem, z) = ρ̄_(s, z) * (s.Cᵥ * s.T₀ + s.g * z)
δp_(s::AbstractProblem, z, δE, δρ) = (s.R / s.Cᵥ) * (δE - s.g * z * δρ)
c̄_(s::AbstractProblem, _) = sqrt(s.γ * s.p₀ / s.ρ₀)

A_(s, z) = [
    0 1 0
    -(s.R / s.Cᵥ)*s.g*z 0 (s.R/s.Cᵥ)
    0 (Ē_(s, z) + p̄_(s, z))/ρ̄_(s, z) 0
]

function element_operators(N, T = Float64)
    ξ, ω = lglpoints(T, N)
    D = spectralderivative(ξ)
    return (ξ = ξ, ω = ω, D = D)
end

function create_mesh(K, elem, z0, z1)
    ξ = elem.ξ

    # Number of pointss in element and polynomial order
    Nq = length(ξ)
    N = Nq - 1

    # Element size
    Δz = (z1 - z0) / K
    J = Δz / (ξ[Nq] - ξ[1])

    # shift ξ to go (0, 1)
    ξ01 = (ξ[1:(end - 1)] .- ξ[1]) / (ξ[Nq] - ξ[1])

    # cg to dg scatter matrix
    Q = scatter_matrix(N, K)

    # CG DOF locations
    zcg = [z0 .+ Δz * (ξ01 .+ (0:(K - 1))')[:]; z1]

    # DG DOF locations
    zdg = reshape(Q * zcg, Nq, K)

    W = Diagonal(elem.ω)
    I_KK = sparse(I, K, K)
    Wcg = Array(diag(Q' * (I_KK ⊗ (J * W)) * Q))

    return (zcg = zcg, zdg = zdg, scatter = Q, J = J, Δz = Δz, Wcg = Wcg)
end

"""
    element_tendency!(∂q, q, elem_operator, problem, z_elem)

Evaluate the element rhs `∂q` associated with `q` for the
`problem::AbstractProblem` using the DG `elem_operator`.
`q` and `∂q` are taken to be `NamedTuple`s with fields `(δρ, δw, δE)`
"""
function element_tendency!(∂q, q, O, J, s, z)
    ρ̄ = ρ̄_.(Ref(s), z)
    Ē = Ē_.(Ref(s), z)
    p̄ = p̄_.(Ref(s), z)

    δp = δp_.(Ref(s), z, q.δE, q.δρ)

    # ∫ ϕ ∂δρ = ∫ (∂_ξ ϕ) ρ̄ δw
    ∂q.δρ .= (O.D' * (O.ω .* ρ̄ .* q.δw)) ./ (J .* O.ω)

    #  ∫ ϕ ∂δw = -∫ (ϕ / ρ̄) (∂_ξ δp + g δρ)
    ∂q.δw .= -s.g * q.δρ ./ ρ̄ - (O.D * δp) ./ (J .* ρ̄)

    #  ∫ ϕ ∂δE = ∫ (∂_ξ ϕ) (Ē + p̄) δw
    ∂q.δE .= (O.D' * (O.ω .* (Ē + p̄) .* q.δw)) ./ (J .* O.ω)
end

function dg2cg_scatter(q, mesh, elem)
    Q = mesh.scatter
    ω = elem.ω
    for i in 1:length(q)
        q[i][:] .= Q * ((Q' * (mesh.J * ω .* q[i])[:]) ./ mesh.Wcg)
    end
end

function tendency!(∂q, q, t, problem, mesh, elem, bc, forcing = nothing)
    # Evaluate the volume terms
    element_tendency!(∂q, q, elem, mesh.J, problem, mesh.zdg)

    # Add in MMS formcing
    if !isnothing(forcing)
        ∂q.δρ .+= forcing.ρ.(mesh.zdg, t)
        ∂q.δw .+= forcing.w.(mesh.zdg, t)
        ∂q.δE .+= forcing.E.(mesh.zdg, t)
    end

    # Use upwind boundary treatment
    if !isnothing(bc[1])
        # Get the boundary velocity and reference density
        δw0 = bc[1].δw(t)
        ρ̄0 = ρ̄_(problem, mesh.zdg[1])

        # Set the plus and minus states
        q⁻ = [q.δρ[1], ρ̄0 * q.δw[1], q.δE[1]]
        q⁺ = [q⁻[1], -q⁻[2] + 2ρ̄0 * δw0, q⁻[3]]

        # Evaluate the upwind flux
        f0 = bc[1].A⁻ * q⁻ + bc[1].A⁺ * q⁺

        # Evaluate the minus-side flux for w (strong derivative)
        flw = bc[1].A[2, :]' * q⁻

        # Apply lift
        ∂q.δρ[1] += f0[1] / (mesh.J * elem.ω[1])
        ∂q.δw[1] += (f0[2] - flw) / (mesh.J * elem.ω[1] * ρ̄0)
        ∂q.δE[1] += f0[3] / (mesh.J * elem.ω[1])
    end

    if !isnothing(bc[2])
        # Get the boundary velocity and reference density
        δw1 = bc[2].δw(t)
        ρ̄1 = ρ̄_(problem, mesh.zdg[end])

        # Set the plus and minus states
        q⁻ = [q.δρ[end], ρ̄1 * q.δw[end], q.δE[end]]
        q⁺ = [q⁻[1], -q⁻[2] + 2ρ̄1 * δw1, q⁻[3]]

        # Evaluate the upwind flux
        f1 = bc[2].A⁻ * q⁻ + bc[2].A⁺ * q⁺

        # Evaluate the minus-side flux for w (strong derivative)
        flw = bc[2].A[2, :]' * q⁻

        # Apply lift
        ∂q.δρ[end] += f1[1] / (mesh.J * elem.ω[end])
        ∂q.δw[end] += (f1[2] - flw) / (mesh.J * elem.ω[end] * ρ̄1)
        ∂q.δE[end] += f1[3] / (mesh.J * elem.ω[end])
    end

    # DG -> CG projection and then scatter for storage
    dg2cg_scatter(∂q, mesh, elem)
end


text:
## generic value

# TODO: this is a _very_ ugly wrapper, but hard to improve since we can't deduce the type
#       of a GenericValue, and need to pass concrete LLVM type objects to the API

export GenericValue, dispose,
       intwidth

@checked struct GenericValue
    ref::API.LLVMGenericValueRef
end
reftype(::Type{GenericValue}) = API.LLVMGenericValueRef

GenericValue(typ::IntegerType, N::Signed) =
    GenericValue(
        API.LLVMCreateGenericValueOfInt(ref(typ),
                                        reinterpret(Culonglong, convert(Int64, N)), True))

GenericValue(typ::IntegerType, N::Unsigned) =
    GenericValue(
        API.LLVMCreateGenericValueOfInt(ref(typ),
                                        reinterpret(Culonglong, convert(UInt64, N)), False))

intwidth(val::GenericValue) = API.LLVMGenericValueIntWidth(ref(val))

Base.convert(::Type{T}, val::GenericValue) where {T<:Signed} =
    convert(T, reinterpret(Clonglong, API.LLVMGenericValueToInt(ref(val), True)))

Base.convert(::Type{T}, val::GenericValue) where {T<:Unsigned} =
    convert(T, API.LLVMGenericValueToInt(ref(val), False))

GenericValue(typ::FloatingPointType, N::AbstractFloat) =
    GenericValue(API.LLVMCreateGenericValueOfFloat(ref(typ), convert(Cdouble, N)))

# NOTE: this ugly three-arg convert is needed to match the C API,
#       which uses the type to call the correct C++ function.

Base.convert(::Type{T}, val::GenericValue, typ::LLVMType) where {T<:AbstractFloat} =
    convert(T, API.LLVMGenericValueToFloat(ref(typ), ref(val)))

GenericValue(ptr::Ptr) =
    GenericValue(API.LLVMCreateGenericValueOfPointer(convert(Ptr{Cvoid}, ptr)))

Base.convert(::Type{Ptr{T}}, val::GenericValue) where {T} =
    convert(Ptr{T}, API.LLVMGenericValueToPointer(ref(val)))

dispose(val::GenericValue) = API.LLVMDisposeGenericValue(ref(val))

## execution engine

export ExecutionEngine, Interpreter, JIT,
       run

@checked struct ExecutionEngine
    ref::API.LLVMExecutionEngineRef
end
reftype(::Type{ExecutionEngine}) = API.LLVMExecutionEngineRef

# NOTE: these takes ownership of the module
function ExecutionEngine(mod::Module)
    out_ref = Ref{API.LLVMExecutionEngineRef}()
    out_error = Ref{Cstring}()
    status = convert(Core.Bool, API.LLVMCreateExecutionEngineForModule(out_ref, ref(mod),
                                                                 out_error))

    if status
        error = unsafe_string(out_error[])
        API.LLVMDisposeMessage(out_error[])
        throw(LLVMException(error))
    end

    return ExecutionEngine(out_ref[])
end
function Interpreter(mod::Module)
    API.LLVMLinkInInterpreter()

    out_ref = Ref{API.LLVMExecutionEngineRef}()
    out_error = Ref{Cstring}()
    status = convert(Core.Bool, API.LLVMCreateInterpreterForModule(out_ref, ref(mod),
                                                             out_error))

    if status
        error = unsafe_string(out_error[])
        API.LLVMDisposeMessage(out_error[])
        throw(LLVMException(error))
    end

    return ExecutionEngine(out_ref[])
end
function JIT(mod::Module, optlevel::API.LLVMCodeGenOptLevel=API.LLVMCodeGenLevelDefault)
    API.LLVMLinkInMCJIT()

    out_ref = Ref{API.LLVMExecutionEngineRef}()
    out_error = Ref{Cstring}()
    status = convert(Core.Bool, API.LLVMCreateJITCompilerForModule(out_ref, ref(mod),
                                                                   optlevel, out_error))

    if status
        error = unsafe_string(out_error[])
        API.LLVMDisposeMessage(out_error[])
        throw(LLVMException(error))
    end

    return ExecutionEngine(out_ref[])
end

dispose(engine::ExecutionEngine) = API.LLVMDisposeExecutionEngine(ref(engine))

for x in [:ExecutionEngine, :Interpreter, :JIT]
    @eval function $x(f::Core.Function, args...)
        engine = $x(args...)
        try
            f(engine)
        finally
            dispose(engine)
        end
    end
end

Base.push!(engine::ExecutionEngine, mod::Module) = API.LLVMAddModule(engine.ref, mod.ref)

# up to user to free the deleted module
function Base.delete!(engine::ExecutionEngine, mod::Module)
    out_ref = Ref{API.LLVMModuleRef}()
    API.LLVMRemoveModule(engine.ref, mod.ref, out_ref, Ref{Cstring}()) # out string is not used
    return Module(out_ref[])
end

Base.run(engine::ExecutionEngine, f::Function, args::Vector{GenericValue}=GenericValue[]) =
    GenericValue(API.LLVMRunFunction(ref(engine), ref(f),
                                     length(args), ref.(args)))

# ExectutionEngine function lookup

export functions

struct ExecutionEngineFunctionSet
    engine::ExecutionEngine
end

functions(engine::ExecutionEngine) = ExecutionEngineFunctionSet(engine)

Base.IteratorSize(::ExecutionEngineFunctionSet) = Base.SizeUnknown()
Base.iterate(::ExecutionEngineFunctionSet) =
    error("Iteration of functions in the execution engine is not supported")

function Base.get(functionset::ExecutionEngineFunctionSet, name::String, default)
    out_ref = Ref{API.LLVMValueRef}()
    API.LLVMFindFunction(functionset.engine.ref, name, out_ref)
    status = convert(Core.Bool, API.LLVMFindFunction(functionset.engine.ref, name, out_ref))
    return status == 0 ? Function(out_ref[]) : default
end

function Base.haskey(functionset::ExecutionEngineFunctionSet, name::String)
    f = get(functionset, name, nothing)
    return f != nothing
end

function Base.getindex(functionset::ExecutionEngineFunctionSet, name::String)
    f = get(functionset, name, nothing)
    return f == nothing ? throw(KeyError(name)) : f
end


text:
module SortBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools

const SUITE = BenchmarkGroup()
const LIST_SIZE = 50000
const LISTS = (
    ("ascending", collect(1:LIST_SIZE)),
    ("descending", collect(LIST_SIZE:-1:1)),
    ("ones", ones(LIST_SIZE)),
    ("random", samerand(LIST_SIZE))
)

#####################################
# QuickSort/MergeSort/InsertionSort #
#####################################

for (group, Alg) in (("quicksort", QuickSort), ("mergesort", MergeSort), ("insertionsort", InsertionSort))
    g = addgroup!(SUITE, group)
    for (kind, list) in LISTS
        ix = collect(1:length(list))
        g["sort forwards", kind] = @benchmarkable sort($list; alg = $Alg)
        g["sort reverse", kind] = @benchmarkable sort($list; alg = $Alg, rev = true)
        g["sortperm forwards", kind] = @benchmarkable sortperm($list; alg = $Alg)
        g["sortperm reverse", kind] = @benchmarkable sortperm($list; alg = $Alg, rev = true)
        g["sort! forwards", kind] = @benchmarkable sort!(x; alg = $Alg) setup=(x = copy($list))
        g["sort! reverse", kind] = @benchmarkable sort!(x; alg = $Alg, rev = true) setup=(x = copy($list))
        g["sortperm! forwards", kind] = @benchmarkable sortperm!(x, $list; alg = $Alg) setup=(x = copy($ix))
        g["sortperm! reverse", kind] = @benchmarkable sortperm!(x, $list; alg = $Alg, rev = true) setup=(x = copy($ix))
    end
    for b in values(g)
        b.params.time_tolerance = 0.30
    end
end

############
# issorted #
############

g = addgroup!(SUITE, "issorted")

for (kind, list) in LISTS
    g["forwards", kind] = @benchmarkable issorted($list)
    g["reverse", kind] = @benchmarkable issorted($list; rev = true)
end

for b in values(g)
    b.params.time_tolerance = 0.30
end

end # module


text:
module SearchLightSQLite

import SQLite, DataFrames, Logging
import SearchLight
import DBInterface

#
# Setup
#

const COLUMN_NAME_FIELD_NAME = :name

function SearchLight.column_field_name()
  COLUMN_NAME_FIELD_NAME
end

const DatabaseHandle = SQLite.DB
const ResultHandle   = Union{Vector{Any}, DataFrames.DataFrame, Vector{Tuple}, Vector{Tuple{Int64}}}

const TYPE_MAPPINGS = Dict{Symbol,Symbol}( # Julia => SQLite
  :char       => :TEXT,
  :string     => :TEXT,
  :text       => :TEXT,
  :integer    => :INTEGER,
  :int        => :INTEGER,
  :float      => :FLOAT,
  :decimal    => :DECIMAL,
  :datetime   => :DATETIME,
  :timestamp  => :INTEGER,
  :time       => :TIME,
  :date       => :DATE,
  :binary     => :BLOB,
  :boolean    => :BOOLEAN,
  :bool       => :BOOLEAN
)

const SELECT_LAST_ID_QUERY_START = "; SELECT CASE WHEN last_insert_rowid() = 0 THEN"
const SELECT_LAST_ID_QUERY_END = "ELSE last_insert_rowid() END AS $(SearchLight.LAST_INSERT_ID_LABEL)"

const CONNECTIONS = DatabaseHandle[]

#
# Connection
#

"""
    connect(conn_data::Dict)::DatabaseHandle
    function connect()::DatabaseHandle

Connects to the database defined in conn_data["database"] or conn_data["host"] and returns a handle.
If no conn_data is provided, a temporary, in-memory database will be used.
"""
function SearchLight.connect(conn_data::Dict = SearchLight.config.db_config_settings) :: DatabaseHandle
    dbname =  if haskey(conn_data, "host") && conn_data["host"] !== nothing
                conn_data["host"]
              elseif haskey(conn_data, "database") && conn_data["database"] !== nothing
                conn_data["database"]
              end

  isempty(dirname(dbname)) || mkpath(dirname(dbname))

  push!(CONNECTIONS, SQLite.DB(dbname))[end]
end

"""
    disconnect(conn::DatabaseHandle)::Nothing

Disconnects from database.
"""
function SearchLight.disconnect(conn::DatabaseHandle = SearchLight.connection()) :: Nothing
  conn = nothing
end

function SearchLight.connection()
  isempty(CONNECTIONS) && throw(SearchLight.Exceptions.NotConnectedException())

  CONNECTIONS[end]
end

#
# Data sanitization
#

"""
    escape_column_name(c::String, conn::DatabaseHandle)::String

Escapes the column name using native features provided by the database backend.

# Examples
```julia
julia>
```
"""
function SearchLight.escape_column_name(c::String, conn::DatabaseHandle = SearchLight.connection()) :: String
  join([SQLite.esc_id(c) for c in split(c, '.')], '.')
end

function SearchLight.escape_value(v::T, _ = nothing)::T where {T} # TODO: deprecate in next major version
  isa(v, Number) ? v : "'$(replace(string(v), "'"=>"''"))'"
end

#
# Query execution
#

"""
    query(sql::String, suppress_output::Bool, conn::DatabaseHandle)::DataFrames.DataFrame

Executes the `sql` query against the database backend and returns a DataFrame result.

# Examples:
```julia
julia> query(SearchLight.to_fetch_sql(Article, SQLQuery(limit = 5)), false, Database.connection)

2017-01-16T21:36:21.566 - info: SQL QUERY: SELECT \"articles\".\"id\" AS \"articles_id\", \"articles\".\"title\" AS \"articles_title\", \"articles\".\"summary\" AS \"articles_summary\", \"articles\".\"content\" AS \"articles_content\", \"articles\".\"updated_at\" AS \"articles_updated_at\", \"articles\".\"published_at\" AS \"articles_published_at\", \"articles\".\"slug\" AS \"articles_slug\" FROM \"articles\" LIMIT 5

  0.000985 seconds (16 allocations: 576 bytes)

5×7 DataFrames.DataFrame
...
```
"""
function SearchLight.query(sql::String, conn::DatabaseHandle = SearchLight.connection()) :: DataFrames.DataFrame
  parts::Vector{String} = if occursin(SELECT_LAST_ID_QUERY_START, sql)
                            split(sql, SELECT_LAST_ID_QUERY_START)
                          else
                            String[sql]
                          end

  if length(parts) == 2
    @info parts[1]
    DBInterface.execute(conn, parts[1])

    parts[2] = string(SELECT_LAST_ID_QUERY_START, parts[2])
    @info parts[2]

    DBInterface.execute(conn, parts[2]) |> DataFrames.DataFrame
  else
    @info parts[1]
    DBInterface.execute(conn, parts[1]) |> DataFrames.DataFrame
  end
end

function SearchLight.to_find_sql(m::Type{T}, q::SearchLight.SQLQuery, joins::Union{Nothing,Vector{SearchLight.SQLJoin}} = nothing)::String where {T<:SearchLight.AbstractModel}
  sql::String = ( string("$(SearchLight.to_select_part(m, q.columns, joins)) $(SearchLight.to_from_part(m)) $(SearchLight.to_join_part(m, joins)) $(SearchLight.to_where_part(q.where)) ",
                      "$(SearchLight.to_group_part(q.group)) $(SearchLight.to_having_part(q.having)) $(SearchLight.to_order_part(m, q.order)) ",
                      "$(SearchLight.to_limit_part(q.limit)) $(SearchLight.to_offset_part(q.offset))")) |> strip
  replace(sql, r"\s+"=>" ")
end

function SearchLight.to_store_sql(m::T; conflict_strategy = :error)::String where {T<:SearchLight.AbstractModel}
  uf = SearchLight.persistable_fields(typeof(m))
  insert_condition = (! SearchLight.ispersisted(m) || (SearchLight.ispersisted(m) && conflict_strategy == :update))

  sql = if insert_condition
    pos = findfirst(x -> x == SearchLight.primary_key_name(typeof(m)), uf)
    pos > 0 && splice!(uf, pos)

    fields = SearchLight.SQLColumn(uf)
    vals = join( map(x -> string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))), uf), ", ")

    "INSERT $(conflict_strategy == :ignore ? " OR IGNORE" : "") INTO $(SearchLight.table(typeof(m))) ($fields) VALUES ($vals)"
  else
    "UPDATE $(SearchLight.table(typeof(m))) SET $(SearchLight.update_query_part(m))"
  end

  id_sql = if insert_condition
    string(
      SELECT_LAST_ID_QUERY_START, raw" ",
      getfield(m, Symbol(SearchLight.primary_key_name(m))).value === nothing ? -1 : getfield(m, Symbol(SearchLight.primary_key_name(m))).value, raw" ",
      SELECT_LAST_ID_QUERY_END
    )
  else
    "; SELECT $( getfield(m, Symbol(SearchLight.primary_key_name(m))).value ) AS LAST_INSERT_ID"
  end

  string( sql, raw" ", id_sql)
end

function SearchLight.delete_all(m::Type{T}; truncate::Bool = true, reset_sequence::Bool = true, cascade::Bool = false)::Nothing where {T<:SearchLight.AbstractModel}
  "DELETE FROM $(SearchLight.table(m))" |> SearchLight.query

  nothing
end

function SearchLight.delete(m::T)::T where {T<:SearchLight.AbstractModel}
  SearchLight.ispersisted(m) || throw(SearchLight.Exceptions.NotPersistedException(m))

  "DELETE FROM $(SearchLight.table(typeof(m))) WHERE $(SearchLight.primary_key_name(typeof(m))) = '$(m.id.value)'" |> SearchLight.query

  m.id = SearchLight.DbId()

  m
end

function Base.count(m::Type{T}, q::SearchLight.SQLQuery = SearchLight.SQLQuery())::Int where {T<:SearchLight.AbstractModel}
  count_column = SearchLight.SQLColumn("COUNT(*) AS __cid", raw = true)
  q = SearchLight.clone(q, :columns, push!(q.columns, count_column))

  SearchLight.DataFrame(m, q)[1, Symbol("__cid")]
end

function SearchLight.update_query_part(m::T)::String where {T<:SearchLight.AbstractModel}
  update_values = join(map(x -> "$(string(SearchLight.SQLColumn(x))) = $(string(SearchLight.to_sqlinput(m, Symbol(x), getfield(m, Symbol(x)))) )", SearchLight.persistable_fields(typeof(m))), ", ")

  " $update_values WHERE $(SearchLight.table(typeof(m))).$(SearchLight.primary_key_name(typeof(m))) = '$(m.id.value)'"
end

function SearchLight.column_data_to_column_name(column::SearchLight.SQLColumn, column_data::Dict{Symbol,Any}) :: String
  "$(SearchLight.to_fully_qualified(column_data[:column_name], column_data[:table_name])) AS $(isempty(column_data[:alias]) ? SearchLight.to_sql_column_name(column_data[:column_name], column_data[:table_name]) : column_data[:alias] )"
end

function SearchLight.to_from_part(m::Type{T})::String where {T<:SearchLight.AbstractModel}
  string("FROM ", SearchLight.escape_column_name(SearchLight.table(m), SearchLight.connection()))
end

function SearchLight.to_where_part(w::Vector{SearchLight.SQLWhereEntity}) :: String
  where = isempty(w) ?
          "" :
          string("WHERE ",
                  (string(first(w).condition) == "AND" ? "TRUE " : "FALSE "),
                  join(map(wx -> string(wx), w), " "))

  replace(where, r"WHERE TRUE AND "i => "WHERE ")
end

function SearchLight.to_order_part(m::Type{T}, o::Vector{SearchLight.SQLOrder})::String where {T<:SearchLight.AbstractModel}
  isempty(o) ?
    "" :
    string("ORDER BY ",
            join(map(x -> string((! SearchLight.is_fully_qualified(x.column.value) ?
                                    SearchLight.to_fully_qualified(m, x.column) :
                                    x.column.value), " ", x.direction),
                      o), ", "))
end

function SearchLight.to_group_part(g::Vector{SearchLight.SQLColumn}) :: String
  isempty(g) ?
    "" :
    string(" GROUP BY ", join(map(x -> string(x), g), ", "))
end

function SearchLight.to_limit_part(l::SearchLight.SQLLimit) :: String
  l.value != "ALL" ? string("LIMIT ", string(l)) : ""
end

function SearchLight.to_offset_part(o::Int) :: String
  o != 0 ? string("OFFSET ", string(o)) : ""
end

function SearchLight.to_having_part(h::Vector{SearchLight.SQLWhereEntity}) :: String
  having =  isempty(h) ?
            "" :
            string("HAVING ",
                    (string(first(h).condition) == "AND" ? "TRUE " : "FALSE "),
                    join(map(w -> string(w), h), " "))

  replace(having, r"HAVING TRUE AND "i => "HAVING ")
end

function SearchLight.to_join_part(m::Type{T}, joins::Union{Nothing,Vector{SearchLight.SQLJoin}} = nothing)::String where {T<:SearchLight.AbstractModel}
  joins === nothing && return ""

  join(map(x -> string(x), joins), " ")
end

function Base.rand(m::Type{T}; limit = 1)::Vector{T} where {T<:SearchLight.AbstractModel}
  SearchLight.find(m, SearchLight.SQLQuery(limit = SearchLight.SQLLimit(limit), order = [SearchLight.SQLOrder("random()", raw = true)]))
end

#### MIGRATIONS ####

"""
    create_migrations_table(table_name::String)::Nothing

Runs a SQL DB query that creates the table `table_name` with the structure needed to be used as the DB migrations table.
The table should contain one column, `version`, unique, as a string of maximum 30 chars long.
"""
function SearchLight.Migration.create_migrations_table(table_name::String = SearchLight.config.db_migrations_table_name) :: Nothing
  SearchLight.query(
    "CREATE TABLE `$table_name` (
      `version` varchar(30) NOT NULL DEFAULT '',
      PRIMARY KEY (`version`)
    )")

  @info "Created table $table_name"

  nothing
end

function SearchLight.Migration.create_table(f::Function, name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: Nothing
  SearchLight.query(create_table_sql(f, string(name), options))

  nothing
end

function create_table_sql(f::Function, name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: String
  "CREATE TABLE $name (" * join(f()::Vector{String}, ", ") * ") $options" |> strip
end

function SearchLight.Migration.column(name::Union{String,Symbol}, column_type::Union{String,Symbol}, options::Any = "";
                  default::Any = nothing, limit::Union{Int,Nothing,String} = nothing, not_null::Bool = false) :: String
  "$name $(TYPE_MAPPINGS[column_type] |> string) " *
    (isa(limit, Int) && (! in(TYPE_MAPPINGS[column_type], [:TEXT])) ? "($limit)" : "") *
    (default === nothing ? "" : " DEFAULT $default ") *
    (not_null ? " NOT NULL " : "") *
    " " * string(options)
end

function SearchLight.Migration.column_id(name::Union{String,Symbol} = "id", options::Union{String,Symbol} = "") :: String
  "$name INTEGER PRIMARY KEY $options"
end

function SearchLight.Migration.add_index(table_name::Union{String,Symbol}, column_name::Union{String,Symbol}; name::Union{String,Symbol} = "", unique::Bool = false) :: Nothing
  name = isempty(name) ? SearchLight.index_name(table_name, column_name) : name
  SearchLight.query("CREATE $(unique ? "UNIQUE" : "") INDEX $(name) ON $table_name ($column_name)")

  nothing
end

function SearchLight.Migration.add_column(table_name::Union{String,Symbol}, name::Union{String,Symbol}, column_type::Union{String,Symbol}; default::Union{String,Symbol,Nothing} = nothing, limit::Union{Int,Nothing} = nothing, not_null::Bool = false) :: Nothing
  SearchLight.query("ALTER TABLE $table_name ADD $(SearchLight.Migration.column(name, column_type, default = default, limit = limit, not_null = not_null))")

  nothing
end

function SearchLight.Migration.remove_column(table_name::Union{String,Symbol}, name::Union{String,Symbol}, options::Union{String,Symbol} = "") :: Nothing
  @warn("SQLite does not support DROP COLUMN")
end

function SearchLight.Migration.drop_table(name::Union{String,Symbol}) :: Nothing
  SearchLight.query("DROP TABLE $name")

  nothing
end

function SearchLight.Migration.remove_index(name::Union{String,Symbol}) :: Nothing
  SearchLight.query("DROP INDEX $name")

  nothing
end

#### TRANSACTIONS ####

function SearchLight.Transactions.begin_transaction() :: Nothing
  SearchLight.query("BEGIN")

  nothing
end

function SearchLight.Transactions.commit_transaction() :: Nothing
  SearchLight.query("COMMIT")

  nothing
end

function SearchLight.Transactions.rollback_transaction() :: Nothing
  SearchLight.query("ROLLBACK")

  nothing
end

#### GENERATOR ####

function SearchLight.Generator.FileTemplates.adapter_default_config()
  """
  $(SearchLight.config.app_env):
    adapter:  SQLite
    filename: db/$(SearchLight.config.app_env).sqlite3
  """
end

end


text:
struct BinaryFuseFilter{T <: Union{UInt8, UInt16, UInt32}}
    seed::UInt64
    segment_length::UInt32
    segment_length_mask::UInt32
    segment_count::UInt32
    segment_count_length::UInt32
	array_length::UInt32
    fingerprints::Vector{T}
    n_keys::UInt32
end

function calculate_segment_length(arity::UInt32, n_keys::UInt32)
	# These parameters are very sensitive. Replacing 'floor' by 'round' can
    # substantially affect the construction time.
	if arity == 3
		return UInt32(1) << Int(floor(log(Float64(n_keys)) / log(3.33) + 2.25))
	elseif arity == UInt32(4)
		return UInt32(1) << Int(floor(log(Float64(n_keys)) / log(2.91) - 0.5))
	else
		return UInt32(65536)
    end
end

function calculate_size_factor(arity::UInt32, n_keys::UInt32)
	if arity == 3
		return Float64(max(1.125, 0.875 + 0.25 * log(1000000)/log(Float64(n_keys))))
	elseif arity == 4
		return Float64(max(1.075, 0.77 + 0.305 * log(600000)/log(Float64(n_keys))))
	else
		return Float64(2.0)
	end
end

function get_hash_from_hash(
	hash::UInt64,
	segment_count_length::UInt32,
	segment_length::UInt32,
	segment_length_mask::UInt32)
	hi = binary_fuse_mulhi(hash, UInt64(segment_count_length))
	h0 = UInt32(hi & 0xFFFFFFFF)
	h1 = h0 + segment_length
	h2 = h1 + segment_length
	h1 ⊻= UInt32((hash >> 18) & segment_length_mask)
	h2 ⊻= UInt32(hash & segment_length_mask)
	return h0, h1, h2
end

"""
`BinaryFuseFilter{T}(
	keys::Vector{UInt64};
	seed = UInt64(0x726b2b9d438b9d4d),
	max_iterations = 10) where T <: Union{UInt8, UInt16, UInt32}`

Build a binary fuse filter with T type fingerprints from
keys in `keys`.

## Arguments

`seed` is used as initial seed value,
which is iteratively changed during structure build step.

`max_iterations` restricts how many rounds of build step to try
before termination.

# Examples
```julia-repl
julia> filter = BinaryFuseFilter{UInt8}(UInt64.[1:10])
BinaryFuseFilter with 10 keys and seed 1198702242888554850.

julia> UInt64(1) in filter
true

julia> UInt64(11) in filter
false
```

"""
function BinaryFuseFilter{T}(
	keys::Vector{UInt64};
	seed = UInt64(0x726b2b9d438b9d4d),
	max_iterations = 10) where T <: Union{UInt8, UInt16, UInt32}

	n_keys = UInt32(length(keys))
	if n_keys == 0
		throw("No keys for construction.")
	end

	arity = UInt32(3)
	segment_length = calculate_segment_length(arity, n_keys)
	if (segment_length > 262144)
    	segment_length = UInt32(262144)
    end

	segment_length_mask = UInt32(segment_length - 1)
	size_factor = calculate_size_factor(arity, n_keys)
	capacity = UInt32(round(Float64(n_keys) * size_factor))
	segment_count = UInt32(ceil((capacity + segment_length - 1) / segment_length - (arity - 1)))
	array_length = (segment_count + arity - 1) * segment_length
	segment_count = UInt32(ceil((array_length + segment_length - 1)/segment_length))
	if (segment_count <= (arity - 1))
		segment_count = 1
	else
		segment_count = segment_count - (arity - 1)
    end
	array_length = (segment_count + arity - 1) * segment_length

	segment_count_length = UInt32(segment_count * segment_length)
	fingerprints = zeros(T, array_length)

	# Adding elements
	rngcounter = seed
	seed, rngcounter = splitmix64(rngcounter)
	capacity = UInt32(array_length)

	alone = zeros(UInt32, capacity)
	# the lowest 2 bits are the h index (0, 1, or 2)
	# so we only have 6 bits for counting;
	# but that's sufficient
	t2count = zeros(UInt8, capacity)
	reverseH = zeros(UInt8, n_keys)
	t2hash = zeros(UInt64, capacity)
	reverseOrder = zeros(UInt64, n_keys + 1)
	reverseOrder[n_keys + 1] = 1 # uneccessary item serving as a blocker

	# the array h0, h1, h2, h0, h1, h2
	h012 = zeros(UInt32, 5)
	iterations = 0

	blockBits = 1
	while UInt32(1 << blockBits) < segment_count
		blockBits += 1
	end
	block = 1 << blockBits

	while true
		iterations += 1
		if iterations > max_iterations
			throw("Too many iterations, you probably have duplicate keys.")
		end

		startPos = zeros(UInt32, block)
		for i in 1:length(startPos)
			# important: we do not want i * n_keys to overflow!!!
			# in java we have (int) ((long) i * size / block)
			startPos[i] = UInt32(((UInt64(i - 1) * UInt64(n_keys)) >> blockBits))
		end

		# counting sort to the h0 region
		maskblock = block - 1
		for key in keys
			hash = mixsplit(key, seed)
			# this is java >>> shifted values are zeros
			# java (int) (hash >>> (64 - blockBits))
			segment_index = hash >> (64 - blockBits)
			while reverseOrder[startPos[segment_index + 1] + 1] != 0
				segment_index += 1
				segment_index &= maskblock
			end
			reverseOrder[startPos[segment_index + 1] + 1] = hash
			startPos[segment_index + 1] += 1
		end

		# We can then scan through the partially sorted keys and update a temporary array of counters—
        # such an array tells us how many keys map to a given location. Because the keys are sorted
        # by region, a forward pass through the keys tends to access the counters in a forward manner,
        # thus reducing the number of cache misses compared to a random approach.
		error = 0
		for i in UInt32.(1:n_keys)
			hash = reverseOrder[i]
			index1, index2, index3 = get_hash_from_hash(
				hash, segment_count_length, segment_length, segment_length_mask)
			t2count[index1 + 1] += 4
			t2hash[index1 + 1] ⊻= hash
			t2count[index2 + 1] += 4
			t2count[index2 + 1] ⊻= 1
			t2hash[index2 + 1] ⊻= hash
			t2count[index3 + 1] += 4
			t2count[index3 + 1] ⊻= 2
			t2hash[index3 + 1] ⊻= hash
			if t2count[index1 + 1] < 4 # no clue what is being checked here...
				error = 1
			end
			if t2count[index2 + 1] < 4
				error = 1
			end
			if t2count[index3 + 1] < 4
				error = 1
			end
		end
		if error == 1
			continue
		end

		Qsize = 0
		# We scan the array of counters to identify the locations corresponding to a single set entry.
        # The entries are added to a stack. Since we scan forward, the later entries in the stack tend to
        # correspond to later locations.
		for i in UInt32.(1:capacity)
			alone[Qsize + 1] = i - 1
			if (t2count[i] >> 2) == 1
				Qsize += 1
			end
		end
		stacksize = UInt32(0)
		while Qsize > 0
			Qsize -= 1
			index = alone[Qsize + 1]
			if (t2count[index + 1] >> 2) == 1
				hash = t2hash[index + 1]
				found = t2count[index + 1] & 3
				reverseH[stacksize + 1] = found
				reverseOrder[stacksize + 1] = hash
				stacksize += 1

				index1, index2, index3 = get_hash_from_hash(
					hash, segment_count_length, segment_length, segment_length_mask)

				h012[1] = index2
				h012[2] = index3
				h012[3] = index1
				h012[4] = h012[1]

				other_index1 = h012[found + 1]
				alone[Qsize + 1] = other_index1
				if (t2count[other_index1 + 1] >> 2) == 2
					Qsize += 1
				end
				t2count[other_index1 + 1] -= 4
				t2count[other_index1 + 1] ⊻= mod3(found + 1)
				t2hash[other_index1 + 1] ⊻= hash

				other_index2 = h012[found + 2]
				alone[Qsize + 1] = other_index2
				if (t2count[other_index2 + 1] >> 2) == 2
					Qsize += 1
				end
				t2count[other_index2 + 1] -= 4
				t2count[other_index2 + 1] ⊻= mod3(found + 2)
				t2hash[other_index2 + 1] ⊻= hash
			end
		end

		if stacksize == n_keys
			break
		end
		for i in 1:n_keys
			reverseOrder[i] = 0
		end
		for i in 1:capacity
			t2count[i] = 0
			t2hash[i] = 0
		end
		rngcounter, seed = splitmix64(rngcounter)
	end

	# Finally, starting from the latter stack, we can construct the binary fuse filter, by making sure
    # that the bitwise exclusive or of the three locations map to the fingerprint of the matched
    # entries. By construction, we tend to go from locations that are at the beginning of the array,
    # working toward the end of the array.
	for i in n_keys:-1:1
		# the hash of the key we insert next
		hash = reverseOrder[i]
		xor2 = fingerprint(T, hash)
		index1, index2, index3 = get_hash_from_hash(
			hash, segment_count_length, segment_length, segment_length_mask)
		found = reverseH[i]
		h012[1] = index1
		h012[2] = index2
		h012[3] = index3
		h012[4] = h012[1]
		h012[5] = h012[2]
		fingerprints[h012[found + 1] + 1] = xor2 ⊻ fingerprints[h012[found + 2] + 1] ⊻ fingerprints[h012[found + 3] + 1]
	end

	return BinaryFuseFilter{T}(
		seed,
		segment_length, segment_length_mask,
    	segment_count, segment_count_length,
		array_length,
    	fingerprints, n_keys)
end

function Base.summary(io::IO, filter::BinaryFuseFilter)
    print(io, "BinaryFuseFilter with $(filter.n_keys) keys and seed $(filter.seed).")
end
Base.show(io::IO, filter::BinaryFuseFilter) = summary(io, filter)

function Base.in(key, filter::BinaryFuseFilter)
    hash = mixsplit(key, filter.seed)
    f = fingerprint(typeof(filter.fingerprints[1]), hash)
	h0, h1, h2 = get_hash_from_hash(
		hash, filter.segment_count_length, filter.segment_length, filter.segment_length_mask)
	f ⊻= (filter.fingerprints[h0 + 1] ⊻ filter.fingerprints[h1 + 1] ⊻ filter.fingerprints[h2 + 1])
	return f == 0
end


text:
"""
    flatfield(img)

Compute the flatfield image from `img` by calculating the median value
along the 3rd dimension which normally corresponds to different
positions. This way the contribution of the optical path can be
determined. Requires an `img` with labeled dimensions `x` and `y`, order
does not matter.
"""
function flatfield(img::AxisArray{T, 3}) where {T}
    ysize, xsize = size(img, Axis{:y}), size(img, Axis{:x})
    med = Array{T}(undef, ysize, xsize)

    R = CartesianIndices((xsize, ysize))
    for i in R
        med[i] = median(img[Axis{:y}(i[1]), Axis{:x}(i[2])])
    end
    med
end

@traitfn function flatfield_correct!(output,
                                     img::AA,
                                     flatfield::AbstractArray{T1, 2},
                                     darkfield::AbstractArray{T2, 2}) where {AA <: AxisArray, T1, T2; HasTimeAxis{AA}}
    for timepoint in timeaxis(img)
        tp = Axis{:time}(timepoint)
        _img = @view img[tp]
        _output = @view output[tp]
        flatfield_correct!(_output, _img, flatfield, darkfield)
    end
end

@traitfn function flatfield_correct!(output,
                                     img::AA,
                                     flatfield::AbstractArray{T1, 2},
                                     darkfield::AbstractArray{T2, 2}) where {AA <: AxisArray, T1, T2; !HasTimeAxis{AA}}
    R = CartesianIndices(Base.axes(img))
    for I in R
        output[I] = (float(img[I]) - darkfield[I]) / (flatfield[I] - darkfield[I])
    end
end

"""
    flatfield_correct(img, flatfield, darkfield)

Computes the flatfield for the given image `img` by dividing the image by the
flatfield after subtracting the darkfield from both[^1]

!!! note

    This function assumes that the flatfield image has *not* had the darkfield
    image subtracted yet.

[^1]: http://nic.ucsf.edu/resources/how-to-acquire-flat-field-correction-images/
"""
function flatfield_correct(img::AxisArray, flatfield::AbstractArray{T1, 2}, darkfield::AbstractArray{T2, 2}) where {T1, T2}
    output = similar(img, Gray{Float64})
    flatfield_correct!(output, img, flatfield, darkfield)
    output
end

flatfield_correct(img::ImageMeta, flatfield, darkfield) = copyproperties(img, flatfield_correct(img.data, flatfield, darkfield))

light_source_contrib(img::AxisArray, seeds::BitArray{3}; dist::Int=30, h::Float64=0.9) =
    light_source_contrib(img, AxisArray(seeds, img.axes); dist=dist, h=h)

"""
    light_source_contrib(img, seeds, dist, h)

Computes the contribution of the light source fluctuations for each time point by first
identifying background pixels by finding pixels that are more than `dist` away from true
values in `seeds`. A kernel density estimate is then fit to the background pixels to get
a continuous distribution and then [`find_peak_center`](@ref) is used to identify the
"average" response of a background pixel.

!!! note

    The distribution of background pixels often ends up being quite non-normal, likely due
    to misidentification. Empirically, I've found the finding the center of the kernel
    density estimate to be more robust than simpler statistics like mean, median, or even
    the maximum value of the KDE.

### Rationale

Certain light sources, especially arc lamps, exhibit substantial variance in
total brightness in time. Additionally, the background signal is dependent on the total
light delivered and we can use that to identify the arc lamp wander and subtract it
from the whole field of view at each timestep. This gives us more stable total
fluorescence values over time.

### Example:

```julia
a = img[Axis{:position}(2), Axis{:channel}(:EPI_mNG)].data
b = img[Axis{:position}(2), Axis{:channel}(:EPI_BFP)].data .> 0.01
SegmentationTools.light_source_contrib(a, b)
```
"""
function light_source_contrib(img::AxisArray{T1, 3},
                              seeds::AxisArray{T2, 3};
                              dist::Int = 30,
                              h::Float64 = 0.9) where {T1, T2 <: Bool}
    n = length(timeaxis(img))
    (n != length(timeaxis(seeds))) && error("The time dimensions of both arrays need to the same")
    bkg_means = Array{Float64}(undef, n)
    for i in 1:n
        seed_slice = seeds[Axis{:time}(i)]
        bkg_mask = distance_transform(feature_transform(seed_slice)) .> dist
        # fit a kernel density estimate to the background pixels
        fitted_kde = kde(vec(Float64.(img[Axis{:time}(i)][bkg_mask])))
        # find the center of the KDE of the background peak
        bkg_means[i] = find_peak_center(collect(fitted_kde.x), fitted_kde.density; h=h)
    end
    bkg_means
end

"""
    find_peak_center(x, y, h)

Given two equal length arrays, `x`, and `y`, computes the half width
at `h` relative to the maximum of the peak in `y`, i.e. for `h=0.5`, this
function returns the value of `x` that corresponds to the half-width
half-max. This is a robust method for identifying the center of a peak.
"""
function find_peak_center(x::Vector{Float64}, y::Vector{Float64}; h::Float64=0.75)
    (length(x) != length(y)) && throw(DimensionMismatch("x and y have to have equal lengths"))
    # value at which we cut the peak
    cutoff = maximum(y)*h

    # find the indices where we "enter" and "exit" the peak
    transitions = findall(i-> i != 0, diff(y .> cutoff))
    (length(transitions) != 2) && error("Dependent variable is not unimodal")

    # translate from indices to the original units of `slice`
    low = x[transitions[1] + 1]
    high = x[transitions[2]]
    low + (high - low) / 2
end

"""
    colorize(img, scheme)

Given a grayscale image `img`, apply a colorscheme `scheme` to the image lazily.
"""
function colorize(img::AbstractArray{Gray{T}, N}; scheme::Symbol=:magma) where {T, N}
    minval, maxval = gray.(extrema(img))
    cscheme = getfield(ColorSchemes, scheme)
    mappedarray(x->RGB{T}(get(cscheme, gray(x), (minval, maxval))), img)
end


text:
using GenericInstruments
GI = GenericInstruments
####################################
resmgr = GI.viOpenDefaultRM()
# Instantiate obj
scope1   = GI.SCOPE.INSTR(:HDO6054A, "USB0::0x05FF::0x1023::4066N51752::INSTR")
GI.SCOPE.connect!(resmgr, scope1)
# Most commands needs Dicts...
# IDN
GI.SCOPE.get_idn(scope1)
# Fetch waveform from channel 1
y, initial_x, x_increm = GI.SCOPE.fetch_wfm(scope1;ch="1")
# Default Setup
GI.SCOPE.default_setup(scope1)
# Vertical range
GI.SCOPE.set_vrange(scope1;ch="1",vrang=12)
# Vertical scale
GI.SCOPE.set_vscale(scope1;ch="1",volt=0.2)
# Vertical offset
GI.SCOPE.set_voffs(scope1;ch="1",voffs=2.2)
# Channel coupling
GI.SCOPE.set_coupling(scope1;ch="1",cpl="ac")
# Channel enabled
GI.SCOPE.set_state(scope1;ch="1",st="off")
# Probe attenuation
GI.SCOPE.set_atten(scope1;ch="1",att="10x")#not working, check VBS
# Screen Gridmode
GI.SCOPE.set_gridmode(scope1;gridmode="single")
# Nr of averages
GI.SCOPE.set_nr_avg(scope1;ch="1",navg=10)
# EnhanceRes
GI.SCOPE.set_eres(scope1;ch="1",bits="0.5b")
# Degauss
GI.SCOPE.degauss(scope1;ch="4",probe="CP030A")
# Clear sweeps
GI.SCOPE.clear_sweeps(scope1) # ??
# Nr of acquisition points
GI.SCOPE.set_nr_pts(scope1;npts=1e6)
# Sampling rate
GI.SCOPE.set_srate(scope1;sr=50e6)
# Acquisition duration
GI.SCOPE.set_hduration(scope1;d=1e-3)
# Horizontal offset origin
GI.SCOPE.set_hoffs_div(scope1;href=1)
# Set Measurement
GI.SCOPE.set_meas(scope1;ch="1",fct="max",par=1,vw=1)#not working, check VBS
# Get Measurement
GI.SCOPE.get_meas(scope1;par=1)
# Set trigger source
GI.SCOPE.set_trig_src(scope1;src="1")
# Set trigger type
GI.SCOPE.set_trig_type(scope1;tp="edge")
# Set trigger mode
GI.SCOPE.set_trig_mode(scope1;mode="single")
# Set trigger slope
GI.SCOPE.set_trig_slope(scope1;slp="pos")
# Set trigger level
GI.SCOPE.set_trig_lev(scope1;lev=1.2) #not working?
# Set trigger holdoff
GI.SCOPE.set_trig_holdoff(scope1;hf=1e-3)#not working?
# Trigger delay time
GI.SCOPE.set_hoffs_t(scope1;href=11e-3)
# Set trigger coupling
GI.SCOPE.set_trig_cpl(scope1;src="1",cpl="dc")
# BW limit
GI.SCOPE.set_bw_lim(scope1;ch="1",bw="full")
# Close
GI.SCOPE.disconnect!(scope1)
# pending...


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 = "MAGEMin"
version = v"1.1.1"

# Collection of sources required to complete build
sources = [GitSource("https://github.com/ComputationalThermodynamics/MAGEMin",
                    "72198b66c94a7026497824b460096d7e034cecca")
        ]

# Bash recipe for building across all platforms
script = raw"""
cd MAGEMin*
if [[ "${target}" == *-mingw* ]]; then
    MPI_LIBS="-lmsmpi"
else
    MPI_LIBS="-lmpi"
fi

CCFLAGS="-O3 -g -fPIC -std=c99"
LIBS="-L${libdir} -lm -lopenblas -lnlopt ${MPI_LIBS}"
INC="-I${includedir}"

# Compile library:
make -j${nproc} CC="${CC}" CCFLAGS="${CCFLAGS}" LIBS="${LIBS}" INC="${INC}" lib

# compile binary
make -j${nproc} CC="${CC}" CCFLAGS="${CCFLAGS}" LIBS="${LIBS}" INC="${INC}" all

install -Dvm 755 libMAGEMin.dylib "${libdir}/libMAGEMin.${dlext}"
install -Dvm 755 MAGEMin* "${bindir}/MAGEMin${exeext}"

# store files
install -vm 644 src/*.h "${includedir}"

install_license LICENSE
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libMAGEMin", :libMAGEMin)
    ExecutableProduct("MAGEMin", :MAGEMin)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="MPICH_jll", uuid="7cb0a576-ebde-5e09-9194-50597f1243b4"); platforms=filter(!Sys.iswindows, platforms))
    Dependency(PackageSpec(name="MicrosoftMPI_jll", uuid="9237b28f-5490-5468-be7b-bb81f5f5e6cf"); platforms=filter(Sys.iswindows, platforms))
    Dependency(PackageSpec(name="NLopt_jll", uuid="079eb43e-fd8e-5478-9966-2cf3e3edb778"))
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"6")


text:
# indicator of the zero cone

export IndZero

"""
**Indicator of the zero cone**

    IndZero()

Returns the indicator function of the set containing the origin, the "zero cone".
"""
struct IndZero <: ProximableFunction end

is_separable(f::IndZero) = true
is_convex(f::IndZero) = true
is_singleton(f::IndZero) = true
is_cone(f::IndZero) = true
is_affine(f::IndZero) = true

function (f::IndZero)(x::AbstractArray{C}) where {R <: Real, C <: Union{R, Complex{R}}}
    for k in eachindex(x)
        if x[k] != C(0)
            return R(Inf)
        end
    end
    return R(0)
end

function prox!(y::AbstractArray{C}, f::IndZero, x::AbstractArray{C}, gamma=R(1)) where {R <: Real, C <: Union{R, Complex{R}}}
    for k in eachindex(x)
        y[k] = C(0)
    end
    return R(0)
end

fun_name(f::IndZero) = "indicator of the zero cone"
fun_dom(f::IndZero) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndZero) = "x ↦ 0 if all(x = 0), +∞ otherwise"
fun_params(f::IndZero) = "none"

function prox_naive(f::IndZero, x::AbstractArray{C}, gamma=R(1)) where {R <: Real, C <: Union{R, Complex{R}}}
    return zero(x), R(0)
end


text:
const ps = SmartGameFormat.Parser

@testset "public interface" begin
    # nothing should be exported
    @test_throws UndefVarError Parser
    @test ParseError === Base.ParseError
    @test parse === Base.parse
    @test tryparse === Base.tryparse
end

# unsupported FF
@testset "FF and other root properties" begin
    @test parse_sgf("(; FF[4])") == [SGFGameTree(SGFNode(:FF=>4))]
    @test parse_sgf("(; CA[UTF-8])") == [SGFGameTree(SGFNode(:CA=>"UTF-8"))]
    @test parse_sgf("(; FF[4])(; CA[UTF-8])") == [SGFGameTree(SGFNode(:FF=>4)), SGFGameTree(SGFNode(:CA=>"UTF-8"))]
    # warning for older versions
    @test @test_warn("FF",parse_sgf("(; FF[3])")) == [SGFGameTree(SGFNode(:FF=>3))]
    # warning for unsupported charsets
    @test @test_warn("unsupported",parse_sgf("(; CA[ABC])")) == [SGFGameTree(SGFNode(:CA=>"ABC"))]
    # illegal version numbers
    @test_throws ps.ParseError parse_sgf("(; FF[5])")
    @test_throws ps.ParseError parse_sgf("(; FF[0])")
    # root property in non root node
    @test_throws ps.ParseError parse_sgf("(; ; FF[4])")
    @test_throws ps.ParseError parse_sgf("(; ; CA[UTF-8])")
end

# illegal value types (not numbers)
@testset "non-numbers where numbers are expected" begin
    # basically checks that a `ps.ParseError` is thrown
    # instead of a `Base.ParseError`
    @test_throws ps.ParseError parse_sgf("(; FF[A])")
    @test_throws ps.ParseError parse_sgf("(; KM[A])")
end

@testset "minimal SGF specifications" begin
    @test_throws ps.ParseError parse_sgf("()")
    @test_throws ps.ParseError parse_sgf("(FF[A])")
    @test_throws ps.ParseError parse_sgf("(;)()")
    @test parse_sgf("(;)") == [SGFGameTree(SGFNode())]
    @test parse_sgf("(\n;)") == [SGFGameTree(SGFNode())]
    @test parse_sgf("( ;)") == [SGFGameTree(SGFNode())]
    @test parse_sgf("(;)(;)") == [SGFGameTree(SGFNode()),SGFGameTree(SGFNode())]

end


text:
using ChainRulesCore
const CRC = ChainRulesCore
using ChainRulesTestUtils

import MathOptSetDistances
const MOD = MathOptSetDistances
const MOI = MathOptSetDistances.MOI
using FiniteDifferences
using Test
using Random
using LinearAlgebra

# type piracy here to use FiniteDiff
function FiniteDifferences.to_vec(s::S) where {S <: Union{MOI.EqualTo, MOI.LessThan, MOI.GreaterThan}}
    function set_from_vec(v)
        return S(v[1])
    end
    return [MOI.constant(s)], set_from_vec
end

# avoid random finite diff fails because of rounding
Random.seed!(42)

"""
Used to test all distances where the set does not take differentiable parameters.
"""
function test_rrule_analytical(x, s; distance = MOD.DefaultDistance(), ntrials = 10, rng=Random.GLOBAL_RNG, atol=1e-4,rtol=1e-4, test_fdiff=true)
    y = MOD.projection_on_set(distance, x, s)
    dΠ = MOD.projection_gradient_on_set(distance, x, s)
    (yprimal, pullback) = CRC.rrule(MOD.projection_on_set, distance, x, s)
    @test yprimal ≈ y
    for _ in 1:ntrials
        xb = randn(rng, length(x))
        yb = randn(rng, length(x))
        (_, _, Δx, _) = pullback(yb)
        @test Δx ≈ dΠ' * yb
        if test_fdiff
            ChainRulesTestUtils.test_rrule(
                MOD.projection_on_set,
                ChainRulesTestUtils.PrimalAndTangent(distance, CRC.NoTangent()),
                ChainRulesTestUtils.PrimalAndTangent(x, xb),
                ChainRulesTestUtils.PrimalAndTangent(s, CRC.NoTangent()),
                atol=atol,rtol=rtol,
                output_tangent=yb,
            )
        end
    end
end

@testset "rrules multivariate" begin
    for n in (1, 2, 10)
        x = randn(n)
        s = MOI.Reals(n)
        test_rrule_analytical(x, s, atol=1e-5, rtol=1e-5)
        s = MOI.Zeros(n)
        test_rrule_analytical(x, s, atol=1e-5, rtol=1e-5, test_fdiff=false)
        @testset "Orthant $s" for s in (MOI.Nonpositives(n), MOI.Nonnegatives(n))
            test_rrule_analytical(x, s)
        end
        @testset "SOC $n" begin
            s = MOI.SecondOrderCone(n+1)
            for _ in 1:10
                x = 10 * randn(n+1)
                if norm(x[2:end]) ≈ x[1]
                    if rand() > 0.8
                        x[1] *= 2
                    else
                        x[1] /= 2
                    end
                end
                test_rrule_analytical(x, s)
            end
        end
    end
end

@testset "rrules univariate" begin
    # note: we do not use the generic test function since sets have derivatives
    @testset "Scalar $ST" for ST in (MOI.EqualTo, MOI.LessThan, MOI.GreaterThan)
        for _ in 1:10
            s = ST(10 * randn())
            x = 10 * randn()
            # avoid non-differentiable points
            if isapprox(x, MOI.constant(s), atol=1e-5)
                δ = 2rand()
                if rand(Bool)
                    x += δ
                else
                    x -= δ
                end
            end
            dΠ = MOD.projection_gradient_on_set(MOD.DefaultDistance(), x, s)
            y = MOD.projection_on_set(MOD.DefaultDistance(), x, s)
            (yprimal, pullback) = CRC.rrule(MOD.projection_on_set, MOD.DefaultDistance(), x, s)
            @test yprimal ≈ y
            for _ in 1:5
                xb = ChainRulesTestUtils.rand_tangent(x)
                yb = ChainRulesTestUtils.rand_tangent(y)
                sb = ChainRulesTestUtils.rand_tangent(s)
                ChainRulesTestUtils.test_rrule(
                    MOD.projection_on_set,
                    ChainRulesTestUtils.PrimalAndTangent(MOD.DefaultDistance(), CRC.NoTangent()),
                    ChainRulesTestUtils.PrimalAndTangent(x, xb),
                    ChainRulesTestUtils.PrimalAndTangent(s, sb),
                    output_tangent=yb,
                    atol=1e-4,
                )
                (_, _, Δx, _) = pullback(yb)
                ChainRulesTestUtils.test_approx(Δx, dΠ' * yb)
            end
        end
    end
end

@testset "frule" begin
    for n in (1, 2, 10)
        @testset "$s" for s in (MOI.Nonnegatives(n), MOI.Nonpositives(n))
            for _ in 1:10
                v = 50 * safe_randn(n)
                for _ in 1:5
                    Δv = ChainRulesTestUtils.rand_tangent(v)
                    ChainRulesTestUtils.test_frule(
                        MOD.projection_on_set,
                        ChainRulesTestUtils.PrimalAndTangent(DD, CRC.NoTangent()),
                        ChainRulesTestUtils.PrimalAndTangent(v, Δv),
                        ChainRulesTestUtils.PrimalAndTangent(s, CRC.NoTangent()),
                        atol=1e-5,
                    )
                end
            end
        end

        s = MOI.PositiveSemidefiniteConeTriangle(n)
        @testset "$s" begin
            for _ in 1:5
                L = 3 * tril(rand(n, n))
                M = L * L'
                v0 = MOD.vec_symm(M)
                v = Vector{Float64}(undef, length(v0))
                Π = Vector{Float64}(undef, length(v0))
                Δv = Vector{Float64}(undef, length(v0))
                dΠ = Matrix{Float64}(undef, length(v0), length(v0))
                @testset "Positive and negative definite" begin
                    for _ in 1:3
                        Δv .= ChainRulesTestUtils.rand_tangent(v)
                        v .= v0
                        (vproj, Δvproj) = CRC.frule((CRC.NoTangent(), CRC.NoTangent(), Δv, CRC.NoTangent()), MOD.projection_on_set, DD, v, s)
                        @test Δvproj ≈ Δv
                        @test vproj ≈ v
                        v .= -v0
                        dΠ .= MOD.projection_gradient_on_set(DD, v, s)
                        Π .= MOD.projection_on_set(DD, v, s)
                        (vproj, Δvproj) = CRC.frule((CRC.NoTangent(), CRC.NoTangent(), Δv, CRC.NoTangent()), MOD.projection_on_set, DD, v, s)
                        @test dΠ * Δv ≈ Δvproj
                        @test vproj ≈ Π
                    end
                end
            end
        end
    end
    @testset "Indefinite matrix" begin
        s = MOI.PositiveSemidefiniteConeTriangle(2)
        A = Matrix{Float64}(undef, 2, 2)
        Π = Matrix{Float64}(undef, 2, 2)
        v = Vector{Float64}(undef, 3)
        DΠ = Matrix{Float64}(undef, 3, 3)
        Xd = Matrix{Float64}(undef, 2, 2)
        Q = [
            1 0
            0 -1
        ]
        Qi = Q'
        B = [
            0 1/2
            1/2 1
        ]
        for _ in 1:10
            # scale factor
            f = 20 * rand() + 5
            A .= [
                -f 0
                0 f
            ]
            Λ = Diagonal([-f, f])
            Λp = Diagonal([0, f])
            v .= MOD.vec_symm(A)
            vproj = MOD.projection_on_set(DD, v, s)
            Π .= MOD.unvec_symm(vproj, 2)
            @test Π ≈ Q * Λp * Qi
            DΠ .= MOD.projection_gradient_on_set(DD, v, s)
            for _ in 1:20
                Xd .= ChainRulesTestUtils.rand_tangent(Π)
                xd = MOD.vec_symm(Xd)
                dir_deriv_theo = MOD.vec_symm(
                    Q * (B .* (Q' * Xd * Q)) * Q'
                )
                @test DΠ * xd ≈ dir_deriv_theo
                (vproj_frule, Δvproj) = CRC.frule((CRC.NoTangent(), CRC.NoTangent(), xd, CRC.NoTangent()), MOD.projection_on_set, DD, v, s)
                @test DΠ * xd ≈ Δvproj
                @test vproj ≈ vproj_frule
            end
        end
    end
    @testset "Exp Cone" begin
        function det_case_exp_cone(v; dual=false)
            v = dual ? -v : v
            if MOD.distance_to_set(DD, v, MOI.ExponentialCone()) < 1e-8
                return 1
            elseif MOD.distance_to_set(DD, -v, MOI.DualExponentialCone()) < 1e-8
                return 2
            elseif v[1] <= 0 && v[2] <= 0 #TODO: threshold here??
                return 3
            else
                return 4
            end
        end
        s = MOI.ExponentialCone()
        case_p = zeros(4)
        tol = 1e-6
        for ii in 1:100
            v = 5*randn(3)
            vproj0 = MOD.projection_on_set(DD, v, s)
            @testset "Primal Cone" begin
                case_p[det_case_exp_cone(v; dual=false)] += 1
                dΠ = MOD.projection_gradient_on_set(DD, v, s)
                grad_fdm1 = FiniteDifferences.jacobian(ffdm, x -> MOD.projection_on_set(DD, x, s), v)[1]'
                grad_fdm2 = FiniteDifferences.jacobian(bfdm, x -> MOD.projection_on_set(DD, x, s), v)[1]'
                @test size(grad_fdm1) == size(grad_fdm2) == size(dΠ)
                @test ≈(dΠ, grad_fdm1,atol=tol) || ≈(dΠ, grad_fdm2, atol=tol)
                for _ in 1:20
                    Δv = 5*randn(3)
                    (vproj, Δvproj) = CRC.frule((CRC.NoTangent(), CRC.NoTangent(), Δv, CRC.NoTangent()), MOD.projection_on_set, DD, v, s)
                    @test dΠ * Δv ≈ Δvproj atol=tol
                    @test vproj ≈ vproj0
                end
            end
        end
        @test all(case_p .> 0)
    end
end


text:
"""
    platform_check_with_freq()
---
- Calculate frequencies
- Calculate correlation matrix of the 11 platforms
"""
function platform_check_with_freq()
    platforms = ["dutch-d1",
                 "dutch-d2",
                 "dutch-d3",
                 "dutch-v2",
                 "dutch-v3",
                 "dutch-v7",
                 "german-v2",
                 "german-v3",
                 "norge-v1",
                 "norge-v2",
                 "norge-v7"]
    title("Correlation between 11 platforms")
    # do_all_shared_alleles_have_the_same_alternatives(platforms)

    subtitle("Using raw data")
    # raw_frq(platforms)

    subtitle("Using imputed data")
    imp_frq(platforms)
end

"""
    function imp_frq(platforms)
---
Split country data on ID of different platforms.  Then calculate new frequencies,
and new correlation like with their raw genotypes.
"""

function imp_frq(platforms)
    raw = joinpath(work_dir, "data/genotypes/step-5.plk")
    fra = joinpath(work_dir, "data/genotypes/step-8.plk")
    til = joinpath(work_dir, "notebooks/frq")
    tmp = joinpath(work_dir, "tmp")
    empty_dir(tmp)
    countries = ["dutch", "german", "norge"]

    item("Merge 3-country data")
    open(joinpath(tmp, "country.lst"), "w") do io
        for country in countries
            write(io, "$fra/$country\n")
        end
    end
    merge_beds(joinpath(tmp, "country.lst"),
               joinpath(tmp, "all")
               )
    done()

    item("Extract all ID from imputed and QC-ed country data")
    ID = begin
        fam = CSV.read(joinpath(tmp, "all.fam"), header = [:a, :id, :b, :c, :d, :c])
        Set(fam.id)
    end
    done()

    item("Platform specific frequencies")
    for pf in platforms
        id = begin
            fam = CSV.read(joinpath(raw, "$pf.fam"), header = [:a, :id, :b, :c, :d, :c])
            Set(fam.id)
        end
        open(joinpath(tmp, "keep.id"), "w") do io
            for i in id
                write(io, "0 $i\n")
            end
        end
        plink_keep_id(joinpath(tmp, "all"),
                      joinpath(tmp, "keep.id"),
                      joinpath(tmp, "sub"))
        allele_maf_stats(joinpath(tmp, "sub"),
                         joinpath(tmp, pf))
    end
    done()

    item("Correlation of frequencies")
    freq = begin
        DF = Any[]
        for pf in platforms
            df = CSV.read(joinpath(tmp, "$pf.frq"), delim = ' ', ignorerepeated = true)
            push!(DF, select(df, :SNP, :MAF))
        end
        innerjoin(DF..., on = :SNP, makeunique = true)
    end
    corrs = DataFrame(cor(Matrix(select(freq, r"MAF"))))
    names = ["SNP", "dd1", "dd2", "dd3", "dv2", "dv3", "dv7", "gv2", "gv3", "nv1", "nv2", "nv7"]
    rename!(freq, names)
    CSV.write(joinpath(til, "imp-frq.csv"), freq)
    rename!(corrs, names[2:end])
    CSV.write(joinpath(til, "imp-cor.csv"), corrs)
    done()
end

"""
    raw_frq(fra, platforms, out)
---
Calculate platform specific allele frequencies.  Then compare them on shared loci.
Calculate before and after imputation.
"""
function raw_frq(platforms)
    item("Platform specific allele frequencies")
    cd(work_dir)
    fra = joinpath(work_dir, "data/genotypes/step-5.plk")
    tmp = joinpath(work_dir, "tmp")
    empty_dir(tmp)
    for pf in platforms
        allele_maf_stats(joinpath(fra, pf), joinpath(tmp, pf))
    end
    done()

    item("Calculate correlation matrix of the 11 platforms")
    til = joinpath(work_dir, "notebooks/frq")
    isdir(til) || mkdir(til)

    shrfreq = begin
        DF = Any[]
        for pf in platforms
            df = CSV.read(joinpath(tmp, "$pf.frq"), delim = ' ', ignorerepeated = true)
            push!(DF, select(df, :SNP, r"MAF"))
        end
        innerjoin(DF..., on = :SNP, makeunique = true)
    end                         # frequencies of shared SNP
    pfcor = DataFrame(cor(Matrix(select(shrfreq, r"MAF"))))

    names = ["SNP", "dd1", "dd2", "dd3", "dv2", "dv3", "dv7", "gv2", "gv3", "nv1", "nv2", "nv7"]
    rename!(shrfreq, names)
    CSV.write(joinpath(til, "raw-frq.csv"), shrfreq)

    rename!(pfcor, names[2:end])
    CSV.write(joinpath(til, "raw-cor.csv"), pfcor)
    done()
end

"""
    do_all_shared_alleles_have_the_same_alternatives(platforms)
---
This function was originally to swap alleles, as order of the alternatives
affects the genotype coding.  But `plink` calculates MAF.  So, swapping the
alternative alleles won't affect the MAF.  The swapping was not successful in
the first place.  I took the apportunity to check if all platforms are using
`Top` results.
"""
function do_all_shared_alleles_have_the_same_alternatives(platforms)
    subtitle("Using raw shared SNP")
    item("Check if they are all of Top results")
    fra = joinpath(work_dir, "data/genotypes/step-5.plk")

    names = ["CHR", "SNP", "Gap", "BP", "A", "B"]
    message("Shared loci and their alternative alleles")
    snpmap = begin
        df = Any[]
        for pf in platforms
            dd = CSV.read(joinpath(fra, "$pf.bim"), header = names)
            select!(dd, [:SNP, :A, :B])
            push!(df, dd)
        end
        innerjoin(df..., on = :SNP, makeunique = true)
    end
    # To write the SNP for genotype subset.
    # CSV.write(joinpath(tmp, "keep.snp"),
    #           select(snpmap, :SNP),
    #           writeheader = false)
    #
    # bed_snp_subset(joinpath(fra, platforms[1]),
    #                joinpath(tmp, "keep.snp"),
    #                joinpath(tmp, platforms[1])
    #                )

    message("The rest platforms, need some swapping")
    ref = Matrix(select(snpmap, [:A, :B]))
    i = 1
    for pf in platforms[2:end]  # Not successful
        tst = Matrix(select(snpmap, ["A_$i", "B_$i"]))
        i += 1
        same = true
        for k in 1:size(ref)[1]
            if ref[k, 1] == tst[k, 1]
                if ref[k, 2] ≠ tst[k, 2]
                    same = false
                end
            else
                if ref[k, 1] ≠ tst[k, 2] || ref[k, 2] ≠ tst[k, 1]
                    same = false
                end
            end
        end
        if !same
            printly("Platform ", platforms[i], " is different.")
        end
    end
    done()
end

"""
    check_abnormal_elements()
---
# Issues
- dv_2 has very high diagonals, treat them as a separate breed.
"""
function check_abnormal_elements()
    title("Check dv_2")
    fra = joinpath(work_dir, "data/genotypes/grm/all")
    til = joinpath(work_dir, "data/genotypes/grm/378")
    isdir(til) || mkdir(til)
    gt = joinpath(til, "genotypes.txt")
    isfile(gt) || symlink(joinpath(til, "genotypes.txt"), gt)

    idx = open(joinpath(til, "country.idx",), "w")
    open(joinpath(fra, "country.idx"), "r") do io
        for _ in 1:378
            line = readline(io)
            write(idx, line, " 0\n")
        end
        for line in eachline(io)
            id, b, c, d = split(line)
            write(idx, "$id 0 $b $c $d\n")
        end
    end
    close(idx)
    open(joinpath(til, "calc_grm.inp"), "w") do io
        write(io, join([
            "44037",
            "genotypes.txt",
            "genotypes",
            "4 country.idx",
            "vanraden",
            "giv 0.00",
            "G ASReml",
            "print_giv=asc",
            "print_geno=no genotypes.dat",
            "12"], '\n'))
    end

    cd(til)
    _ = run(`calc_grm`)
    _ = read(run(pipeline("G.grm", `gawk '{if($1==$2) print $3}'`, "diag.txt")), String)
    cd(work_dir)
end


text:
## variables (general)

function MOI.add_variable(o::Optimizer)
    allow_modification(o)
    vr = add_variable(o.inner)
    var::Ptr{SCIP_VAR} = o.inner.vars[vr][]
    register!(o, var, vr)
    return MOI.VariableIndex(vr.val)
end

MOI.add_variables(o::Optimizer, n) = [MOI.add_variable(o) for i=1:n]
MOI.get(o::Optimizer, ::MOI.NumberOfVariables) = length(o.inner.vars)
MOI.get(o::Optimizer, ::MOI.ListOfVariableIndices) = sort!([VI(k.val) for k in keys(o.inner.vars)], by=v->v.value)
MOI.is_valid(o::Optimizer, vi::VI) = haskey(o.inner.vars, VarRef(vi.value))

function MOI.get(o::Optimizer, ::MOI.VariableName, vi::VI)::String
    return GC.@preserve o unsafe_string(SCIPvarGetName(var(o, vi)))
end

function MOI.set(o::Optimizer, ::MOI.VariableName, vi::VI, name::String)
    @SCIP_CALL SCIPchgVarName(o, var(o, vi), name)
    return nothing
end

function MOI.delete(o::Optimizer, vi::VI)
    # Note: SCIP does not currently support deletion of variables that are still
    # present in some constraint. We don't want the overhead of keeping track of
    # the variable-in-constraint relation, so, to be conservative, we only allow
    # to delete a variable when there are no constraints in the model.
    if length(o.inner.conss) > 0
        error("Can not delete variable while model contains constraints!")
    end
    allow_modification(o)
    if !haskey(o.inner.vars, VarRef(vi.value))
        throw(MOI.InvalidIndex(vi))
    end
    # TODO still necessary?
    if !haskey(o.reference, var(o, vi))
        throw(MOI.InvalidIndex(vi))
    end
    haskey(o.binbounds, vi) && delete!(o.binbounds, vi)
    delete!(o.reference, var(o, vi))
    delete(o.inner, VarRef(vi.value))
    return nothing
end

## variable types (binary, integer)

MOI.supports_constraint(o::Optimizer, ::Type{VI}, ::Type{<:VAR_TYPES}) = true

scip_vartype(::Type{MOI.ZeroOne}) = SCIP_VARTYPE_BINARY
scip_vartype(::Type{MOI.Integer}) = SCIP_VARTYPE_INTEGER

function MOI.add_constraint(o::Optimizer, vi::VI, set::S) where {S <: VAR_TYPES}
    allow_modification(o)
    v = var(o, vi)
    infeasible = Ref{SCIP_Bool}()
    @SCIP_CALL SCIPchgVarType(o, v, scip_vartype(S), infeasible)
    # TODO: warn if infeasible[] == TRUE?
    if S <: MOI.ZeroOne
        # Need to adjust bounds for SCIP, which fails with an error otherwise.
        # Check for conflicts with existing bounds first:
        lb, ub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
        if lb == -SCIPinfinity(o) && ub == SCIPinfinity(o)
            @debug "Implicitly setting bounds [0,1] for binary variable at $(vi.value)!"
            @SCIP_CALL SCIPchgVarLb(o, v, 0.0)
            @SCIP_CALL SCIPchgVarUb(o, v, 1.0)
        else
            # Store old bounds for later recovery.
            o.binbounds[vi] = MOI.Interval(lb, ub)
            if ub > 1.0
                @debug "Tightening upper bound $ub to 1 for binary variable at $(vi.value)!"
                @SCIP_CALL SCIPchgVarUb(o, v, 1.0)
            end
            if lb < 0.0
                @debug "Tightening lower bound $lb to 0 for binary variable at $(vi.value)!"
                @SCIP_CALL SCIPchgVarLb(o, v, 0.0)
            end
        end
    end
    # use var index for cons index of this type
    i = vi.value
    return register!(o, CI{VI, S}(i))
end

function MOI.delete(o::Optimizer, ci::CI{VI,S}) where {S <: VAR_TYPES}
    _throw_if_invalid(o, ci)
    allow_modification(o)

    vi = VI(ci.value)
    v = var(o, vi)

    # Reset bounds from constraint if this was a binary, see below.
    reset_bounds = SCIPvarGetType(v) == SCIP_VARTYPE_BINARY

    # don't actually delete any SCIP constraint, just reset type
    infeasible = Ref{SCIP_Bool}()
    @SCIP_CALL SCIPchgVarType(o, var(o, vi), SCIP_VARTYPE_CONTINUOUS, infeasible)
    # TODO: warn if infeasible[] == TRUE?

    # Can only change bounds after changing the var type.
    if reset_bounds
        bounds = get(o.binbounds, vi, nothing)
        if bounds !== nothing
            @SCIP_CALL SCIPchgVarLb(o, v, bounds.lower)
            @SCIP_CALL SCIPchgVarUb(o, v, bounds.upper)
        end
    end

    # but do delete the constraint reference
    delete!(o.constypes[VI, S], ConsRef(ci.value))

    return nothing
end

## variable bounds

MOI.supports_constraint(o::Optimizer, ::Type{VI}, ::Type{<:BOUNDS}) = true

function MOI.add_constraint(o::Optimizer, vi::VI, set::S) where {S <: BOUNDS}
    allow_modification(o)
    v = var(o, vi)
    inf = SCIPinfinity(o)

    newlb, newub = bounds(set)
    newlb = newlb === nothing ? -inf : newlb
    newub = newub === nothing ?  inf : newub

    # Check for existing bounds first.
    oldlb, oldub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
    if (oldlb != -inf && newlb != -inf || oldub != inf && newub != inf)
        # special case: implicit bounds for binaries
        if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
            # Store new bounds
            o.binbounds[vi] = MOI.Interval(newlb, newub)

            if newlb < 0.0
                @debug "Ignoring relaxed lower bound $newlb for binary variable at $(vi.value)!"
                newlb = oldlb
            end

            if newub > 1.0
                @debug "Ignoring relaxed upper bounds $newub for binary variable at $(vi.value)!"
                newub = oldub
            end

            if newlb != oldlb || newub != oldub
                @debug "Overwriting existing bounds [0.0,1.0] with [$newlb,$newub] for binary variable at $(vi.value)!"
            end
        else # general case (non-binary variable)
            # TODO find initially-set bound constraint
            throw(MOI.LowerBoundAlreadySet{S,S}(vi))
        end
    end

    if newlb != -inf
        @SCIP_CALL SCIPchgVarLb(o, v, newlb)
    end

    if newub != inf
        @SCIP_CALL SCIPchgVarUb(o, v, newub)
    end

    # use var index for cons index of this type
    i = vi.value
    return register!(o, CI{VI, S}(i))
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{VI, S}) where S <: Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}
    @SCIP_CALL SCIPchgVarLb(o, v, lb)
    @SCIP_CALL SCIPchgVarUb(o, v, ub)
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{VI, MOI.GreaterThan{Float64}})
    @SCIP_CALL SCIPchgVarLb(o, v, lb)
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{VI, MOI.LessThan{Float64}})
    @SCIP_CALL SCIPchgVarUb(o, v, ub)
end

function MOI.delete(o::Optimizer, ci::CI{VI,S}) where {S <: BOUNDS}
    _throw_if_invalid(o, ci)
    allow_modification(o)

    # Don't actually delete any SCIP constraint, just reset bounds
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        reset_bounds(o, v, 0.0, 1.0, ci)
    else
        inf = SCIPinfinity(o)
        reset_bounds(o, v, -inf, inf, ci)
    end

    # but do delete the constraint reference
    haskey(o.binbounds, vi) && delete!(o.binbounds, vi)
    delete!(o.constypes[VI, S], ConsRef(ci.value))

    return nothing
end

function MOI.set(o::SCIP.Optimizer, ::MOI.ConstraintSet, ci::CI{VI,S}, set::S) where {S <: BOUNDS}
    allow_modification(o)
    v = var(o, VI(ci.value)) # cons index is actually var index
    lb, ub = bounds(set)
    lb = lb === nothing ? -SCIPinfinity(o) : lb
    ub = ub === nothing ?  SCIPinfinity(o) : ub
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        o.binbounds[vi] = MOI.Interval(lb, ub)
        lb = max(lb, 0.0)
        ub = min(lb, 1.0)
    end
    @SCIP_CALL SCIPchgVarLb(o, v, lb)
    @SCIP_CALL SCIPchgVarUb(o, v, ub)
    return nothing
end

# TODO: is actually wrong for unbounded variables?
function MOI.is_valid(o::Optimizer, ci::CI{VI,S}) where {S <: BOUNDS}
    if !haskey(o.inner.vars, VarRef(ci.value))
        return false
    end
    cons_set = get(o.constypes, (VI, S), nothing)
    if cons_set === nothing
        return false
    end
    return ConsRef(ci.value) in o.constypes[VI, S]
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VI, S}) where {S <: BOUNDS}
    if !MOI.is_valid(o, ci)
        throw(MOI.InvalidIndex(ci))
    end
    return VI(ci.value)
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VI, S}) where {S <: BOUNDS}
    if !MOI.is_valid(o, ci)
        throw(MOI.InvalidIndex(ci))
    end
    vi = VI(ci.value)
    v = var(o, vi)
    lb, ub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        bounds = o.binbounds[vi]
        lb, ub = bounds.lower, bounds.upper
    end
    return from_bounds(S, lb, ub)
end

function MOI.is_valid(o::Optimizer, ci::CI{VI,S}) where {S <: Union{MOI.ZeroOne, MOI.Integer}}
    if !haskey(o.inner.vars, VarRef(ci.value))
        return false
    end
    return ConsRef(ci.value) in o.constypes[VI, S]
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VI, MOI.ZeroOne})
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        return MOI.ZeroOne()
    end
    throw(MOI.InvalidIndex(ci))
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VI, MOI.Integer})
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_INTEGER
        return MOI.Integer()
    end
    throw(MOI.InvalidIndex(ci))
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VI, MOI.ZeroOne})
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        return vi
    end
    throw(MOI.InvalidIndex(ci))
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VI, MOI.Integer})
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_INTEGER
        return vi
    end
    throw(MOI.InvalidIndex(ci))
end

# (partial) warm starts

MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOI.get(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI)
    return if haskey(o.start, vi)
        o.start[vi]
    else
        nothing
    end
end

function MOI.set(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Float64)
    o.start[vi] = value
    return
end

function MOI.set(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Nothing)
    if haskey(o.start, vi)
        delete!(o.start, vi)
    end
    return
end

function MOI.set(::Optimizer, ::MOI.ConstraintFunction, ::CI{VI}, ::VI)
    throw(MOI.SettingVariableIndexNotAllowed())
end

function MOI.get(o::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    ptr = SCIPfindVar(o, name)
    if ptr == C_NULL
        return nothing
    end
    var_ref = get(o.reference, ptr, nothing)
    if var_ref === nothing
        return var_ref
    end
    return MOI.VariableIndex(var_ref.val)
end


text:
#!/usr/bin/env julia
using Test

include("test_cycle.jl")
include("test_demo3x7.jl")
include("test_lines.jl")
include("test_cxx.jl")
include("test_queues.jl")
include("test_dual.jl")
include("test_tree_dp.jl")
include("test_tree0.jl")
include("test_mgt.jl")
include("test_gaplas.jl")
include("test_gaplax.jl")
include("test_cyclgaplas.jl")
include("test_cyclegap.jl")
include("test_momentum.jl")

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text: Dict(zip(keys, values))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# 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.
using JSON

# =========================================================================================
"""
    mutable struct Languages
        version::Int64
        lastcheck::String
        server::String
        files::Dict{String, GitHubFile}
    end

This structure holds details about the language files used by Tesseract

__Values:__

| Name       | Description
| :--------- | :----------
| version    | The version number of the file on disk.
| lastcheck  | The time when we last checked with the server.
| server     | The URL of the server.
| files      | The details about all the training data files.

__Constructors:__

    Languages()

Create a new instance of the object.
"""
mutable struct Languages
    version::Int64
    lastcheck::String
    server::String
    files::Dict{String, GitHubFile}
    Languages() = new(INDEX_VERSION, "", "", Dict{String, GitHubFile}())
    Languages(version, lastcheck, server, files) = new(version, lastcheck, server, files)
end

# =========================================================================================
"""
    get_language_file(
        langs::Languages,
        language::AbstractString
    )::GitHubFile

Retrieve the file information for the specified file.  If the file is not found then a new
record is created an returned.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | langs    |         | The `Languages` structure to look for the file in.
| R | language |         | The language to get the file for.
"""
function get_language_file(
            langs::Languages,
            language::AbstractString
        )::GitHubFile
    return get!(langs.files, language) do
        GitHubFile(string(language, DATA_FILE_EXT), 0.0, "", "", "")
    end
end

# =========================================================================================
"""
has_all_languages(
        langs::Languages,
        languages::AbstractVector
    )::Bool

Check if we have all the specified languages defined.  If we don't have a definition for
one of the languages `false` is returned.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | langs     |         | The `Languages` object to look for the languages in.
| R | languages |         | The languages we're looking for.
"""
function has_all_languages(
            langs::Languages,
            languages::AbstractVector
        )::Bool
    for language in languages
        if haskey(langs.files, language) == false
            return false
        end
    end
    return true
end

# =========================================================================================
"""
    has_language(
        langs::Languages,
        language::AbstractString
    )::Bool

Check if we know about the specified language.  Returns `false` if the language file is is
not in the index.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | langs    |         | The `Languages` object to look for the language in.
| R | language |         | The language to check for.
"""
function has_language(
            langs::Languages,
            language::AbstractString
        )::Bool
    return haskey(langs.files, language)
end

# =========================================================================================
"""
    is_expired(
        langs::Languages,
        frequency::Integer
    )::Bool

Check if our languages index data is older than `frequency` days.  Returns `true` if the
data is older.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | langs     |         | The object to check.
| R | frequency |         | How old our data should be before we synchronize again.
"""
function is_expired(
            langs::Languages,
            frequency::Integer
        )::Bool
    if isempty(langs.lastcheck)
        return true
    elseif langs.server == ""
        return true
    elseif now() > DateTime(langs.lastcheck, DATE_FORMAT) + Day(frequency)
        return true
    end

    return false
end

# =========================================================================================
"""
    update_languages(
        lang::Languages
        source::GitHubProject,
    )::Nothing

Grab the latest URLs and SHA1 codes for the language data files from GitHub and update the
index file with them.

__Arguments:__

| T | Name   | Default | Description
|:--| :----- | :------ | :----------
| R | langs  |         | The language data to sync with the server.
| R | source |         | The GitHub project to sync the files with.

__Exceptions:__

| Name         | Reason
| :----------- | :-----
| NetworkError | There was a network error retrieving the contents from GitHub.
"""
function update_languages(
            langs::Languages,
            source::GitHubProject
        )::Nothing
    local data = gh_contents(source)

    if isa(data, AbstractArray) == false
        @warn "Unexpected response from the server: $(JSON.json(data))"
        throw(NetworkError(gh_url(source)))
    end

    local regex = r"^(.+)\.traineddata$"

    for file in data
        local name = get_json_safe(String, file, "name")
        local sha  = get_json_safe(String, file, "sha")
        local url  = get_json_safe(String, file, "download_url")

        if name === nothing
            @warn "Directory entry has no name!"
            continue
        elseif sha === nothing
            @warn "Directory entry has no sha!"
            continue
        end
        local m = match(regex, name)

        if m !== nothing && url !== nothing
            local language = m.captures[1]
            local details  = get_language_file(langs, language)

            details.remotesha1 = sha
            details.url        = url

            found = true
        elseif m !== nothing && url === nothing

            @warn "Data file has no download_url!"

        end
    end

    langs.server = gh_url(source)
    langs.lastcheck = Dates.format(now(), DATE_FORMAT)

    nothing
end

# =========================================================================================
"""
    refresh_languages(
        langs::Languages,
        source::GitHubProject,
        frequency::Integer,
        force::Bool
    )::Bool

Refresh our language data with the server if we are out of date.  Returns true if we
refreshed with the server.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | langs     |         | The `Languages` object to refresh.
| R | source    |         | The GitHub project/directory to synchronize with.
| R | frequency |         | How old our data should be before we synchronize again.
| R | force     |         | Flag that indicates we should refresh no mater what.

__Details:__

This method returns true under the following conditions:

  * We've never synchronized with the server.
  * The URL for the server is different from the URL we synchronized with last.
  * The last synchronization was more than `frequency` days ago.
"""
function refresh_languages(
            langs::Languages,
            source::GitHubProject,
            frequency::Integer,
            force::Bool
        )::Bool
    local result = false

    try
        # ---------------------------------------------------------------------------------
        # Update the training information if needed.
        if force || is_expired(langs, frequency)
            update_languages(langs, source)
            result = true
        end

    catch ex
        # ---------------------------------------------------------------------------------
        # If the exception is NOT a NetworkError rethrow it.
        if isa(ex, NetworkError) == false
            rethrow(ex)
        end
    end

    return result
end

# =========================================================================================
"""
    parse_languages(
        json::AbstractDict
    )::Languages

Create and populate a Languages object from a JSON object.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | json |         | The JSON object to extract the values from.

__Details:__

If JSON contains invalid information then an empty Languages structure is returned.
"""
function parse_languages(
            json::AbstractDict
        )::Languages

    # -------------------------------------------------------------------------------------
    # Grab the values we're interested in from the JSON.
    local version   = get_json_safe(Int64, json, "version")
    local lastcheck = get_json_safe(AbstractString, json, "lastcheck")
    local server    = get_json_safe(AbstractString, json, "server")
    local files     = get_json_safe(AbstractVector, json, "files")

    version === nothing && return Languages()
    lastcheck == nothing && return Languages()
    server == nothing && return Languages()
    files == nothing && return Languages()

    # -------------------------------------------------------------------------------------
    # Check the version.
    if version != INDEX_VERSION
        return Languages()
    end

    # -------------------------------------------------------------------------------------
    # Create the file map.
    local filemap = Dict{String, GitHubFile}()
    local regex   = r"^(.+)\.traineddata$"

    for f in files
        local file = parse_github_file(f)

        if file === nothing
            @warn "$LANGUAGES_FILE corrupt."
            return Languages()
        end
        local m = match(regex, file.filename)

        if m !== nothing
            filemap[m.captures[1]] = file
        end
    end

    return Languages(version, lastcheck, server, filemap)
end

# =========================================================================================
"""
    load_languages(
        path::AbstractString
    )::Languages

Load the languages data file from the specified directory.  If the file doesn't exist or is
corrupt a new empty object is returned.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | path |         | The path to load the languages file from.
"""
function load_languages(
            path::AbstractString
        )::Languages
    local result   = nothing
    local filename = joinpath(path, LANGUAGES_FILE)

    if create_data_dir(path) == true && isfile(filename) == true
        try
            # -----------------------------------------------------------------------------
            # Load and parse the file.
            local file = read(filename, String)
            local data = JSON.parse(file)

            if isa(data, AbstractDict) == true
                result = parse_languages(data)
            else
                # -------------------------------------------------------------------------
                # Bad JSON file.
                @warn "Could not parse: $filename"
                result = Languages()
            end
        catch ex
            # -----------------------------------------------------------------------------
            # Error reading the file.
            if isa(ex, SystemError) == true
                @warn "Could not read: $filename"
            elseif isa(ex, ErrorException) == false
                @warn "Could not parse: $filename"
            else
                rethrow(ex)
            end
            result = Languages()
        end
    else
        # ---------------------------------------------------------------------------------
        # File does not exist.
        result = Languages()
    end

    return result
end

# =========================================================================================
"""
    save_languages(
        langs::Languages,
        path::AbstractString
    )::Bool

Save languages file to the specified directory.  Returns `false` if the the operations
fails.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | langs |         | The languages data to save.
| R | path  |         | The path to save the file to.
"""
function save_languages(
            langs::Languages,
            path::AbstractString
        )::Bool
    local result   = false
    local filename = joinpath(path, LANGUAGES_FILE)
    local data     = Dict{String, Any}()

    data["version"]    = INDEX_VERSION
    data["lastcheck"]  = langs.lastcheck
    data["server"]     = langs.server
    data["files"]      = collect(values(langs.files))

    try
        write(filename, JSON.json(data))
        result = true
    catch ex
        if isa(ex, SystemError) == false
            rethrow(ex)
        end
    end

    return result
end


text:
function transform(input::AbstractDict)
    Dict(c=>score for (score, list) in input for c in lowercase.(list))
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism