text:
# DEBUG SEVERAL SOLUTIONS #####################################################

function DEBUG_correct_solutions(solutions::Vector{Sol})
    for sol in solutions
        DEBUG_correct_solution(sol)
    end
end

# DEBUG ONE MOMKP SOLUTION

function DEBUG_correct_solution(sol::Sol)

    momkp = sol.prob
    n = momkp.n

    @assert length(sol.x) == n "sol.x not of the right size"

    z_true = [0.0, 0.0]
    w_true = zeros(Float64, momkp.m)
    for var in 1:n
        @assert sol.x[var] == 0 || sol.x[var] == 1 "sol.x has values different from 0 or 1, so not a knapsack solution."
        if sol.x[var] == 1
            z_true[1] += momkp.P[1,var]
            z_true[2] += momkp.P[2,var]
            for dim in 1:momkp.m
                w_true[dim] += momkp.W[dim,var]
            end
        end
    end

    @assert z_true[1] == sol.z[1] && z_true[2] == sol.z[2] "sol.z != z_true, $(sol.z) != $z_true"
    @assert w_true == sol.w "sol.w != w_true, $(sol.w) != $w_true"
end

# DEBUG SEVERAL FEASIBLE SOLUTIONS #####################################################

function DEBUG_feasible_solutions(solutions::Vector{Sol})
    for sol in solutions
        DEBUG_feasible_solution(sol)
    end
end

# Debug functions are used to check assertions on the variables at runtime
# They are time consuming! Use only for debugging purposes!

# DEBUG ONE FEASIBLE MOMKP SOLUTION

function DEBUG_feasible_solution(sol::Sol)
    DEBUG_correct_solution(sol)

    for dim in 1:sol.prob.m
        @assert sol.w[dim] <= sol.prob.ω[dim] "Solution not feasible, sol.w[$dim] > sol.prob.ω[$dim] : $(sol.w[dim]) > $(sol.prob.ω[dim])"
    end
end

# VARIOUS DEBUG TESTS #####################################################

# check if the solutions are mutually non dominated
function DEBUG_is_front(front::Vector{Sol})
    n = length(front)

    for i in 1:n
        for j in i+1:n
            # if one dominates the other, error
            @assert !((front[i].z[1] >= front[j].z[1] && front[i].z[2] > front[j].z[2]) || (front[i].z[1] > front[j].z[1] && front[i].z[2] >= front[j].z[2])) "All solutions in front aren't mutually non-dominated"
            @assert !((front[j].z[1] >= front[i].z[1] && front[j].z[2] > front[i].z[2]) || (front[j].z[1] > front[i].z[1] && front[j].z[2] >= front[i].z[2])) "All solutions in front aren't mutually non-dominated"
        end
    end
end

function DEBUG_two_sols_sorted(sol1::Sol, sol2::Sol)
    @assert sol1.z[1] <= sol2.z[1] "Solutions not sorted wrt z1"
end

function DEBUG_correct_evaluations(population::Vector{Sol}; crowding = false, hv = false)

    pop_copy = copy(population)

    for sol in pop_copy # init values checked
        sol.rank = -1
        if crowding sol.crowding = -1 end
        if hv sol.hv_contrib = -1 end
    end

    F = INDICATOR_fast_non_dominated_sort!(pop_copy)
    for front in F
        crowding && INDICATOR_crowding_distance!(front)
        hv       && INDICATOR_hypervolume_contribution!(front)
    end

    for idx_sol in 1:length(population)
        @assert population[idx_sol].rank == pop_copy[idx_sol].rank "Population not well evaluated (wrong ranks)"
        crowding && @assert population[idx_sol].crowding == pop_copy[idx_sol].crowding "Population not well evaluated (wrong crowding)"
        hv       && @assert population[idx_sol].hv_contrib == pop_copy[idx_sol].hv_contrib "Population not well evaluated (wrong hypervolume contribution)"
    end
end

function DEBUG_unicity(population::Vector{Sol})

    for i in 1:length(population)
        for j in i+1:length(population)
            @assert population[i] != population[j] "No unicity in population"
            @assert population[i].z != population[j].z "No unicity : $(population[i]), $(population[j])"
        end
    end

end


text:
# These are keywords present in the TEXT section that are guaranteed by the spec
const required_keywords = [
    "\$BEGINANALYSIS", # Byte-offset to the beginning of the ANALYSIS segment.
    "\$BEGINDATA", # Byte-offset to the beginning of the DATA segment.
    "\$BEGINSTEXT", # Byte-offset to the beginning of a supplemental TEXT segment.
    "\$BYTEORD", # Byte order for data acquisition computer.
    "\$DATATYPE", # Type of data in DATA segment (ASCII, integer, floating point).
    "\$ENDANALYSIS", # Byte-offset to the last byte of the ANALYSIS segment.
    "\$ENDDATA", # Byte-offset to the last byte of the DATA segment.
    "\$ENDSTEXT", # Byte-offset to the last byte of a supplemental TEXT segment.
    "\$MODE", # Data mode (list mode - preferred, histogram - deprecated).
    "\$NEXTDATA", # Byte offset to next data set in the file.
    "\$PAR", # Number of parameters in an event.
    "\$PnB", # Number of bits reserved for parameter number n.
    "\$PnE", # Amplification type for parameter n.
    "\$PnN", # Short name for parameter n.
    "\$PnR", # Range for parameter number n.
    "\$TOT"  # Total number of events in the data set.
    ]

"""
    grab_word(iter, state, delimiter) -> word, state

Grabs the next word from the iterator `iter`. Takes care to handle escaped
delimiters specified by `delimiter`. Returns a string containing the word
and the state of the iterator.
"""
function grab_word(iter, state, delimiter::Char)
    word = Char[]
    prev = ' '
    iter_result = iterate(iter, state)
    while iter_result !== nothing
        i, state = iter_result

        # only add character if the current and previous are both
        # delimiters (i.e. escaped) or neither are
        if !xor((prev == delimiter), (i == delimiter))
            push!(word, i)
            prev = i
        else
            break
        end
        iter_result = iterate(iter, state)
    end
    join(word), state
end

"""
    verify_text(text_mappings) -> Void

Checks that all required keywords are present in the text_mapping dictionary
returned by `parse_text`
"""
function verify_text(text_mappings::Dict{String, String})
    # get all parameterized keywords $P1N, $P2N, etc
    is_param = [occursin("n", keyword) for keyword in required_keywords]

    # verify that all non-parameterized keywords are present in the mapping
    for non_param in required_keywords[.~is_param]
        if !haskey(text_mappings, non_param)
            error("FCS file is corrupted. It is missing required keyword $non_param in its TEXT section")
        end
    end

    # TODO: Add support for modes other than list
    (text_mappings["\$MODE"] != "L") && error("Non list mode FCS files are not supported yet")

    n_params = parse(Int, text_mappings["\$PAR"])

    for params in required_keywords[is_param]
        for i in 1:n_params
            if !haskey(text_mappings, replace(params, "n"=>i))
                error("FCS file is corrupted. It is missing required keyword $non_param in its TEXT section")
            end
        end
    end
end

function get_endian_func(text_mappings::Dict{String, String})
    byte_order = text_mappings["\$BYTEORD"]
    if byte_order == "1,2,3,4" # least significant byte first
        return ltoh
    elseif byte_order == "4,3,2,1" # most significant byte first
        return ntoh
    else
        error("FCS file is malformed. '$(byte_order)' is not a valid byte order.")
    end

end

ext: jl


text:
using SampleProjectTest
using SampleProject
using Test

@testset "Sum of series" begin
    @bench @testset "Geometric series" begin
        r = 0.9999
        s = geometric_series(1, r, 1_000_000)
        @test isapprox(s, 1 / (1 - r); atol = 1e-6)
    end

    @testset "Pi" begin
        @bench @testset "Gregory–Leibniz" begin
            s = pi_gregory_leibniz(1_000_000)
            @test isapprox(s, π; atol = 1e-6)
        end

        @bench @testset "Nilakantha" begin
            s = pi_nilakantha(100)
            @test isapprox(s, π; atol = 1e-6)
        end
    end
end

max_stars_repo_name: maxbennedich/julia-regression-analysis

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
using ProjAssistant
@quickactivate

## ----------------------------------------------------------------------------
@time begin
    using Serialization

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014
    const iJR = ChF.iJR904

    using Plots
    import GR
    !isinteractive() && GR.inline("png")

    using Serialization
    using Base.Threads

    using SimTools
    const SimT = SimTools

end

## ----------------------------------------------------------------------------
let
    monfile = cachedir(iJR, :ME_MONITOR)
    while !isfile(monfile); @warn("Mon file not found", monfile); sleep(10); end
    mon = SimT.OnDiskMonitor(monfile)
    live_proves = Dict()

    mean_min_val, mean_max_val = -1e2, 1e2

    fun(v) = clamp(v, mean_min_val, mean_max_val)

    SimT.watch(mon; wt = 15.0) do ddat

        isempty(ddat) && return

        for (exp, tdat) in ddat
            method = get(tdat, :method, "")

            # Check activity
            old_lprove = get!(live_proves, exp, -1.0)
            new_lprove = get(tdat, :live_prove, -1.0)
            old_lprove == new_lprove && continue
            @info("Doing", exp, old_lprove, new_lprove)
            live_proves[exp] = new_lprove

            for datk in [:round, :gd]
                kdat = get!(tdat, datk, Dict())
                ps = Plots.Plot[]

                # means
                for (avk, limk) in [
                            (:vg_avPME, :cgD_X),
                            (:biom_avPME, :exp_growth),
                        ]

                    avdat = get(kdat, avk, [])
                    p = plot(fun.(avdat);  title = string(avk),
                        xlabel = "iter", ylabel = string(avk),
                        lw = 3, label = string(avk)
                    )

                    limdat = get(tdat, limk, 0.0)
                    hline!(p, [fun(limdat)]; label = string(limk),
                        lw = 3, ls = :dash, color = :black,
                    )
                    plot!(p; legend = :topleft)
                    push!(ps, p)
                end

                # betas
                for bk in [:vg_beta, :biom_beta]
                    dat = get(kdat, bk, [])
                    p = plot(dat;  title = string(bk),
                        xlabel = "iter", ylabel = string(bk),
                        lw = 3, label = string(bk)
                    )
                    plot!(p; legend = :topleft)
                    push!(ps, p)
                end

                sfig(iJR, ps,
                    @fileid, "maxent_monitor",
                    (;datk, exp, method),
                    ".png"
                ); @info("Done", exp, datk)

            end # for datk in

        end # for (exp, tdat)

    end
end


text:
using Publish
using DataLoaders

p = Publish.Project(DataLoaders)

# needed to prevent error when overwriting
rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

# build documentation
deploy(DataLoaders; root = "/DataLoaders.jl", force = true, label = "dev")


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

using Markdown
using InteractiveUtils

# ╔═╡ 5fed24f3-0f30-49d9-857a-52546d0b5a7b
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.Registry.update()
	Pkg.add("Yao")
	Pkg.add("YaoPlots")
end

# ╔═╡ fe142e4e-46dd-11eb-3b8b-79a5e76789a5
using Yao, YaoPlots

# ╔═╡ fb8094ea-46da-11eb-3653-4b38fefb019c
md"# qRAM and Uncomputation"

# ╔═╡ 259b2164-46db-11eb-1969-19fbb8a1d18c
md"
*A random access memory (RAM) uses n bits to randomly address  N = `` 2^n `` distinct memory cells. A quantum random access memory (qRAM) uses n qubits to address any quantum superposition of N memory cells.*
Consider encoding data in qubits. Assume the following array."

# ╔═╡ 99efc04a-46dd-11eb-1084-23f01b24d6e0
a = [0, 0, 1, 0, 0]

# ╔═╡ cdab262c-46dd-11eb-1d74-d579e68851ee
md"It'll take 5 qubits to encode this."

# ╔═╡ dc1b6eec-46dd-11eb-324e-5dc8a639b138
ArrayReg(bit"00100") |> r->measure(r)

# ╔═╡ 0e7ff754-46de-11eb-13d5-35b80cd84409
md"or"

# ╔═╡ 160c0030-46de-11eb-0970-593c919558e4
zero_state(5) |> put(5, 3=>X) |> r->measure(r)

# ╔═╡ 4eeb3288-46de-11eb-18d5-59ff9c23a0b0
md"Either way, lets assume we've to encode 4 such arrays, in qubits."

# ╔═╡ 67c39c3a-46de-11eb-009f-9396ff599539
b = [0, 1, 1, 0, 1]

# ╔═╡ 868c2a06-46de-11eb-1b0b-f15ec79de048
c = [1, 1, 0, 0, 0]

# ╔═╡ be7b3a6a-46de-11eb-1ecd-4b300e1c4735
d = [1, 0, 1, 1, 1]

# ╔═╡ 95399610-46de-11eb-3876-253ee5fea025
md"To encode these 4 arrays, it'd take 20 qubits, judging by the above approach. But using QRAMS, we can use 7 qubits, to encode all the 4 arrays."

# ╔═╡ f7efa09c-46de-11eb-293a-e18468e3854d
begin
	f(x) = chain(7, [control(1:2, (k+2)=>X) for k in findall(isone, x)])
	QRAM = chain(7, repeat(H, 1:2), repeat(X, 1:2), f(a), repeat(X, 1:2), put(1=>X), f(b), put(1=>X), put(2=>X), f(c), put(2=>X), f(d))
	plot(QRAM)
end

# ╔═╡ 1811259c-46e1-11eb-23f2-31d638eb73b3
md"The first two qubits, are called the address qubits.
 - When the address qubits give `` 00 `` or `` 0 `` in decimal, for the next 5 qubits, we get $ $a $.
 - When the address qubits give `` 01 `` or `` 1 `` in decimal, for the next 5 qubits, we get $ $b $.
 - When the address qubits give `` 10 `` or `` 2 `` in decimal, for the next 5 qubits, we get $ $c $.
 - When the address qubits give `` 11 `` or `` 3 `` in decimal, for the next 5 qubits, we get $ $d $."

# ╔═╡ e26bbe6a-46e1-11eb-2fdc-e94b3c7884ee
md"The input to the QRAM is `` |0000000〉 ``."

# ╔═╡ 123ffdb8-46e2-11eb-2c4f-031a675cd551
output = zero_state(7) |> QRAM |> r->measure(r, nshots = 1024)

# ╔═╡ d892e034-46e2-11eb-1e8b-83b6b547961e
begin
	using StatsBase: fit, Histogram
	hist = fit(Histogram, Int.(output), 0:2^7)
	o1 = hist.weights[findall(!iszero, hist.weights)]
	o2 = reverse.(string.(0:(2^7-1), base=2, pad=7)[findall(!iszero, hist.weights)])
end

# ╔═╡ 8d630ef4-46e2-11eb-3476-f7b2bad5da0e
md"The below code records the frequency of measurements."

# ╔═╡ 59c077c4-46e4-11eb-1335-75104c6b0d26
md"
- When address qubits were $ $(o2[1][1:2]) $, we got $ $(o2[1][3:end]) $, with the frequency $ $(o1[1]) $.
- When address qubits were $ $(o2[2][1:2]) $, we got $ $(o2[2][3:end]) $, with the frequency $ $(o1[2]) $.
- When address qubits were $ $(o2[3][1:2]) $, we got $ $(o2[3][3:end]) $, with the frequency $ $(o1[3]) $.
- When address qubits were $ $(o2[4][1:2]) $, we got $ $(o2[4][3:end]) $, with the frequency $ $(o1[4]) $."

# ╔═╡ a431f16a-46e5-11eb-0344-3de152d77226
md"We use Uncomputation to reverse the everything we do in a circuit to reverse it to its former state. It's often used with qRAMs.

An example would be "

# ╔═╡ 83972a30-4756-11eb-3dfa-4d1222bcf547
plot(chain(7, put(1:7 => label(QRAM,"qRAM")), put(1:7 => label(Daggered(QRAM),"qRAM†"))))

# ╔═╡ d8b3b16e-4756-11eb-0e23-37f97d9795c8
begin
	uncomputation = chain(7, put(1:7 => QRAM), put(1:7 => QRAM'))
	plot(uncomputation)
end

# ╔═╡ 35c5548e-4757-11eb-24b0-1d2ea7aa3eb9
zero_state(7) |> uncomputation |> r->measure(r, nshots=1024)

# ╔═╡ 52c7a73a-4757-11eb-3112-e321f2878da4
md"As you can see, we first apply the QRAM circuit, and then its dagger, which undoes the effect."

# ╔═╡ Cell order:
# ╠═5fed24f3-0f30-49d9-857a-52546d0b5a7b
# ╟─fb8094ea-46da-11eb-3653-4b38fefb019c
# ╟─259b2164-46db-11eb-1969-19fbb8a1d18c
# ╠═99efc04a-46dd-11eb-1084-23f01b24d6e0
# ╟─cdab262c-46dd-11eb-1d74-d579e68851ee
# ╠═fe142e4e-46dd-11eb-3b8b-79a5e76789a5
# ╠═dc1b6eec-46dd-11eb-324e-5dc8a639b138
# ╟─0e7ff754-46de-11eb-13d5-35b80cd84409
# ╠═160c0030-46de-11eb-0970-593c919558e4
# ╟─4eeb3288-46de-11eb-18d5-59ff9c23a0b0
# ╠═67c39c3a-46de-11eb-009f-9396ff599539
# ╠═868c2a06-46de-11eb-1b0b-f15ec79de048
# ╠═be7b3a6a-46de-11eb-1ecd-4b300e1c4735
# ╟─95399610-46de-11eb-3876-253ee5fea025
# ╠═f7efa09c-46de-11eb-293a-e18468e3854d
# ╟─1811259c-46e1-11eb-23f2-31d638eb73b3
# ╟─e26bbe6a-46e1-11eb-2fdc-e94b3c7884ee
# ╠═123ffdb8-46e2-11eb-2c4f-031a675cd551
# ╟─8d630ef4-46e2-11eb-3476-f7b2bad5da0e
# ╠═d892e034-46e2-11eb-1e8b-83b6b547961e
# ╟─59c077c4-46e4-11eb-1335-75104c6b0d26
# ╟─a431f16a-46e5-11eb-0344-3de152d77226
# ╠═83972a30-4756-11eb-3dfa-4d1222bcf547
# ╠═d8b3b16e-4756-11eb-0e23-37f97d9795c8
# ╠═35c5548e-4757-11eb-24b0-1d2ea7aa3eb9
# ╟─52c7a73a-4757-11eb-3112-e321f2878da4


text:
module OdeIntegrators

using LinearAlgebra
using RigidBodyDynamics
using StaticArrays
using DocStringExtensions
using LoopThrottle

export runge_kutta_4,
    MuntheKaasIntegrator,
    ButcherTableau,
    OdeResultsSink,
    RingBufferStorage,
    ExpandingStorage,
    integrate,
    step

"""
$(TYPEDEF)

A [Butcher tableau](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Explicit_Runge.E2.80.93Kutta_methods).
"""
struct ButcherTableau{N, T, L}
    a::SMatrix{N, N, T, L}
    b::SVector{N, T}
    c::SVector{N, T}
    explicit::Bool

    function ButcherTableau(a::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
        N = size(a, 1)
        L = N^2
        @assert N > 0
        @assert size(a, 2) == N
        @assert length(b) == N
        c = vec(sum(a, dims=2))
        explicit = all(triu(a) .== 0)
        new{N, T, L}(SMatrix{N, N}(a), SVector{N}(b), SVector{N}(c), explicit)
    end
end
num_stages(::ButcherTableau{N}) where {N} = N
isexplicit(tableau::ButcherTableau) = tableau.explicit

"""
$(SIGNATURES)

Return the Butcher tableau for the standard fourth order Runge-Kutta integrator.
"""
function runge_kutta_4(scalar_type::Type{T}) where {T}
    a = zeros(T, 4, 4)
    a[2, 1] = 1/2
    a[3, 2] = 1/2
    a[4, 3] = 1
    b = T[1/6, 1/3, 1/3, 1/6]
    ButcherTableau(a, b)
end

"""
$(TYPEDEF)

Does 'something' with the results of an ODE integration (e.g. storing results,
visualizing, etc.). Subtypes must implement:
* `initialize(sink, state)`: called with the initial state when integration begins.
* `process(sink, t, state)`: called at every integration time step with the current state and time.
"""
abstract type OdeResultsSink end
initialize(::OdeResultsSink, state) = error("concrete subtypes must implement")
process(::OdeResultsSink, t, state) = error("concrete subtypes must implement")

"""
$(TYPEDEF)

An `OdeResultsSink` that stores the state at each integration time step in a
ring buffer.
"""
mutable struct RingBufferStorage{T, Q<:AbstractVector, V<:AbstractVector} <: OdeResultsSink
    ts::Vector{T}
    qs::Vector{Q}
    vs::Vector{V}
    last_index::Int64

    function RingBufferStorage{T}(state, n) where {T}
        Q = typeof(configuration(state))
        V = typeof(velocity(state))
        ts = Vector{T}(undef, n)
        qs = [similar(configuration(state)) for i in 1 : n]
        vs = [similar(velocity(state)) for i in 1 : n]
        new{T, Q, V}(ts, qs, vs, 0)
    end
end
Base.eltype(storage::RingBufferStorage{T}) where {T} = T
Base.length(storage::RingBufferStorage) = length(storage.ts)

function initialize(storage::RingBufferStorage, t, state)
    process(storage, t, state)
end

function process(storage::RingBufferStorage, t, state)
    index = storage.last_index % length(storage) + 1
    storage.ts[index] = t
    copyto!(storage.qs[index], configuration(state))
    copyto!(storage.vs[index], velocity(state))
    storage.last_index = index
    nothing
end

"""
$(TYPEDEF)

An `OdeResultsSink` that stores the state at each integration time step in
`Vectors` that may expand.
"""
mutable struct ExpandingStorage{T, Q<:AbstractVector, V<:AbstractVector} <: OdeResultsSink
    ts::Vector{T}
    qs::Vector{Q}
    vs::Vector{V}

    function ExpandingStorage{T}(state, n) where {T}
        Q = typeof(configuration(state))
        V = typeof(velocity(state))
        ts = T[]; sizehint!(ts, n)
        qs = Q[]; sizehint!(qs, n)
        vs = V[]; sizehint!(vs, n)
        new{T, Q, V}(ts, qs, vs)
    end
end
Base.eltype(storage::ExpandingStorage{T}) where {T} = T
Base.length(storage::ExpandingStorage) = length(storage.ts)

initialize(storage::ExpandingStorage, t, state) = process(storage, t, state)

function process(storage::ExpandingStorage, t, state)
    push!(storage.ts, t)
    q = similar(configuration(state))
    copyto!(q, configuration(state))
    push!(storage.qs, q)
    v = similar(velocity(state))
    copyto!(v, velocity(state))
    push!(storage.vs, v)
    nothing
end

mutable struct MuntheKaasStageCache{N, T, Q<:AbstractVector, V<:AbstractVector, S<:AbstractVector}
    q0::Q # global coordinates
    vs::SVector{N, V} # velocity vector for each stage
    vds::SVector{N, V} # time derivatives of vs
    ϕs::SVector{N, V} # local coordinates around q0 for each stage
    ϕds::SVector{N, V} # time derivatives of ϕs
    ss::SVector{N, S} # additional state for each stage
    sds::SVector{N, S} # time derivatives of ss
    ϕstep::V # local coordinates around q0 after complete step
    vstep::V # velocity after complete step
    sstep::S # additional state after complete step

    function MuntheKaasStageCache{N, T}(state) where {N, T}
        q0 = similar(configuration(state), T)
        vs = SVector{N}((similar(velocity(state), T) for i in 1 : N)...)
        vds = SVector{N}((similar(velocity(state), T) for i in 1 : N)...)
        ϕs = SVector{N}((similar(velocity(state), T) for i in 1 : N)...)
        ϕds = SVector{N}((similar(velocity(state), T) for i in 1 : N)...)
        ss = SVector{N}((similar(additional_state(state), T) for i in 1 : N)...)
        sds = SVector{N}((similar(additional_state(state), T) for i in 1 : N)...)
        ϕstep = similar(velocity(state), T)
        vstep = similar(velocity(state), T)
        sstep = similar(additional_state(state), T)
        new{N, T, typeof(q0), typeof(ϕstep), typeof(sstep)}(q0, vs, vds, ϕs, ϕds, ss, sds, ϕstep, vstep, sstep)
    end
end

"""
$(TYPEDEF)

A Lie-group-aware ODE integrator.

`MuntheKaasIntegrator` is used to properly integrate the dynamics of globally
parameterized rigid joints (Duindam, Port-Based Modeling and Control for
Efficient Bipedal Walking Robots, 2006, Definition 2.9). Global parameterizations of e.g. ``SO(3)``
are needed to avoid singularities, but this leads to the problem that the tangent
space no longer has the same dimension as the ambient space of the global parameterization.
A Munthe-Kaas integrator solves this problem by converting back and forth
between local and global coordinates at every integration time step.

The idea is to do the dynamics and compute the stages of the integration scheme
in terms of local coordinates centered around the global parameterization of
the configuration at the end of the previous time step (e.g. exponential coordinates),
combine the stages into a new set of local coordinates as usual for Runge-Kutta methods,
and then convert the local coordinates back to global coordinates.

From [Iserles et al., 'Lie-group methods' (2000)](https://hal.archives-ouvertes.fr/hal-01328729).

Another useful reference is [Park and Chung, 'Geometric Integration on Euclidean Group with Application to Articulated Multibody Systems' (2005)](http://www.ent.mrt.ac.lk/iml/paperbase/TRO%20Collection/TRO/2005/october/7.pdf).
"""
struct MuntheKaasIntegrator{N, T, F, S<:OdeResultsSink, X, L, M<:MuntheKaasStageCache{N, T}}
    dynamics!::F # dynamics!(vd, sd, t, state), sets vd (time derivative of v) and sd (time derivative of s) given time t and state
    tableau::ButcherTableau{N, T, L}
    sink::S
    state::X
    stages::M

    @doc """
    $(SIGNATURES)

    Create a `MuntheKaasIntegrator` given:

    * a callable `dynamics!(vd, t, state)` that updates the joint acceleration vector `vd` at time `t` and in state `state`;
    * a [`ButcherTableau`](@ref) `tableau`, specifying the integrator coefficients;
    * an [`OdeResultsSink`](@ref) `sink` which processes the results of the integration procedure at each time step.

    `state` must be of a type for which the following functions are defined:

    * `configuration(state)`, returns the configuration vector in global coordinates;
    * `velocity(state)`, returns the velocity vector;
    * `additional_state(state)`, returns the vector of additional states;
    * `set_velocity!(state, v)`, sets velocity vector to `v`;
    * `set_additional_state!(state, s)`, sets vector of additional states to `s`;
    * `global_coordinates!(state, q0, ϕ)`, sets global coordinates in state based on local coordinates `ϕ` centered around global coordinates `q0`;
    * `local_coordinates!(ϕ, ϕd, state, q0)`, converts state's global configuration `q` and velocity `v` to local coordinates centered around global coordinates `q0`.
    """ ->
    function MuntheKaasIntegrator(state::X, dynamics!::F, tableau::ButcherTableau{N, T, L}, sink::S) where {N, T, F, S<:OdeResultsSink, X, L}
        @assert isexplicit(tableau)
        stages = MuntheKaasStageCache{N, T}(state)
        new{N, T, F, S, X, L, typeof(stages)}(dynamics!, tableau, sink, state, stages)
    end
end

num_stages(::MuntheKaasIntegrator{N}) where {N} = N
Base.eltype(::MuntheKaasIntegrator{N, T}) where {N, T} = T

"""
$(SIGNATURES)

Take a single integration step.
"""
function step(integrator::MuntheKaasIntegrator, t::Real, Δt::Real)
    tableau = integrator.tableau
    stages = integrator.stages
    n = num_stages(integrator)

    # Use current configuration as the configuration around which the local coordinates for this step will be centered.
    q0, v0, s0 = stages.q0, stages.vstep, stages.sstep
    state = integrator.state
    q0 .= configuration(state)
    v0 .= velocity(state)
    s0 .= additional_state(state)

    # Compute integrator stages.
    for i = 1 : n
        # Update local coordinates and velocities
        ϕ = stages.ϕs[i]
        v = stages.vs[i]
        s = stages.ss[i]
        ϕ .= 0
        v .= v0
        s .= s0
        for j = 1 : i - 1
            aij = tableau.a[i, j]
            if aij != zero(aij)
                weight = Δt * aij
                ϕ .= muladd.(weight, stages.ϕds[j], ϕ)
                v .= muladd.(weight, stages.vds[j], v)
                s .= muladd.(weight, stages.sds[j], s)
            end
        end

        # Convert from local to global coordinates and set state
        # TODO: multiple setdirty! calls when using MechanismState:
        global_coordinates!(state, q0, ϕ)
        set_velocity!(state, v)
        set_additional_state!(state, s)

        # Dynamics in global coordinates
        vd = stages.vds[i]
        sd = stages.sds[i]
        integrator.dynamics!(vd, sd, muladd(tableau.c[i], Δt, t), state)

        # Convert back to local coordinates
        ϕd = stages.ϕds[i] # TODO: ϕ not actually needed, just ϕd!
        local_coordinates!(ϕ, ϕd, state, q0)
    end

    # Combine stages
    ϕ = stages.ϕstep
    ϕ .= 0
    v = stages.vstep # already initialized to v0
    s = stages.sstep # already initialized to s0
    for i = 1 : n
        weight = tableau.b[i] * Δt
        ϕ .= muladd.(weight, stages.ϕds[i], ϕ)
        v .= muladd.(weight, stages.vds[i], v)
        s .= muladd.(weight, stages.sds[i], s)
    end

    # Convert from local to global coordinates
    # TODO: multiple setdirty! calls when using MechanismState:
    global_coordinates!(state, q0, ϕ)
    set_velocity!(state, v)
    set_additional_state!(state, s)

    nothing
end

"""
$(SIGNATURES)

Integrate dynamics from the initial state at time ``0`` to `final_time`
using step size `Δt`.
"""
function integrate(integrator::MuntheKaasIntegrator, final_time, Δt; max_realtime_rate::Float64 = Inf)
    t = zero(eltype(integrator))
    state = integrator.state
    initialize(integrator.sink, t, state)
    @throttle t while t < final_time
        step(integrator, t, Δt)
        t += Δt
        process(integrator.sink, t, state)
    end max_rate=max_realtime_rate min_sleep_time=1/60.
end

end # module


text:
test_group("Operations on DataFrames that have column groupings")

x = DataFrame(quote
  a = [1,2,3]
  b = [4,5,6]
end)
y = DataFrame(quote
  c = [1,2,3]
  d = [4,5,6]
end)

# set_group(x, "group1", ["a", "b"])
# set_group(y, "group2", ["c", "d"])

z = deepcopy(x)
# @test is_group(z, "group1")

z = cbind(x, y)
# @test is_group(z, "group1")
# @test is_group(z, "group2")

v = DataFrame(quote
  a = [5,6,7]
  b = [8,9,10]
end)
z = rbind({v, x})
# @test is_group(z, "group1")

z = rbind(v,x)
# @test is_group(z, "group1")

# Deleting columns removes any mention from groupings
delete!(x, "a")
# @test colnames(x) == ["b"]
# @test get_groups(x)["group1"] == ["b"]

## del calls ref, which properly deals with groupings
y = without(y, "c")
@test colnames(y) == ["d"]
#@test get_groups(y)["group2"] == ["d"]
z1 = z[[1]]
@test colnames(z1) == ["a"]
#@test get_groups(z1)["group1"] == ["a"]

z2 = z[:,[1,1,2]]
@test colnames(z2) == ["a", "a_1", "b"]
#@test get_groups(z2)["group1"] == ["a_1", "b"]

test_group("DataFrame assignment")
df1 = DataFrame(quote
    a = 1:5
    b2 = letters[1:5]
    v2 = randn(5)
end)
df2 = DataFrame(quote
    a = reverse([1:5])
    b2 = reverse(letters)[1:5]
    v2 = randn(5)
end)
df1[1:2,:] = df2[4:5,:]
@test df1[1:2,:] == df2[4:5,:]

test_group("Empty DataFrame constructors")
df = DataFrame(10, 5)
@assert nrow(df) == 10
@assert ncol(df) == 5
@assert typeof(df[:, 1]) == DataVector{Float64}

df = DataFrame(Int64, 10, 3)
@assert nrow(df) == 10
@assert ncol(df) == 3
@assert typeof(df[:, 1]) == DataVector{Int64}
@assert typeof(df[:, 2]) == DataVector{Int64}
@assert typeof(df[:, 3]) == DataVector{Int64}

df = DataFrame({Int64, Float64, ASCIIString}, 100)
@assert nrow(df) == 100
@assert ncol(df) == 3
@assert typeof(df[:, 1]) == DataVector{Int64}
@assert typeof(df[:, 2]) == DataVector{Float64}
@assert typeof(df[:, 3]) == DataVector{ASCIIString}

df = DataFrame({Int64, Float64, ASCIIString}, ["A", "B", "C"], 100)
@assert nrow(df) == 100
@assert ncol(df) == 3
@assert typeof(df[:, 1]) == DataVector{Int64}
@assert typeof(df[:, 2]) == DataVector{Float64}
@assert typeof(df[:, 3]) == DataVector{ASCIIString}

df = DataFrame(zeros(10, 5))
@assert nrow(df) == 10
@assert ncol(df) == 5
@assert typeof(df[:, 1]) == DataVector{Float64}

df = DataFrame(ones(10, 5))
@assert nrow(df) == 10
@assert ncol(df) == 5
@assert typeof(df[:, 1]) == DataVector{Float64}

df = DataFrame(eye(10, 5))
@assert nrow(df) == 10
@assert ncol(df) == 5
@assert typeof(df[:, 1]) == DataVector{Float64}

test_group("Other DataFrame constructors")
df = DataFrame([{"a"=>1, "b"=>'c'}, {"a"=>3, "b"=>'d'}, {"a"=>5}])
@assert nrow(df) == 3
@assert ncol(df) == 2
@assert typeof(df[:,"a"]) == DataVector{Int}
@assert typeof(df[:,"b"]) == DataVector{Char}

df = DataFrame([{"a"=>1, "b"=>'c'}, {"a"=>3, "b"=>'d'}, {"a"=>5}], ["a", "b"])
@assert nrow(df) == 3
@assert ncol(df) == 2
@assert typeof(df[:,"a"]) == DataVector{Int}
@assert typeof(df[:,"b"]) == DataVector{Char}

data = {"A" => [1, 2], "C" => ["1", "2"], "B" => [3, 4]}
df = DataFrame(data)
# Specify column_names
df = DataFrame(data, ["C", "A", "B"])

# This assignment was missing before
df = DataFrame(quote Column = ["A"] end)
df[1, "Column"] = "Testing"

# flipud() tests
df = DataFrame(quote A = 1:4 end)
@assert isequal(flipud(df), DataFrame(quote A = [4, 3, 2, 1] end))
@assert isequal(df, DataFrame(quote A = [1, 2, 3, 4] end))
flipud!(df)
@assert isequal(df, DataFrame(quote A = [4, 3, 2, 1] end))


text:
"""
The PhilipsHue module provides functions for controlling the lights on a PhilipsHue bridge.

Quick start:
```
B = PhilipsHueBridge("192.168.1.2")
initialize(B; devicetype="juliascript#user1") # now run to the bridge!
testlights(B)
```
"""
module PhilipsHue

using JSON, HTTP, Colors

export  PhilipsHueBridge,
        getIP,
        getbridgeinfo,
        getbridgeconfig,
        isinitialized,
        getlights,
        getlight,
        getlightnumbers,
        setlight,
        setlights,
        randomcolors,
        testlights,
        register,
        initialize

mutable struct PhilipsHueBridge
    ip::AbstractString
    username:: AbstractString
    function PhilipsHueBridge(ip, username = "")
        fields = split(ip,'.')
        if length(fields) != 4
            throw(ArgumentError("IP address must have exactly four components."))
        end
        for f in fields
            fint = parse(Int, f)
            if fint < 0 || fint > 255
                throw(ArgumentError("IP address components must be between 0 and 255."))
            end
        end
        new(ip, username)
    end
end

"""
Initialize a bridge for the first time, supplying a devicetype (app name).
Registering this script with the bridge may require you to run to the bridge
and press the button.

The returned 'username' is stored in the bridge.username. This is needed for
future use, so you should remember it.

Returns true or false.

For example:

    B = PhilipsHueBridge("192.168.1.2")
    initialize(B; devicetype="juliascript#user1")
    testlights(B)

    # B.username is something like "2e4bdae26d734a73aeec4c21d4fd6"

then in a future Julia session you can do:

    B = PhilipsHueBridge("192.168.1.2", "2e4bdae26d734a73aeec4c21d4fd6")
    testlights(B)

"""
function initialize(bridge::PhilipsHueBridge; devicetype="juliascript#user1")
    println("initialize(): Trying to get the IP address of the Philips bridge.")
    ipaddress = "192.168.1.2"
    try
        ipaddress = getIP()
    catch e
        println("error was $e")
        return false
    end
    bridge.ip = ipaddress
    println("initialize(): Found bridge at $(bridge.ip).")
    println("initialize(): Trying to register $devicetype with the bridge at $(bridge.ip)...")
    username = register(bridge.ip, devicetype=devicetype)
    if ! isempty(username)
        println("initialize(): Registration successful")
        # save username in bridge
        bridge.username = username
        println("your 'username' is $username")
        return true
    else
        @warn("initialize(): Registration failed")
        return false
    end
end

"""
    isinitialized(bridge::PhilipsHueBridge)

Return true if the bridge has been initialized, and there is a connection to the portal.
"""
function isinitialized(bridge::PhilipsHueBridge)
    result = false
    try
        if getbridgeconfig(bridge)["portalconnection"] == "connected"
            result = true
        else
            result = false
        end
    catch e
        println("Bridge isn't initialized: error was $e")
        result = false
    end
    return result
end
"""
    getIP()

Read the bridge's IP settings from the [meethue.com]("https://discovery.meethue.com/") website.
"""
function getIP()
    response = HTTP.request("GET", "https://discovery.meethue.com/")
    # this url sometimes redirect, we should follow...
    if response.status == 302
        println("trying curl instead, in case of redirects")
        bridgeinfo = JSON.parse(read(`curl -sL https://discovery.meethue.com/`, String))
    else
        bridgeinfo = JSON.parse(String(response.body))
    end
    return bridgeinfo[1]["internalipaddress"]
end

"""
    B = PhilipsHueBridge("192.168.1.90", "username")
    getbridgeconfig(B)

Read the current bridge configuration. For example:
"""
getbridgeconfig(bridge::PhilipsHueBridge) = getbridgeinfo(bridge, "config")

"""
    getbridgeinfo(bridge[, category::String])

Get some information from the bridge. Category can be one of:

- "lights" resource which contains all the light resources
- "groups" resource which contains all the groups
- "config" resource which contains all the configuration items
- "schedules" which contains all the schedules
- "scenes" which contains all the scenes
- "sensors" which contains all the sensors
- "rules" which contains all the rules

Default is "config".
"""
function getbridgeinfo(bridge::PhilipsHueBridge, category::String="config")
    categories = ["lights", "groups", "config", "schedules", "scenes", "sensors", "rules"]
    if category in categories
        response = HTTP.request("GET", "http://$(bridge.ip)/api/$(bridge.username)/$category")
        return JSON.parse(String(response.body))
    else
        error("category must be one of $categories")
    end
end

"""
    getlightnumbers(bridge::PhilipsHueBridge)

Return the numbers of lights connected to the bridge.

Returns eg:

    [1, 3, 5, 6]
"""
function getlightnumbers(bridge::PhilipsHueBridge)
  return map(f -> parse(Int, f), collect(keys(getlights(bridge))))
end

"""
    getlights(bridge::PhilipsHueBridge)

Return the current settings of all lights connected to the bridge.
"""
function getlights(bridge::PhilipsHueBridge)
    response = HTTP.request("GET", "http://$(bridge.ip)/api/$(bridge.username)/lights")
    return JSON.parse(String(response.body))
end

"""
    getlight(bridge::PhilipsHueBridge, light=1)

Return the settings of the specified light. Note that if you have four lights,
they are not necessarily going to be numbered 1, 2, 3, 4.

Returns response.

"""
function getlight(bridge::PhilipsHueBridge, light)
    if ! in(light, getlightnumbers(bridge))
        return("no light with number $(light). Try using getlightnumbers().")
    end
    response = HTTP.request("GET", "http://$(bridge.ip)/api/$(bridge.username)/lights/$(string(light))")
    return JSON.parse(String(response.body))
end

"""
    setlight(B, 1, Dict("on" => true))
    setlight(B, 3, Dict("on" => false))
    setlight(B, 2, Dict("on" => true, "sat" => 123, "bri" => 243, "hue" => 123)

Set a light by passing a dictionary of settings.

eg Dict{Any,Any}("on" => true, "sat" => 123, "bri" => 123, "hue" => 123),
"hue" is from 0 to 65535, "sat" and "bri" are saturation and brightness from 1 to 254,
0 is red, yellow is 12750, green is 25500, blue is 46920, etc.

If keys are omitted, that aspect of the light won't be changed.

Keys are AbstractStrings, values can be numeric and will get converted to AbstractStrings
"""
function setlight(bridge::PhilipsHueBridge, light::Int, settings::Dict)
    if ! in(light, getlightnumbers(bridge))
        return("no light with number $(light). Try using getlightnumbers().")
    end
    state = AbstractString[]
    for (k, v) in settings
        push!(state, ("\"$k\": $(string(v))"))
    end
    state = "{" * join(state, ",") * "}"
    response = HTTP.request("PUT", "http://$(bridge.ip)/api/$(bridge.username)/lights/$(string(light))/state", body="$(state)")
    return JSON.parse(String(response.body))
end

"""
    setlight(bridge::PhilipsHueBridge, light::Int, col::ColorTypes.Colorant)
    setlight(B, 1, Colors.RGB(0.75, 0.25, 0.75))
    setlight(B, 1, colorant"Pink")

Set color of a light using Colors.jl style colors. (You might need `using Colors`.)
"""
function setlight(bridge::PhilipsHueBridge, light::Int, col::Color)
    if ! in(light, getlightnumbers(bridge))
        return("no light with number $(light). Try using getlightnumbers().")
    end
    c = convert(Colors.HSV, col)
    h, s, v = round(Int, (c.h / 360) * 65535), round(Int, c.s * 255), round(Int, c.v * 255)
    setlight(bridge, light, Dict("on" => true, "sat" => s, "bri" => v, "hue" => h))
end

"""
    setlights(B, Dict("on" => true))
    setlights(B, Dict("on" => false))
    setlights(B, Dict("on" => true, "sat" => 123, "bri" => 243, "hue" => 123)

Set all lights in a group by passing a dictionary of settings.

eg Dict{Any,Any}("on" => true, "sat" => 123, "bri" => 123, "hue" => 123),
"hue" is from 0 to 65535 (?), "sat" and "bri" are saturation and brightness from 1 to 254,
0 is red, yellow is 12750, green is 25500, blue is 46920, etc.

If keys are omitted, that aspect of the light won't be changed.

Keys are AbstractStrings, values can be numeric and will get converted to AbstractStrings
"""
function setlights(bridge::PhilipsHueBridge, settings::Dict)
    state = AbstractString[]
    for (k, v) in settings
        push!(state,("\"$k\": $(string(v))"))
    end
    state = "{" * join(state, ",") * "}"
    response = HTTP.request("PUT", "http://$(bridge.ip)/api/$(bridge.username)/groups/0/action", body="$(state)")
    return JSON.parse(String(response.body))
end

"""
    register(bridge_ip; devicetype="juliascript", blankusername="")

Register the devicetype and 'username' with the bridge.

Quoth Philips: If the username is not provided, a random key will be
generated and returned in the response. Important! The
username will soon be deprecated in the bridge. It is
strongly recommended not to use this and use the randomly
generated bridge username.

So we'll return the randomly generated key, or "" on failure.
"""
function register(bridge_ip; devicetype="juliascript", blankusername="")
    response     = HTTP.request("POST", "http://$(bridge_ip)/api/", body="{\"devicetype\":\"$(devicetype)#$(blankusername)\"}")
    responsedata = JSON.parse(String(response.body))
    # responsedata is probably:
    # 1-element Array{Any,1}:
    # ["error"=>["type"=>101,"description"=>"link button not pressed","address"=>"/"]]
    if responsedata[1][first(keys(responsedata[1]))]["description"] == "link button not pressed"
        println("register(): Quick, you have ten seconds to press the button on the bridge!")
        sleep(10)
        response = HTTP.post("http://$(bridge_ip)/api/", body="{\"devicetype\":\"$(devicetype)#$(blankusername)\"}")
        responsedata = JSON.parse(String(response.body))
        if first(keys(responsedata[1])) == "success"
            println("register(): Successfully registered $devicetype with the bridge at $bridge_ip")
            # returns username which is randomly generated key
            username = responsedata[1]["success"]["username"]
            println("register(): username is $username")
            return username
        else
            warn("register(): Failed to register $devicetype#$blankusername with the bridge at $bridge_ip")
            return ""
        end
    end
end

"""
    testlights(bridge::PhilipsHueBridge, total=5)

Test all lights. Since not all Hue lights do color, not testing color here.
"""
function testlights(bridge::PhilipsHueBridge, total=5)
    setlights(bridge, Dict("on" => true, "hue" => 10000, "sat" => 0, "bri" => 254))
    for i in 1:total
        setlights(bridge, Dict("on" => false))
        sleep(1)
        setlights(bridge, Dict("on" => true))
        sleep(1)
    end
    for light in enumerate(getlights(bridge))
        setlight(bridge, light[1], Dict("on"=>false))
        sleep(1)
        setlight(bridge, light[1], Dict("on"=> true))
        sleep(1)
    end
    setlights(bridge, Dict("hue" => 10000, "sat" => 0, "bri" => 254))
end

"""
  randomcolors(bridge::PhilipsHueBridge, delay = 1, repetitions=10; showstatus=false)

Do a set of repetitions of random color settings for all lights.
"""
function randomcolors(bridge::PhilipsHueBridge, delay = 1, repetitions=10; showstatus=false)
    setlights(bridge, Dict("on" => true, "hue" => 10000, "sat" => 0, "bri" => 254))
    for r in 1:repetitions
        for light in enumerate(getlights(bridge))
            response = getlight(bridge, light[1])
            if length(response) != 2 # don't do this with white only bulbs
                setlight(bridge, light[1], Dict("on" => true, "sat" => rand(128:254), "hue" => rand(0:65200)))
            end
            sleep(delay)
        end
    end
    response = setlights(bridge, Dict("on" => true, "hue" => 10000, "sat" => 0, "bri" => 254))
end

end


text:
######## THIS FILE IS THE SAME AS `haplotype_pair.jl`
######## except it solves the least squares objective on only `keep (default 100)` unique
######## haplotypes, where the top haplotypes are selected by minimizing dist(x, h).

function haplopair_thin_BLAS2!(
    X::AbstractMatrix,
    H::AbstractMatrix;
    allele_freq::Union{AbstractVector, Nothing} = nothing,
    keep::Int = 100,
    # preallocated vectors
    happair1:: AbstractVector{Int32}   = ones(Int32, size(X, 2)),      # length n
    happair2:: AbstractVector{Int32}   = ones(Int32, size(X, 2)),      # length n
    hapscore:: AbstractVector{Float32} = Vector{Float32}(undef, size(X, 2)), # length n
    maxindx :: AbstractVector{Int32}   = Vector{Int32}(undef, keep),         # length keep
    maxgrad :: AbstractVector{Float32} = Vector{Float32}(undef, keep),       # length keep
    Xi      :: AbstractVector{Float32} = Vector{Float32}(undef, size(H, 1)), # length p
    N       :: AbstractVector{Float32} = Vector{Float32}(undef, keep),       # length keep
    # preallocated matrices
    Hk    :: AbstractMatrix{Float32}=Matrix{Float32}(undef,size(H,1),keep),     # p × keep
    M     :: AbstractMatrix{Float32}=Matrix{Float32}(undef,keep,keep),          # keep × keep
    Xwork :: AbstractMatrix{Float32}=Matrix{Float32}(undef,size(X,1),size(X,2)),# p × n
    Hwork :: AbstractMatrix{Float32}=convert(Matrix{Float32},H),                # p × d
    R     :: AbstractMatrix{Float32}=Matrix{Float32}(undef,size(H,2),size(X,2)),# d × n
    )

    p, n = size(X)
    d    = size(H, 2)

    # do global search for small problems
    if keep ≥ d
        return haplopair!(Xw_aligned, Hw_aligned, happair1=happair1,
            happair2=happair2, hapscore=hapscore, Xwork=Xwork, Hwork=Hwork,
            M=M)
    end

    # create views for matrices
    t6 = @elapsed begin
        if size(Hk, 1) != size(H, 1)
            Hk = Matrix{Float32}(undef, size(H, 1), keep)
            Xi = Vector{Float32}(undef, size(H, 1))
            Xwork = Matrix{Float32}(undef, p, n)
            Hwork = Matrix{Float32}(undef, p, d)
        end
        Rwork = view(R, 1:d, :)
        Hwork_view = view(Hwork, :, 1:d)
        copyto!(Hwork_view, H)
    end

    # initialize missing
    t5 = @elapsed initXfloat!(Xwork, X)

    # compute distances between each column of H and each column of X
    t1 = @elapsed begin
        if !isnothing(allele_freq)
            # Rij = wk * || H[:, i] - X[:, j] ||²
            pairwise!(Rwork, WeightedSqEuclidean(allele_freq), Hwork_view,
                Xwork, dims=2)
        else
            # Rij = || H[:, i] - X[:, j] ||²
            pairwise!(Rwork, SqEuclidean(), Hwork_view, Xwork, dims=2)
        end
    end

    t2 = t3 = 0.0
    @inbounds for i in 1:n
        # find top matching haplotypes for sample i
        t1 += @elapsed findbotr!(Rwork, i, maxgrad, maxindx)

        # sync Hk, Xi, M, N
        t2 += @elapsed begin
            for k in 1:keep
                col = maxindx[k]
                for j in 1:p
                    Hk[j, k] = Hwork_view[j, col]
                end
            end
            for j in eachindex(Xi)
                Xi[j] = Xwork[j, i]
            end
            update_M!(M, Hk)
            update_N!(N, Xi, Hk)
        end

        # search routine
        t3 += @elapsed begin
            hapscore[i], h1, h2 = haplopair!(M, N)
            happair1[i], happair2[i] = maxindx[h1], maxindx[h2]
        end

        # supply constant term in objective
        t2 += @elapsed hapscore[i] += dot(Xi, Xi)
    end

    t4 = 0.0 # no haplotype rescreening

    return t1, t2, t3, t4, t5, t6
end

"""
    findbotr!(A, col, maxval, index)

Find the smallest `length(index)` elements of column `A[:, col]`. `val` is
filled with the found smallest `r` elements in sorted order and `idx` is filled
with their corresponding indices.
"""
@inline function findbotr!(
    A      :: AbstractMatrix{T},
    col    :: Integer,
    minval :: AbstractVector{T},
    index  :: AbstractVector{<:Integer}
    ) where T <: Real
    fill!(minval, typemax(T))
    @inbounds for row in 1:size(A, 1)
        a = A[row, col]
        k = searchsortedlast(minval, a)
        if k < length(index)
            popinsert_rev!(minval, k+1, a)
            popinsert_rev!(index, k+1, row)
        end
    end
    nothing
end

"""
    popinsert_rev!(v, k, vk)

Move elements in `v[k:end-1]` to `v[k+1:end]` and insert `vk` at position `k` of
vector `v`.
"""
@inline function popinsert_rev!(v::AbstractVector, k::Integer, vk)
    @inbounds for i in Iterators.reverse(k+1:length(v))
        v[i] = v[i-1]
    end
    v[k] = vk
    v
end

function update_M!(M::AbstractMatrix, H::AbstractMatrix)
    d = size(M, 1)
    mul!(M, Transpose(H), H)
    for j in 1:d, i in 1:(j - 1) # off-diagonal
        M[i, j] = 2M[i, j] + M[i, i] + M[j, j]
    end
    for j in 1:d # diagonal
        M[j, j] *= 4
    end
    return nothing
end

function update_N!(N::AbstractVector, x::AbstractVector, H::AbstractMatrix)
    mul!(N, Transpose(H), x)
    @simd for I in eachindex(N)
        N[I] *= 2
    end
    return nothing
end

function haplopair!(
    M::AbstractMatrix{T},
    N::AbstractVector{T},
    ) where T <: Real

    d = length(N)
    hapmin, hap1, hap2 = typemax(T), 1, 1

    @inbounds for k in 1:d
        Nk = N[k]
        for j in 1:k
            score = M[j, k] - N[j] - Nk
            if score < hapmin
                hapmin, hap1, hap2 = score, j, k
            end
        end
    end

    return hapmin, hap1, hap2
end

# """
# Same as `haplopair_thin_BLAS2` but internally it uses only BLAS 3 calls and does not recompute `M`
# for every sample. But this implies the search routine in `haplopair!` is not cache aware.
# """
# function haplopair_thin_BLAS3!(
#     X::AbstractMatrix,
#     H::AbstractMatrix;
#     alt_allele_freq::Union{Nothing, AbstractVector{Float32}} = nothing,
#     keep::Int = 100,
#     # N::ElasticArray{Float32} = ElasticArray{Float32}(undef, size(X, 2), size(H, 2)), # n × d
#     happair1::AbstractVector = ones(Int, size(X, 2)),     # length n
#     happair2::AbstractVector = ones(Int, size(X, 2)),     # length n
#     hapscore::AbstractVector = zeros(Float32, size(X, 2)) # length n
#     )

#     Xwork = zeros(Float32, size(X, 1), size(X, 2))
#     Hwork = convert(Matrix{Float32}, H)
#     initXfloat!(Xwork, X)

#     p, n     = size(X)
#     d        = size(H, 2)
#     M        = zeros(Float32, d, d)
#     N        = zeros(Float32, n, d)
#     happairs = ones(Int, n), ones(Int, n)
#     hapscore = zeros(Float32, n)

#     t1, t2, t3 = haplopair_thin2!(Xwork, Hwork, M, N, happairs, hapscore, alt_allele_freq, keep)
#     t4 = 0 # no haplotype rescreening

#     return happairs, hapscore, t1, t2, t3, t4
# end

# function haplopair_thin_BLAS3(
#     X::AbstractMatrix{Float32},
#     H::AbstractMatrix{Float32},
#     M::AbstractMatrix{Float32},
#     N::AbstractMatrix{Float32},
#     happairs::Tuple{AbstractVector, AbstractVector},
#     hapscore::AbstractVector,
#     alt_allele_freq::Union{Nothing, AbstractVector{Float32}},
#     keep::Int
#     )

#     p, n, d = size(X, 1), size(X, 2), size(H, 2)

#     # compute distances between each column of H and each column of X
#     t1 = @elapsed begin
#         if !isnothing(alt_allele_freq)
#             map!(x -> 1 / (2 * x * (1 - x)), alt_allele_freq, alt_allele_freq) # scale by 1 / 2p(1-p)
#             R = pairwise(WeightedEuclidean(alt_allele_freq), H, X, dims=2)
#         else
#             R = pairwise(Euclidean(), H, X, dims=2) # Rij = d(H[:, i], X[:, j])
#         end
#     end

#     t2 = @elapsed begin
#         # assemble M (upper triangular only)
#         mul!(M, Transpose(H), H)
#         for j in 1:d, i in 1:(j - 1) # off-diagonal
#             M[i, j] = 2M[i, j] + M[i, i] + M[j, j]
#         end
#         for j in 1:d # diagonal
#             M[j, j] *= 4
#         end

#         # assemble N
#         mul!(N, Transpose(X), H)
#         @simd for I in eachindex(N)
#             N[I] *= 2
#         end
#     end

#     # computational routine
#     t3 = @elapsed haplopair_thin2!(happairs[1], happairs[2], hapscore, M, N, R, keep)

#     # supplement the constant terms in objective
#     t3 += @elapsed begin
#         @inbounds for j in 1:n
#             @simd for i in 1:p
#                 hapscore[j] += abs2(X[i, j])
#             end
#         end
#     end

#     return t1, t2, t3
# end

# function haplopair_thin_BLAS3!(
#     happair1::AbstractVector{Int},
#     happair2::AbstractVector{Int},
#     hapmin::AbstractVector{T},
#     M::AbstractMatrix{T},
#     N::AbstractMatrix{T},
#     R::AbstractMatrix{T},
#     keep::Int
#     ) where T <: Real

#     n, d = size(N)
#     perm = zeros(Int, d)
#     fill!(hapmin, typemax(T))
#     keep > d && (keep = d)

#     @inbounds for i in 1:n
#         # find top matching haplotypes for sample i
#         partialsortperm!(perm, view(R, :, i), keep) # perm[1:keep] = hap indices that best matches xi

#         for iter1 in 1:keep
#             k = perm[iter1]
#             for iter2 in 1:keep
#                 j = perm[iter2]
#                 k > j && continue # since M upper triangular

#                 score = M[k, j] - N[i, k] - N[i, j]
#                 if score < hapmin[i]
#                     hapmin[i], happair1[i], happair2[i] = score, k, j
#                 end
#             end
#         end
#     end

#     return nothing
# end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
using Dates

export def_parametre_input

"""
Def_Parametre_Input(name_folder,name_file_graph)

    Function used to define a file with all of the input data used in the modeling.
    This function creates a txt file named Input_Data.txt

    Input:
    - name_folder::String (name of the folder where the file is created)
    - name_file_graph::String (name of the graph used in the modeling)
    Output:
    - NaN

"""
function def_parametre_input(name_folder,name_file_graph::String)

    open("Res/$(lpad(name_folder,1))/Input_Data.txt", "a+") do io
        # Parameters of the modeling
        writedlm(io, ["Name_graph_file" name_file_graph])
        writedlm(io, ["ParamSimu.nb_day" ParamSimu.NB_DAY])
        writedlm(io, ["ParamSimu.nb_cont" ParamSimu.NB_CONT])
        writedlm(io, ["ParamSimu.nb_MC" ParamSimu.NB_MC])
        # Parameters of the covid
        writedlm(io, ["ParametreCovid.mu_incub" ParametreCovid.MU_INCUB])
        writedlm(io, ["ParametreCovid.sigma_incub" ParametreCovid.SIGMA_INCUB])
        writedlm(io, ["ParametreCovid.mu_dt_cont" ParametreCovid.MU_DT_CONT])
        writedlm(io, ["ParametreCovid.sigma_dt_cont" ParametreCovid.SIGMA_DT_CONT])
        writedlm(io, ["ParametreCovid.mu_deb" ParametreCovid.MU_DEB])
        writedlm(io, ["ParametreCovid.sigma_deb" ParametreCovid.SIGMA_DEB])
        writedlm(io, ["ParametreCovid.mu_prob" ParametreCovid.MU_PROB])
        writedlm(io, ["ParametreCovid.sigma_prob" ParametreCovid.SIGMA_PROB])
        # Parameters of the population
        writedlm(io, ["ParamPop.nb_out" ParamPop.NB_OUT])
        writedlm(io, ["ParamPop.nb_in" ParamPop.NB_IN])
        writedlm(io, ["ParamPop.taux_etud" ParamPop.TAUX_ETUD])
        writedlm(io, ["ParamPop.nbperclass" ParamPop.NBPERCLASS])
        writedlm(io, ["ParamPop.nbperamphi" ParamPop.NBPERAMPHI])
        writedlm(io, ["ParamPop.lambda_tissu" ParamPop.LAMBDA_TISSU])
        writedlm(io, ["ParamPop.m2_employ" ParamPop.M2_EMPLOY])
        # Module DistanceC
        writedlm(io, ["DistanceC.minD_tissu" DistanceC.MIN_TISSU])
        writedlm(io, ["DistanceC.maxD_tissu" DistanceC.MAX_TISSU])
        writedlm(io, ["DistanceC.minD_htissu" DistanceC.MIN_HTISSU])
        writedlm(io, ["DistanceC.maxD_htissu" DistanceC.MAX_HTISSU])
        writedlm(io, ["DistanceC.minD_menage" DistanceC.MIN_MENAGE])
        writedlm(io, ["DistanceC.maxD_menage" DistanceC.MAX_MENAGE])
        writedlm(io, ["DistanceC.minD_resto" DistanceC.MIN_RESTO])
        writedlm(io, ["DistanceC.maxD_resto" DistanceC.MAX_RESTO])
        # Module Journee
        #writedlm(io, ["Journee.Semaine.min_dep",Journee.Semaine.min_dep])
        #writedlm(io, ["Journee.Semaine.max_dep",Journee.Semaine.max_dep])
        #writedlm(io, ["Journee.Semaine.min_arr",Journee.Semaine.min_arr])
        #writedlm(io, ["Journee.Semaine.max_arr",Journee.Semaine.min_arr])
        #writedlm(io, ["Journee.Semaine.Type",Journee.Semaine.Type])
        #writedlm(io, ["Journee.Semaine.Deb_Type",Journee.Semaine.Deb_Type])
        #writedlm(io, ["Journee.Semaine.Fin_Type",Journee.Semaine.Fin_Type])
        #writedlm(io, ["Journee.Proba.sortieJournee",Journee.Proba.sortieJournee])
        #writedlm(io, ["Journee.Proba.resto",Journee.Proba.resto])
        #writedlm(io, ["Journee.Proba.sortieSoir",Journee.Proba.sortieSoir])
        #writedlm(io, ["Journee.Proba.sortieCourses",Journee.Proba.sortieCourses])

    end

end

export parallelisation
"""
Parallelisation()

    Function used to define the parallelization. The function returns the rank of the processus

    Input:
    - NaN
    Output:
    - comm
    - rank::Int

"""
function parallelisation()

    if(ParamSimu.vMPI)
        MPI.Init() ;
        comm = MPI.COMM_WORLD ;
        # Nombre de processus alloué à la parallélisation
        # nbproc = MPI.Comm_size(comm) ;
        # Numéro du rang du processus en cours
        rank = MPI.Comm_rank(comm) ;
        #println("Je suis le processus ", rank, " parmi ", nbproc , " processus")
    else
        rank = 0 ;
        println("La simulation ne tourne que sur un processus !")
        comm = 0 ;
    end

    return(comm,rank)

end

"""
Function DistanceResto
Fonction permettant de définir la distance la plus courte entre une entreprise et l'ensemble des services de restaurations.
"""
function distance_resto(id_aliment::Array{Int32,1},id_work::Array{Int32,1})

    simple_graph = loadgraph("Data/SaintVicNoChildrenLeftBehind.graphml", GraphIO.GraphML.GraphMLFormat()) ;

    chemin = zeros(Int32,(length(id_work),length(id_aliment))) ;

    @time for i = 1:length(id_work)
        for j=1:length(id_aliment)
            ch = a_star(simple_graph,id_work[i],id_aliment[j]) ;
            chemin[i,j] = length(ch) ;
        end

        println("J\'ai fini la",i,"eme entreprise")

    end

    f = writedlm("Data/Distance_resto.txt", chemin) ;

    return()

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const COMPILER_TEMP_SYM = Symbol("#temp#")

# add the real backedges
function finalize_backedges(frame::InferenceState)
    toplevel = !isa(frame.linfo.def, Method)
    if !toplevel && frame.cached && frame.max_valid == typemax(UInt)
        caller = frame.linfo
        for edges in frame.stmt_edges
            i = 1
            while i <= length(edges)
                to = edges[i]
                if isa(to, MethodInstance)
                    ccall(:jl_method_instance_add_backedge, Cvoid, (Any, Any), to, caller)
                    i += 1
                else
                    typeassert(to, Core.MethodTable)
                    typ = edges[i + 1]
                    ccall(:jl_method_table_add_backedge, Cvoid, (Any, Any, Any), to, typ, caller)
                    i += 2
                end
            end
        end
    end
end

# at the end, all items in b's cycle
# will now be added to a's cycle
function union_caller_cycle!(a::InferenceState, b::InferenceState)
    callers_in_cycle = b.callers_in_cycle
    b.parent = a.parent
    b.callers_in_cycle = a.callers_in_cycle
    contains_is(a.callers_in_cycle, b) || push!(a.callers_in_cycle, b)
    if callers_in_cycle !== a.callers_in_cycle
        for caller in callers_in_cycle
            if caller !== b
                caller.parent = a.parent
                caller.callers_in_cycle = a.callers_in_cycle
                push!(a.callers_in_cycle, caller)
            end
        end
    end
    return
end

function merge_call_chain!(parent::InferenceState, ancestor::InferenceState, child::InferenceState)
    # add backedge of parent <- child
    # then add all backedges of parent <- parent.parent
    # and merge all of the callers into ancestor.callers_in_cycle
    # and ensure that walking the parent list will get the same result (DAG) from everywhere
    while true
        add_backedge!(child, parent, parent.currpc)
        union_caller_cycle!(ancestor, child)
        child = parent
        parent = child.parent
        child === ancestor && break
    end
end

# Walk through `linfo`'s upstream call chain, starting at `parent`. If a parent
# frame matching `linfo` is encountered, then there is a cycle in the call graph
# (i.e. `linfo` is a descendant callee of itself). Upon encountering this cycle,
# we "resolve" it by merging the call chain, which entails unioning each intermediary
# frame's `callers_in_cycle` field and adding the appropriate backedges. Finally,
# we return `linfo`'s pre-existing frame. If no cycles are found, `nothing` is
# returned instead.
function resolve_call_cycle!(linfo::MethodInstance, parent::InferenceState)
    frame = parent
    uncached = false
    while isa(frame, InferenceState)
        uncached |= !frame.cached # ensure we never add an uncached frame to a cycle
        if frame.linfo === linfo
            uncached && return true
            merge_call_chain!(parent, frame, frame)
            return frame
        end
        for caller in frame.callers_in_cycle
            if caller.linfo === linfo
                uncached && return true
                merge_call_chain!(parent, frame, caller)
                return caller
            end
        end
        frame = frame.parent
    end
    return false
end

# build (and start inferring) the inference frame for the linfo
function typeinf_frame(linfo::MethodInstance,
                       optimize::Bool, cached::Bool, params::Params)
    frame = InferenceState(linfo, optimize, cached, params)
    frame === nothing && return nothing
    cached && (linfo.inInference = true)
    typeinf(frame)
    return frame
end

# compute (and cache) an inferred AST and return the current best estimate of the result type
function typeinf_edge(method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    code = code_for_method(method, atypes, sparams, caller.params.world)
    code === nothing && return Any, nothing
    code = code::MethodInstance
    if isdefined(code, :inferred)
        # return rettype if the code is already inferred
        # staged functions make this hard since they have two "inferred" conditions,
        # so need to check whether the code itself is also inferred
        inf = code.inferred
        if !isa(inf, CodeInfo) || (inf::CodeInfo).inferred
            if isdefined(code, :inferred_const)
                return AbstractEvalConstant(code.inferred_const), code
            else
                return code.rettype, code
            end
        end
    end
    if !caller.cached && caller.parent === nothing
        # this caller exists to return to the user
        # (if we asked resolve_call_cyle, it might instead detect that there is a cycle that it can't merge)
        frame = false
    else
        frame = resolve_call_cycle!(code, caller)
    end
    if frame === false
        # completely new
        code.inInference = true
        frame = InferenceState(code, #=optimize=#true, #=cached=#true, caller.params) # always optimize and cache edge targets
        if frame === nothing
            # can't get the source for this, so we know nothing
            code.inInference = false
            return Any, nothing
        end
        if caller.cached # don't involve uncached functions in cycle resolution
            frame.parent = caller
        end
        typeinf(frame)
        return frame.bestguess, frame.inferred ? frame.linfo : nothing
    elseif frame === true
        # unresolvable cycle
        return Any, nothing
    end
    frame = frame::InferenceState
    return frame.bestguess, nothing
end

#### entry points for inferring a MethodInstance given a type signature ####

# compute an inferred AST and return type
function typeinf_code(method::Method, @nospecialize(atypes), sparams::SimpleVector,
                      optimize::Bool, cached::Bool, params::Params)
    code = code_for_method(method, atypes, sparams, params.world)
    code === nothing && return (nothing, nothing, Any)
    return typeinf_code(code::MethodInstance, optimize, cached, params)
end
function typeinf_code(linfo::MethodInstance, optimize::Bool, cached::Bool,
                      params::Params)
    for i = 1:2 # test-and-lock-and-test
        i == 2 && ccall(:jl_typeinf_begin, Cvoid, ())
        if cached && isdefined(linfo, :inferred)
            # see if this code already exists in the cache
            # staged functions make this hard since they have two "inferred" conditions,
            # so need to check whether the code itself is also inferred
            if min_world(linfo) <= params.world <= max_world(linfo)
                inf = linfo.inferred
                if linfo.jlcall_api == 2
                    method = linfo.def::Method
                    tree = ccall(:jl_new_code_info_uninit, Ref{CodeInfo}, ())
                    tree.code = Any[ Expr(:return, quoted(linfo.inferred_const)) ]
                    tree.signature_for_inference_heuristics = nothing
                    tree.slotnames = Any[ COMPILER_TEMP_SYM for i = 1:method.nargs ]
                    tree.slotflags = UInt8[ 0 for i = 1:method.nargs ]
                    tree.slottypes = nothing
                    tree.ssavaluetypes = 0
                    tree.inferred = true
                    tree.pure = true
                    tree.inlineable = true
                    i == 2 && ccall(:jl_typeinf_end, Cvoid, ())
                    return svec(linfo, tree, linfo.rettype)
                elseif isa(inf, CodeInfo)
                    if inf.inferred
                        i == 2 && ccall(:jl_typeinf_end, Cvoid, ())
                        return svec(linfo, inf, linfo.rettype)
                    end
                end
            end
        end
    end
    frame = typeinf_frame(linfo, optimize, cached, params)
    ccall(:jl_typeinf_end, Cvoid, ())
    frame === nothing && return svec(nothing, nothing, Any)
    frame = frame::InferenceState
    frame.inferred || return svec(nothing, nothing, Any)
    frame.cached || return svec(nothing, frame.src, widenconst(frame.bestguess))
    return svec(frame.linfo, frame.src, widenconst(frame.bestguess))
end

# compute (and cache) an inferred AST and return the inferred return type
function typeinf_type(method::Method, @nospecialize(atypes), sparams::SimpleVector,
                      cached::Bool, params::Params)
    if contains_is(unwrap_unionall(atypes).parameters, Union{})
        return Union{}
    end
    code = code_for_method(method, atypes, sparams, params.world)
    code === nothing && return nothing
    code = code::MethodInstance
    for i = 1:2 # test-and-lock-and-test
        i == 2 && ccall(:jl_typeinf_begin, Cvoid, ())
        if cached && isdefined(code, :inferred)
            # see if this rettype already exists in the cache
            # staged functions make this hard since they have two "inferred" conditions,
            # so need to check whether the code itself is also inferred
            inf = code.inferred
            if !isa(inf, CodeInfo) || (inf::CodeInfo).inferred
                i == 2 && ccall(:jl_typeinf_end, Cvoid, ())
                return code.rettype
            end
        end
    end
    frame = typeinf_frame(code, cached, cached, params)
    ccall(:jl_typeinf_end, Cvoid, ())
    frame === nothing && return nothing
    frame = frame::InferenceState
    frame.inferred || return nothing
    return widenconst(frame.bestguess)
end

function typeinf_ext(linfo::MethodInstance, world::UInt)
    if isa(linfo.def, Method)
        # method lambda - infer this specialization via the method cache
        return typeinf_code(linfo, true, true, Params(world))
    else
        # toplevel lambda - infer directly
        ccall(:jl_typeinf_begin, Cvoid, ())
        result = InferenceResult(linfo)
        frame = InferenceState(result, linfo.inferred::CodeInfo,
                               true, true, Params(world))
        typeinf(frame)
        ccall(:jl_typeinf_end, Cvoid, ())
        @assert frame.inferred # TODO: deal with this better
        @assert frame.linfo === linfo
        linfo.rettype = widenconst(frame.bestguess)
        return svec(linfo, frame.src, linfo.rettype)
    end
end

#### do the work of inference ####

function typeinf_work(frame::InferenceState)
    @assert !frame.inferred
    frame.dont_work_on_me = true # mark that this function is currently on the stack
    W = frame.ip
    s = frame.stmt_types
    n = frame.nstmts
    while frame.pc´´ <= n
        # make progress on the active ip set
        local pc::Int = frame.pc´´ # current program-counter
        while true # inner loop optimizes the common case where it can run straight from pc to pc + 1
            #print(pc,": ",s[pc],"\n")
            local pc´::Int = pc + 1 # next program-counter (after executing instruction)
            if pc == frame.pc´´
                # need to update pc´´ to point at the new lowest instruction in W
                min_pc = next(W, pc)[2]
                if done(W, min_pc)
                    frame.pc´´ = max(min_pc, n + 1)
                else
                    frame.pc´´ = min_pc
                end
            end
            delete!(W, pc)
            frame.currpc = pc
            frame.cur_hand = frame.handler_at[pc]
            frame.stmt_edges[pc] === () || empty!(frame.stmt_edges[pc])
            stmt = frame.src.code[pc]
            changes = abstract_interpret(stmt, s[pc]::VarTable, frame)
            if changes === ()
                break # this line threw an error and so there is no need to continue
                # changes = s[pc]
            end
            if frame.cur_hand !== () && isa(changes, StateUpdate)
                # propagate new type info to exception handler
                # the handling for Expr(:enter) propagates all changes from before the try/catch
                # so this only needs to propagate any changes
                l = frame.cur_hand[1]
                if stupdate1!(s[l]::VarTable, changes::StateUpdate) !== false
                    if l < frame.pc´´
                        frame.pc´´ = l
                    end
                    push!(W, l)
                end
            end
            if isa(changes, StateUpdate)
                changes_var = changes.var
                if isa(changes_var, SSAValue)
                    # directly forward changes to an SSAValue to the applicable line
                    record_ssa_assign(changes_var.id + 1, changes.vtype.typ, frame)
                end
            elseif isa(stmt, NewvarNode)
                sn = slot_id(stmt.slot)
                changes = changes::VarTable
                changes[sn] = VarState(Bottom, true)
            elseif isa(stmt, GotoNode)
                pc´ = (stmt::GotoNode).label
            elseif isa(stmt, Expr)
                stmt = stmt::Expr
                hd = stmt.head
                if hd === :gotoifnot
                    condt = abstract_eval(stmt.args[1], s[pc], frame)
                    condval = maybe_extract_const_bool(condt)
                    l = stmt.args[2]::Int
                    changes = changes::VarTable
                    # constant conditions
                    if condval === true
                    elseif condval === false
                        pc´ = l
                    else
                        # general case
                        frame.handler_at[l] = frame.cur_hand
                        if isa(condt, Conditional)
                            changes_else = StateUpdate(condt.var, VarState(condt.elsetype, false), changes)
                            changes = StateUpdate(condt.var, VarState(condt.vtype, false), changes)
                        else
                            changes_else = changes
                        end
                        newstate_else = stupdate!(s[l], changes_else)
                        if newstate_else !== false
                            # add else branch to active IP list
                            if l < frame.pc´´
                                frame.pc´´ = l
                            end
                            push!(W, l)
                            s[l] = newstate_else
                        end
                    end
                elseif hd === :return
                    pc´ = n + 1
                    rt = abstract_eval(stmt.args[1], s[pc], frame)
                    if !isa(rt, Const) && !isa(rt, Type)
                        # only propagate information we know we can store
                        # and is valid inter-procedurally
                        rt = widenconst(rt)
                    end
                    if tchanged(rt, frame.bestguess)
                        # new (wider) return type for frame
                        frame.bestguess = tmerge(frame.bestguess, rt)
                        for (caller, caller_pc) in frame.backedges
                            # notify backedges of updated type information
                            if caller.stmt_types[caller_pc] !== ()
                                if caller_pc < caller.pc´´
                                    caller.pc´´ = caller_pc
                                end
                                push!(caller.ip, caller_pc)
                            end
                        end
                    end
                elseif hd === :enter
                    l = stmt.args[1]::Int
                    frame.cur_hand = (l, frame.cur_hand)
                    # propagate type info to exception handler
                    l = frame.cur_hand[1]
                    old = s[l]
                    new = s[pc]::Array{Any,1}
                    newstate_catch = stupdate!(old, new)
                    if newstate_catch !== false
                        if l < frame.pc´´
                            frame.pc´´ = l
                        end
                        push!(W, l)
                        s[l] = newstate_catch
                    end
                    typeassert(s[l], VarTable)
                    frame.handler_at[l] = frame.cur_hand
                elseif hd === :leave
                    for i = 1:((stmt.args[1])::Int)
                        frame.cur_hand = frame.cur_hand[2]
                    end
                end
            end
            pc´ > n && break # can't proceed with the fast-path fall-through
            frame.handler_at[pc´] = frame.cur_hand
            newstate = stupdate!(s[pc´], changes)
            if isa(stmt, GotoNode) && frame.pc´´ < pc´
                # if we are processing a goto node anyways,
                # (such as a terminator for a loop, if-else, or try block),
                # consider whether we should jump to an older backedge first,
                # to try to traverse the statements in approximate dominator order
                if newstate !== false
                    s[pc´] = newstate
                end
                push!(W, pc´)
                pc = frame.pc´´
            elseif newstate !== false
                s[pc´] = newstate
                pc = pc´
            elseif pc´ in W
                pc = pc´
            else
                break
            end
        end
    end
    frame.dont_work_on_me = false
end

function typeinf(frame::InferenceState)
    typeinf_work(frame)

    # If the current frame is part of a cycle, solve the cycle before finishing
    no_active_ips_in_callers = false
    while !no_active_ips_in_callers
        no_active_ips_in_callers = true
        for caller in frame.callers_in_cycle
            caller.dont_work_on_me && return
            if caller.pc´´ <= caller.nstmts # equivalent to `isempty(caller.ip)`
                # Note that `typeinf_work(caller)` can potentially modify the other frames
                # `frame.callers_in_cycle`, which is why making incremental progress requires the
                # outer while loop.
                typeinf_work(caller)
                no_active_ips_in_callers = false
            end
            if caller.min_valid < frame.min_valid
                caller.min_valid = frame.min_valid
            end
            if caller.max_valid > frame.max_valid
                caller.max_valid = frame.max_valid
            end
        end
    end

    # with no active ip's, type inference on frame is done

    if isempty(frame.callers_in_cycle)
        @assert !(frame.dont_work_on_me)
        frame.dont_work_on_me = true
        optimize(frame)
        finish(frame)
        finalize_backedges(frame)
    else # frame is in frame.callers_in_cycle
        for caller in frame.callers_in_cycle
            @assert !(caller.dont_work_on_me)
            caller.dont_work_on_me = true
        end
        # complete the computation of the src optimizations
        for caller in frame.callers_in_cycle
            optimize(caller)
            if frame.min_valid < caller.min_valid
                frame.min_valid = caller.min_valid
            end
            if frame.max_valid > caller.max_valid
                frame.max_valid = caller.max_valid
            end
        end
        # update and store in the global cache
        for caller in frame.callers_in_cycle
            caller.min_valid = frame.min_valid
        end
        for caller in frame.callers_in_cycle
            finish(caller)
        end
        for caller in frame.callers_in_cycle
            finalize_backedges(caller)
        end
    end

    nothing
end

max_stars_repo_name: BenChung/Julia-BTI

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
"""
Module for accessing GRIB files. Open a file with `GribFile(filename)`.

Exports
-------
Message, getbytes, writemessage, missingvalue, clone, data, maskedvalues,
Index, addfile!, keycount, select!, destroy, Nearest, find, findmultiple,
eachkey, keys, eachpoint, GribFile, nomultisupport
"""
module GRIB

using eccodes_jll

export nomultisupport

# Library constants
const CODES_NEAREST_SAME_GRID = 1 << 0
const CODES_NEAREST_SAME_DATA = 1 << 1
const CODES_NEAREST_SAME_POINT = 1 << 2

const CODES_TYPE_UNDEFINED = 0
const CODES_TYPE_LONG = 1
const CODES_TYPE_DOUBLE = 2
const CODES_TYPE_STRING = 3
const CODES_TYPE_BYTES = 4
const CODES_TYPE_SECTION = 5
const CODES_TYPE_LABEL = 6
const CODES_TYPE_MISSING = 7

const CODES_KEYS_ITERATOR_ALL_KEYS = 0

const CODES_END_OF_INDEX = -43

const errors = [
    "No error",		# 0 GRIB_SUCCESS
    "End of resource reached",		# -1 GRIB_END_OF_FILE
    "Internal error",		# -2 GRIB_INTERNAL_ERROR
    "Passed buffer is too small",		# -3 GRIB_BUFFER_TOO_SMALL
    "Function not yet implemented",		# -4 GRIB_NOT_IMPLEMENTED
    "Missing 7777 at end of message",		# -5 GRIB_7777_NOT_FOUND
    "Passed array is too small",		# -6 GRIB_ARRAY_TOO_SMALL
    "File not found",		# -7 GRIB_FILE_NOT_FOUND
    "Code not found in code table",		# -8 GRIB_CODE_NOT_FOUND_IN_TABLE
    "Array size mismatch",		# -9 GRIB_WRONG_ARRAY_SIZE
    "Key/value not found",		# -10 GRIB_NOT_FOUND
    "Input output problem",		# -11 GRIB_IO_PROBLEM
    "Message invalid",		# -12 GRIB_INVALID_MESSAGE
    "Decoding invalid",		# -13 GRIB_DECODING_ERROR
    "Encoding invalid",		# -14 GRIB_ENCODING_ERROR
    "Code cannot unpack because of string too small",		# -15 GRIB_NO_MORE_IN_SET
    "Problem with calculation of geographic attributes",		# -16 GRIB_GEOCALCULUS_PROBLEM
    "Memory allocation error",		# -17 GRIB_OUT_OF_MEMORY
    "Value is read only",		# -18 GRIB_READ_ONLY
    "Invalid argument",		# -19 GRIB_INVALID_ARGUMENT
    "Null handle",		# -20 GRIB_NULL_HANDLE
    "Invalid section number",		# -21 GRIB_INVALID_SECTION_NUMBER
    "Value cannot be missing",		# -22 GRIB_VALUE_CANNOT_BE_MISSING
    "Wrong message length",		# -23 GRIB_WRONG_LENGTH
    "Invalid key type",		# -24 GRIB_INVALID_TYPE
    "Unable to set step",		# -25 GRIB_WRONG_STEP
    "Wrong units for step (step must be integer)",		# -26 GRIB_WRONG_STEP_UNIT
    "Invalid file id",		# -27 GRIB_INVALID_FILE
    "Invalid grib id",		# -28 GRIB_INVALID_GRIB
    "Invalid index id",		# -29 GRIB_INVALID_INDEX
    "Invalid iterator id",		# -30 GRIB_INVALID_ITERATOR
    "Invalid keys iterator id",		# -31 GRIB_INVALID_KEYS_ITERATOR
    "Invalid nearest id",		# -32 GRIB_INVALID_NEAREST
    "Invalid order by",		# -33 GRIB_INVALID_ORDERBY
    "Missing a key from the fieldset",		# -34 GRIB_MISSING_KEY
    "The point is out of the grid area",		# -35 GRIB_OUT_OF_AREA
    "Concept no match",		# -36 GRIB_CONCEPT_NO_MATCH
    "Hash array no match",		# -37 GRIB_HASH_ARRAY_NO_MATCH
    "Definitions files not found",		# -38 GRIB_NO_DEFINITIONS
    "Wrong type while packing",		# -39 GRIB_WRONG_TYPE
    "End of resource",		# -40 GRIB_END
    "Unable to code a field without values",		# -41 GRIB_NO_VALUES
    "Grid description is wrong or inconsistent",		# -42 GRIB_WRONG_GRID
    "End of index reached",		# -43 GRIB_END_OF_INDEX
    "Null index",		# -44 GRIB_NULL_INDEX
    "End of resource reached when reading message",		# -45 GRIB_PREMATURE_END_OF_FILE
    "An internal array is too small",		# -46 GRIB_INTERNAL_ARRAY_TOO_SMALL
    "Message is too large for the current architecture",		# -47 GRIB_MESSAGE_TOO_LARGE
    "Constant field",		# -48 GRIB_CONSTANT_FIELD
    "Switch unable to find a matching case",		# -49 GRIB_SWITCH_NO_MATCH
    "Underflow",		# -50 GRIB_UNDERFLOW
    "Message malformed",		# -51 GRIB_MESSAGE_MALFORMED
    "Index is corrupted",		# -52 GRIB_CORRUPTED_INDEX
    "Invalid number of bits per value",		# -53 GRIB_INVALID_BPV
    "Edition of two messages is different",		# -54 GRIB_DIFFERENT_EDITION
    "Value is different",		# -55 GRIB_VALUE_DIFFERENT
    "Invalid key value",		# -56 GRIB_INVALID_KEY_VALUE
    "String is smaller than requested",		# -57 GRIB_STRING_TOO_SMALL
    "Wrong type conversion",		# -58 GRIB_WRONG_CONVERSION
    "Missing BUFR table entry for descriptor",		# -59 GRIB_MISSING_BUFR_ENTRY
    "Null pointer",		# -60 GRIB_NULL_POINTER
    "Attribute is already present, cannot add",		# -61 GRIB_ATTRIBUTE_CLASH
    "Too many attributes. Increase MAX_ACCESSOR_ATTRIBUTES",		# -62 GRIB_TOO_MANY_ATTRIBUTES
    "Attribute not found.",		# -63 GRIB_ATTRIBUTE_NOT_FOUND
    "Edition not supported.",		# -64 GRIB_UNSUPPORTED_EDITION
    "Value out of coding range",		# -65 GRIB_OUT_OF_RANGE
    "Size of bitmap is incorrect",		# -66 GRIB_WRONG_BITMAP_SIZE
    "Value mismatch",		# 1 GRIB_VALUE_MISMATCH
    "double values are different",		# 2 GRIB_DOUBLE_VALUE_MISMATCH
    "long values are different",		# 3 GRIB_LONG_VALUE_MISMATCH
    "byte values are different",		# 4 GRIB_BYTE_VALUE_MISMATCH
    "string values are different",		# 5 GRIB_STRING_VALUE_MISMATCH
    "Offset mismatch",		# 6 GRIB_OFFSET_MISMATCH
    "Count mismatch",		# 7 GRIB_COUNT_MISMATCH
    "Name mismatch",		# 8 GRIB_NAME_MISMATCH
    "Type mismatch",		# 9 GRIB_TYPE_MISMATCH
    "Type and value mismatch",		# 10 GRIB_TYPE_AND_VALUE_MISMATCH
    "Unable to compare accessors",		# 11 GRIB_UNABLE_TO_COMPARE_ACCESSORS
    "Unable to reset iterator",		# 12 GRIB_UNABLE_TO_RESET_ITERATOR
    "Assertion failure",		# 13 GRIB_ASSERTION_FAILURE
]

@inline function errorcheck(code)
    if code != 0
        throw(ErrorException("Error: $(errors[-code+1])"))
    end
end

mutable struct codes_context
end

# Turn on multi-field support by default
function __init__()
    if Sys.iswindows()
        @warn "You are using GRIB.jl on Windows which is currently is not supported. You will likely hit a ReadOnlyMemoryError, and, if you know how to get around it, please let us know!"
    end

    # Let library know where its own shared data is
    share = joinpath(dirname(eccodes_jll.eccodes_path), "..", "share", "eccodes")
    ENV["ECCODES_DEFINITION_PATH"] = joinpath(Base.@__DIR__, "..", "definitions") * ":" * joinpath(share, "definitions")
    ENV["ECCODES_SAMPLES_PATH"] = joinpath(share, "samples")

    # Turn on multi-field support by default for straightforward use of NCEP files
    ccall((:codes_grib_multi_support_on, eccodes), Cvoid, (Ptr{codes_context},), C_NULL)
end

"""
    nomultisupport

Turn off multi-field message support. Recommended if you are sure no files you will
work with in this session will have multi-field messages.
"""
function nomultisupport()
    ccall((:codes_grib_multi_support_off, eccodes), Cvoid, (Ptr{codes_context},), C_NULL)
end

include("index.jl")
include("gribfile.jl")
include("message.jl")
include("nearest.jl")
include("messageiters.jl")

end


text:
using Pkg
Pkg.add("assert")


text:
function gui_init_all()
	global window=GtkWindow("Foam World",576,512;resizable=false,visible=false)
	signal_connect(gui_leave_all,window,"destroy")
end
function gtk_leave_all()
end
# 工具
RGB_(r::UInt8,g::UInt8,b::UInt8)=RGB(reinterpret(N0f8,r),reinterpret(N0f8,g),reinterpret(N0f8,b))
gtk_showstr(s::String)=GtkTextView(;editable=false,buffer=GtkTextBuffer(;text=s))
# 主菜单
function gui_主菜单_init()
	global winbox=GtkBox(:v)
	push!(window,winbox)
	lc=langs[:choice]
	but_set=GtkButton(lc[:to_setting]);push!(winbox,but_set)
	but_new=GtkButton(lc[:to_newgame]);push!(winbox,but_new)
	but_open=GtkButton(lc[:to_open]);push!(winbox,but_open)
	signal_connect(gui_主菜单_click_set,but_set,"clicked")
	signal_connect(gui_主菜单_click_new,but_new,"clicked")
	signal_connect(gui_主菜单_click_open,but_open,"clicked")
end
function gui_主菜单_click_set()
	global curpage=:设置
	gui_设置_init()
end
# setting
function gui_设置_init()
	set_username_box=Box(:h)
end
# game
function info_clear()
	global infobox=GtkBox(:h)
end
function info_log(s::String,color::RGB=RGB(0,0,0))
	push!(infobox,gtk_showstr(s))
end
function info_error(s::Symbol...;extra::String="")
	t=langs[:error][s...]
	if extra!=""
		t*=": $extra"
	end
	info_log(t,RGB(1,0,0))
end
function info_help(s::Symbol...;extra::String="")
	if gsetting[:show_help]
		t=langs[:help][s...]
		if extra!=""
			t*=": $extra"
		end
		info_log(t,RGB(1,1,0))
	end
end


text:
@testset "Day8.jl" begin
    # Write your tests here.
    lines = read_infile(Day8_1, "$(@__DIR__)/assets/day8_test1.txt")
    @test solve(Day8_1, lines) == 5

    @test solve(Day8_2, lines) == 8
end


text:
using BinaryProvider

# This is where all binaries will get installed
const prefix = Prefix(!isempty(ARGS) ? ARGS[1] : joinpath(@__DIR__,"usr"))
libsnappy = LibraryProduct(prefix, "libsnappy", :libsnappy)
products = [libsnappy]

# Download binaries from hosted location
bin_prefix = "https://github.com/davidanthoff/SnappyBuilder/releases/download/v1.1.7%2Bbuild.1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64) => ("$bin_prefix/SnappyBuilder.aarch64-linux-gnu.tar.gz", "8cb4b9af85da7efa20bc4787ea16af4b76123d790590bfd447cbe092f123f2f8"),
    Linux(:armv7l) => ("$bin_prefix/SnappyBuilder.arm-linux-gnueabihf.tar.gz", "67469c1ce2391e65c3d9f5f8162f0e14f69ce0e19af0389ad8d77226bf79a922"),
    Linux(:i686) => ("$bin_prefix/SnappyBuilder.i686-linux-gnu.tar.gz", "f34512ca97f500cdaeeea19aeb3e3075b3aa2b4342cf3a5a3788c77eb634f0bd"),
    Windows(:i686) => ("$bin_prefix/SnappyBuilder.i686-w64-mingw32.tar.gz", "53a28428e5c3d75bef1b2ce90b8af8905177b7dae313c8a859afe24c7c9571dd"),
    Linux(:powerpc64le) => ("$bin_prefix/SnappyBuilder.powerpc64le-linux-gnu.tar.gz", "32ecdd9185c2f39e2c95db9e8adb2adf8f251d65c0a8e43c0ee6ec46778045cb"),
    MacOS() => ("$bin_prefix/SnappyBuilder.x86_64-apple-darwin14.tar.gz", "8fe643e462799ce141b4d306427f73b2a373651282cdfabcaeb976e705069b7c"),
    Linux(:x86_64) => ("$bin_prefix/SnappyBuilder.x86_64-linux-gnu.tar.gz", "9a45685503de0db44068126d67f49e241f9013970ddc560cfd7ba17b411dac8c"),
    Windows(:x86_64) => ("$bin_prefix/SnappyBuilder.x86_64-w64-mingw32.tar.gz", "1c13a31c23225b812ee477dffa1fc77d6f501c0a3b65a9920edf5aa4a76f4475"),
)

# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=true) for p in products)
    if haskey(download_info, platform_key())
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
    else
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something more even more ambitious here.
        error("Your platform $(triplet(platform_key())) is not supported by this package!")
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
function uniform_grid(N, M, l₀=1.0)
    u = zeros(Float64, (4,N,M))

    for i = 1:N, j=1:M
        u[3,i,j] = j*l₀
        u[4,i,j] = -i*l₀
    end

    return u
end


text:
push!(LOAD_PATH,joinpath(@__DIR__, "../src"))

using Documenter, InstrumentOperator, Literate, Plots
using DocumenterTools: Themes
using Pkg
Pkg.add(url="https://github.com/RadiativeTransfer/RadiativeTransfer.jl")

push!(LOAD_PATH,"./src/")
ENV["GKSwstype"] = "nul"

function build()

    tutorials = ["CrossSection_convolution.jl"] # ,
    tutorials_paths = [joinpath(@__DIR__, "src", "pages", "tutorials", tutorial) for tutorial in tutorials]

    for tutorial in tutorials_paths
        Literate.markdown(tutorial, joinpath(@__DIR__, "src", "pages", "tutorials"))
    end

    tutorials_md = [joinpath("pages", "tutorials", tutorial[1:end-3]) * ".md" for tutorial in tutorials]

    pages = Any[
            "Home"                  => "index.md",
            "Tutorials"             => tutorials_md
        ]

    mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict(),
            ),
        ))

    format = Documenter.HTML(
        #assets = [
        #    asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
        #    ],
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = mathengine,
        collapselevel = 1,
        )
    # download the themes
    for file in ("radiativetransfer-lightdefs.scss", "radiativetransfer-darkdefs.scss", "radiativetransfer-style.scss")
        download("https://raw.githubusercontent.com/RadiativeTransfer/doctheme/main/$file", joinpath(@__DIR__, file))
    end
    # create the themes
    for w in ("light", "dark")
        header = read(joinpath(@__DIR__, "radiativetransfer-style.scss"), String)
        theme = read(joinpath(@__DIR__, "radiativetransfer-$(w)defs.scss"), String)
        write(joinpath(@__DIR__, "radiativetransfer-$(w).scss"), header*"\n"*theme)
    end
    # compile the themes
    Themes.compile(joinpath(@__DIR__, "radiativetransfer-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
    Themes.compile(joinpath(@__DIR__, "radiativetransfer-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

    makedocs(
    sitename = "Instrument Operator",
    format = format,
    clean = false,
    modules = [InstrumentOperator],
    pages = pages
)
end

build()

deploydocs(
    repo = "github.com/RadiativeTransfer/InstrumentOperator.jl.git",
    target = "build",
    devbranch = "main",
    push_preview = true
)


text:
#######################################################################################
#                                   Load Packages
#######################################################################################
cd(@__DIR__)
using Pkg
Pkg.activate("../..")
using StatsPlots, Random, ACTRModels, CSV, FFTDists, DifferentialEvolutionMCMC
include("Count.jl")
Random.seed!(761321)
#######################################################################################
#                                   Generate Data
#######################################################################################
n_trials = 50
s = 0.3
blc = 1.5
data = simulate(n_trials; s, blc)
#######################################################################################
#                                    Define Model
#######################################################################################
priors = (
  blc=(Normal(1.5, 1),),
  s=(Truncated(Normal(.3, .5), 0.0, Inf),)
)
bounds = ((-Inf,Inf),(eps(),Inf))
model = DEModel(;priors, model=loglike, data)
de = DE(;priors, bounds, burnin=1000)
n_iter = 2000
#######################################################################################
#                                 Estimate Parameters
#######################################################################################
chain = sample(model, de, MCMCThreads(), n_iter, progress=true)
#######################################################################################
#                                         Plot
#######################################################################################
pyplot()
post_plot = plot(chain, xaxis=font(5), yaxis=font(5), seriestype=(:pooleddensity),
  grid=false, titlefont=font(5), size=(200,200), color=:gray)
plot!(title = "", xlabel="blc", subplot=1)
plot!(title = "", xlabel="s", subplot=2)
#######################################################################################
#                                  Posterior Predictive
#######################################################################################
preds = posterior_predictive(x->simulate(n_trials; x...), chain, 1000)
rts = vcat(preds...)
post_pred = histogram(rts, xlabel = "RT", ylabel="Frequency", xaxis=font(7), yaxis=font(7),
    grid=false, color=:grey, leg=false, size=(300,150), titlefont=font(7),
    xlims=(.5,2.5))
savefig(post_plot, "Count_posteriors.eps")
savefig(post_pred, "Count_postpred.eps")


text:
# Check consistency with the original WORLD
# Ground truth data was genearted using `generate_groundtruth.sh` in
# https://github.com/r9y9/WORLD
# Note that all of the results of the original WORLD were dumped as %.16lf
# accuracy on Ubuntu 14.04 64bit machine. g++ v4.9.2 was used to compile the
# WORLD.

using WORLD
using Test, DelimitedFiles

@assert WORLD.version >= v"0.3.0"

@info("Check consistency with the original WORLD")

x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "x.txt")))

fs = 22050
period = 5.0
opt = HarvestOption(40.0, 800.0, period)

# Fundamental frequency (f0) estimation by DIO
f0, timeaxis = harvest(x, fs, opt)
f0_org = vec(readdlm(joinpath(dirname(@__FILE__), "data", "f0.txt")))

println("Maximum error in Harvest is $(maximum(abs.(f0-f0_org)))")
@test length(f0) == length(f0_org)
@test isapprox(f0, f0_org, atol=1.0e-10)

# Spectral envelope estimation by CheapTrick
spectrogram = cheaptrick(x, fs, timeaxis, f0; opt=CheapTrickOption(fs))
spectrogram_org = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram.txt"))'

println("Maximum error in CheapTrick is $(maximum(abs.(spectrogram - spectrogram_org)))")
@test size(spectrogram) == size(spectrogram_org)
@test isapprox(spectrogram, spectrogram_org, atol=1.0e-9)

aperiodicity = d4c(x, fs, timeaxis, f0; opt=D4COption())
aperiodicity_org = readdlm(joinpath(dirname(@__FILE__), "data", "aperiodicity.txt"))'
println("Maximum error in D4C is $(maximum(abs.(aperiodicity-aperiodicity_org)))")
@test size(aperiodicity) == size(aperiodicity)
@test isapprox(aperiodicity, aperiodicity_org, atol=1.0e-6)

# Synthesis
y_length = trunc(Int, ((length(f0)-1)*period/1000 * fs) + 1)
y = synthesis(f0, spectrogram, aperiodicity, period, fs, y_length)
y_org = vec(readdlm(joinpath(dirname(@__FILE__), "data", "x_synthesized.txt")))
println("Maximum error in Synthesis is $(maximum(abs.(y-y_org)))")
@test length(y) == length(y_org)
@test isapprox(y, y_org, atol=1.0e-10)