text:
@testset "Stateful max-min filters" begin

    len = 1000
    Random.seed!(0)
    T = Float64
    vec = rand(T, len)
    window = 15

    mx_ref = similar(vec)
    mn_ref = similar(vec)

    for i = 1:len
        mx_ref[i] = maximum(vec[max(1, i-window+1) : i])
        mn_ref[i] = minimum(vec[max(1, i-window+1) : i])
    end
    range_ref = mx_ref .- mn_ref

    @testset "MaxMinFilter" begin
        f = MaxMinFilter{T}(window)
        mx, mn = filt(f, vec)
        @test mx == mx_ref
        @test mn == mn_ref
        # reset state
        reset!(f)
        mx, mn = filt(f, vec[1:len÷2])
        @test mx == mx_ref[1:len÷2]
        @test mn == mn_ref[1:len÷2]
        # maintain state
        mx, mn = filt(f, vec[len÷2+1:end])
        @test mx == mx_ref[len÷2+1:end]
        @test mn == mn_ref[len÷2+1:end]
    end

    @testset "RangeFilter" begin
        f = RangeFilter{T}(window)
        range = filt(f, vec)
        @test range == range_ref
        # reset state
        reset!(f)
        range = filt(f, vec[1:len÷2])
        @test range == range_ref[1:len÷2]
        # maintain state
        range = filt(f, vec[len÷2+1:end])
        @test range == range_ref[len÷2+1:end]
    end

    @testset "MaxFilter" begin
        f = MaxFilter{T}(window)
        mx = filt(f, vec)
        @test mx == mx_ref
        # reset state
        reset!(f)
        mx = filt(f, vec[1:len÷2])
        @test mx == mx_ref[1:len÷2]
        # maintain state
        mx = filt(f, vec[len÷2+1:end])
        @test mx == mx_ref[len÷2+1:end]
    end

    @testset "MinFilter" begin
        f = MinFilter{T}(window)
        mn = filt(f, vec)
        @test mn == mn_ref
        # reset state
        reset!(f)
        mn = filt(f, vec[1:len÷2])
        @test mn == mn_ref[1:len÷2]
        # maintain state
        mn = filt(f, vec[len÷2+1:end])
        @test mn == mn_ref[len÷2+1:end]
    end

end


text:
using StatsBase: entropy

if get_stage() in ("Others", "all")
    @testset "TuringWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dW1 = Wishart(dim + 4, A)
        dW2 = TuringWishart(dim + 4, A)
        dW3 = TuringWishart(dW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, meanlogdet, entropy, cov, var)
            @test F(dW1) == F(dW2) == F(dW3)
        end
        @test Matrix(mode(dW1)) == mode(dW2) == mode(dW3)
        xw = rand(dW2)
        @test insupport(dW1, xw)
        @test insupport(dW2, xw)
        @test insupport(dW3, xw)
        @test logpdf(dW1, xw) == logpdf(dW2, xw) == logpdf(dW3, xw)
    end

    @testset "TuringInverseWishart" begin
        dim = 3
        A = Matrix{Float64}(I, dim, dim)
        dIW1 = InverseWishart(dim + 4, A)
        dIW2 = TuringInverseWishart(dim + 4, A)
        dIW3 = TuringInverseWishart(dIW1)
        mean = Distributions.mean
        @testset "$F" for F in (size, rank, mean, cov, var)
            @test F(dIW1) == F(dIW2) == F(dIW3)
        end
        @test Matrix(mode(dIW1)) == mode(dIW2) == mode(dIW3)
        xiw = rand(dIW2)
        @test insupport(dIW1, xiw)
        @test insupport(dIW2, xiw)
        @test insupport(dIW3, xiw)
        @test logpdf(dIW1, xiw) == logpdf(dIW2, xiw) == logpdf(dIW3, xiw)
    end

    @testset "TuringMvNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                C = Matrix{Float64}(I, 3, 3)
                d1 = TuringMvNormal(m, C)
            elseif TD <: TuringDiagMvNormal
                C = ones(3)
                d1 = TuringMvNormal(m, C)
            else
                C = 1.0
                d1 = TuringMvNormal(m, C)
            end
            d2 = MvNormal(m, C)

            @testset "$F" for F in (length, size)
                @test F(d1) == F(d2)
            end

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test isapprox(logpdf(d1, x1), logpdf(d2, x1), rtol = 1e-6)
            @test isapprox(logpdf(d1, x2), logpdf(d2, x2), rtol = 1e-6)
        end
    end

    @testset "TuringMvLogNormal" begin
        @testset "$TD" for TD in [TuringDenseMvNormal, TuringDiagMvNormal, TuringScalMvNormal]
            m = rand(3)
            if TD <: TuringDenseMvNormal
                C = Matrix{Float64}(I, 3, 3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            elseif TD <: TuringDiagMvNormal
                C = ones(3)
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            else
                C = 1.0
                d1 = TuringMvLogNormal(TuringMvNormal(m, C))
            end
            d2 = MvLogNormal(MvNormal(m, C))

            @test length(d1) == length(d2)

            x1 = rand(d1)
            x2 = rand(d1, 3)
            @test isapprox(logpdf(d1, x1), logpdf(d2, x1), rtol = 1e-6)
            @test isapprox(logpdf(d1, x2), logpdf(d2, x2), rtol = 1e-6)

            x2[:, 1] .= -1
            @test isinf(logpdf(d1, x2)[1])
            @test isinf(logpdf(d2, x2)[1])
        end
    end

    @testset "TuringUniform" begin
        @test logpdf(TuringUniform(), param(0.5)) == 0
    end

    @testset "Semicircle" begin
        @test Tracker.data(logpdf(Semicircle(1.0), param(0.5))) == logpdf(Semicircle(1.0), 0.5)
    end

    @testset "TuringPoissonBinomial" begin
        d1 = TuringPoissonBinomial([0.5, 0.5])
        d2 = PoissonBinomial([0.5, 0.5])
        @test quantile(d1, 0.5) == quantile(d2, 0.5)
        @test minimum(d1) == minimum(d2)
    end

    @testset "Inverse of pi" begin
        @test 1/pi == inv(pi)
    end

    @testset "Others" begin
        @test fill(param(1.0), 3) isa TrackedArray
        x = rand(3)
        @test isapprox(Tracker.data(Tracker.gradient(logsumexp, x)[1]),
            ForwardDiff.gradient(logsumexp, x), atol = 1e-5)
        A = rand(3, 3)'; A = A + A' + 3I;
        C = cholesky(A; check = true)
        factors, info = DistributionsAD.turing_chol(A, true)
        @test factors == C.factors
        @test info == C.info
        B = copy(A)
        @test DistributionsAD.zygote_ldiv(A, B) == A \ B
    end

    @testset "Entropy" begin
        sigmas = exp.(randn(10))
        d1 = TuringDiagMvNormal(zeros(10), sigmas)
        d2 = MvNormal(zeros(10), sigmas)

        @test isapprox(entropy(d1), entropy(d2), rtol = 1e-6)
    end

    @testset "Params" begin
        m = rand(10)
        sigmas = randexp(10)

        d = TuringDiagMvNormal(m, sigmas)
        @test params(d) == (m, sigmas)

        d = TuringScalMvNormal(m, sigmas[1])
        @test params(d) == (m, sigmas[1])
    end
end

ext: jl


text:
# This file was generated, do not modify it. # hide
mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
@show round.(fp.coefs[1:3], sigdigits=3)
@show round(fp.intercept, sigdigits=3)


text:
using Oceananigans.Grids: topology
using CUDA

include("utils_for_runtests.jl")
include("data_dependencies.jl")

archs = test_architectures()

function get_fields_from_checkpoint(filename)
    file = jldopen(filename)

    tracers = keys(file["tracers"])
    tracers = Tuple(Symbol(c) for c in tracers)

    velocity_fields = (u = file["velocities/u/data"],
                       v = file["velocities/v/data"],
                       w = file["velocities/w/data"])

    tracer_fields =
        NamedTuple{tracers}(Tuple(file["tracers/$c/data"] for c in tracers))

    current_tendency_velocity_fields = (u = file["timestepper/Gⁿ/u/data"],
                                        v = file["timestepper/Gⁿ/v/data"],
                                        w = file["timestepper/Gⁿ/w/data"])

    current_tendency_tracer_fields =
        NamedTuple{tracers}(Tuple(file["timestepper/Gⁿ/$c/data"] for c in tracers))

    previous_tendency_velocity_fields = (u = file["timestepper/G⁻/u/data"],
                                         v = file["timestepper/G⁻/v/data"],
                                         w = file["timestepper/G⁻/w/data"])

    previous_tendency_tracer_fields =
        NamedTuple{tracers}(Tuple(file["timestepper/G⁻/$c/data"] for c in tracers))

    close(file)

    solution = merge(velocity_fields, tracer_fields)
    Gⁿ = merge(current_tendency_velocity_fields, current_tendency_tracer_fields)
    G⁻ = merge(previous_tendency_velocity_fields, previous_tendency_tracer_fields)

    return solution, Gⁿ, G⁻
end

include("regression_tests/thermal_bubble_regression_test.jl")
include("regression_tests/rayleigh_benard_regression_test.jl")
include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")
include("regression_tests/hydrostatic_free_turbulence_regression_test.jl")

@testset "Regression" begin
    @info "Running regression tests..."

    for arch in archs
        for grid_type in [:regular, :vertically_unstretched]
            @testset "Thermal bubble [$(typeof(arch)), $grid_type grid]" begin
                @info "  Testing thermal bubble regression [$(typeof(arch)), $grid_type grid]"
                run_thermal_bubble_regression_test(arch, grid_type)
            end

            @testset "Rayleigh–Bénard tracer [$(typeof(arch)), $grid_type grid]]" begin
                @info "  Testing Rayleigh–Bénard tracer regression [$(typeof(arch)), $grid_type grid]"
                run_rayleigh_benard_regression_test(arch, grid_type)
            end

            for closure in (AnisotropicMinimumDissipation(ν=1.05e-6, κ=1.46e-7), SmagorinskyLilly(C=0.23, Cb=1, Pr=1, ν=1.05e-6, κ=1.46e-7))
                closurename = string(typeof(closure).name.wrapper)
                @testset "Ocean large eddy simulation [$(typeof(arch)), $closurename, $grid_type grid]" begin
                    @info "  Testing oceanic large eddy simulation regression [$(typeof(arch)), $closurename, $grid_type grid]"
                    run_ocean_large_eddy_simulation_regression_test(arch, grid_type, closure)
                end
            end
        end

        explicit_free_surface = ExplicitFreeSurface(gravitational_acceleration=1.0)
        implicit_free_surface = ImplicitFreeSurface(gravitational_acceleration = 1.0,
                                                    solver_method = :PreconditionedConjugateGradient,
                                                    tolerance = 1e-15)

        x_bounded_lat_lon_grid  = RegularLatitudeLongitudeGrid(size = (160, 60, 3),
                                                               longitude = (-160, 160),
                                                               latitude = (-60, 60),
                                                               z = (-90, 0),
                                                               halo = (2, 2, 2))

        x_periodic_lat_lon_grid  = RegularLatitudeLongitudeGrid(size = (180, 60, 3),
                                                                longitude = (-180, 180),
                                                                latitude = (-60, 60),
                                                                z = (-90, 0),
                                                                halo = (2, 2, 2))

        for grid in [x_bounded_lat_lon_grid, x_periodic_lat_lon_grid]
            for free_surface in [explicit_free_surface, implicit_free_surface]

                # GPU + ExplicitFreeSurface is broken. See:
                # https://github.com/CliMA/Oceananigans.jl/pull/1985
                # if !(arch isa GPU && topology(grid, 1) === Periodic && free_surface isa ExplicitFreeSurface)

                free_surface_str = string(typeof(free_surface).name.wrapper)

                @testset "Hydrostatic free turbulence regression [$(typeof(arch)), $(topology(grid, 1)) longitude, $free_surface_str]" begin
                    @info "  Testing Hydrostatic free turbulence [$(typeof(arch)), $(topology(grid, 1)) longitude, $free_surface_str]"
                    run_hydrostatic_free_turbulence_regression_test(grid, free_surface, arch)
                end
            end
	    end
	end
end


text:
Vector{UInt8}(
"24\r\n" *
"'Twas brillig, and the slithy toves\n" *
"27\r\n" *
"      Did gyre and gimble in the wabe:\n" *
"1E\r\n" *
"All mimsy were the borogoves,\n" *
"22\r\n" *
"      And the mome raths outgrabe." *
"0\r\n" *
"\r\n"
)


text:
## Inplace definition (the `du` array is passed to the function and changed), the default way to
## define this in Julia.
function lv_comp(du, u, p, t)
    @unpack r1, r2, α12, α21, K1, K2 = p
    du[1] = u[1] * r1 * (1 - (u[1] + α12 * u[2]) / K1)
    du[2] = u[2] * r2 * (1 - (u[2] + α21 * u[1]) / K2)
    return
end

## Make a version that allocates the output `du`, useful for symbolic calculations
function lv_comp(u, p, t)
    du = similar(u)
    lv_comp(du, u, p, t)
    return du
end

# # Model Parameters

@with_kw mutable struct LVPar
    α12 = 0.8
    α21 = 0.3
    r1 = 1.0
    r2 = 2.0
    K1 = 1.5
    K2 = 1.0
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using Neuroimaging, Test, BDF

@testset "GeneralEEG" begin

    fname = joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf")

    s = read_EEG(fname)

    @test isa(s, NeuroimagingMeasurement)
    @test isa(s, EEG)

    @testset "Show" begin
        show(s)
    end

    @testset "Read file" begin

        s = read_EEG(fname)

        @info samplingrate(s)
        @test samplingrate(s) == 8192.0
        @info samplingrate(s)
        @test samplingrate(Int, s) == 8192
        @test isa(samplingrate(s), AbstractFloat) == true
        @test isa(samplingrate(Int, s), Int) == true

        @test isapprox(maximum(s.data[:, 2]), 54.5939; atol = 0.1)
        @test isapprox(minimum(s.data[:, 4]), -175.2503; atol = 0.1)

        s = read_EEG(fname, valid_triggers = [8])
        s = read_EEG(fname, min_epoch_length = 8388)
        s = read_EEG(fname, max_epoch_length = 8389)

        s = read_EEG(fname)
        original_length = length(s.triggers["Index"])
        s = read_EEG(fname, remove_first = 10)
        @test original_length - 10 == length(s.triggers["Index"])

        s = read_EEG(fname, max_epochs = 12)
        @test length(s.triggers["Index"]) == 12

        s = read_EEG(fname)
        s.triggers = extra_triggers(s.triggers, 1, 7, 0.7, samplingrate(s))
        @test length(s.triggers["Index"]) == 56
    end

    @testset "Sensors" begin

        s1 = sensors(deepcopy(s))
        @test length(s1) == 6
        @test length(x(s1)) == 6
        @test length(labels(s1)) == 6

    end

    @testset "Channel names" begin

        s1 = deepcopy(s)
        s1 = channelnames(s1, 1, "A01")
        s1 = channelnames(s1, 2, "A05")
        s1 = channelnames(s1, 3, "A11")
        s1 = channelnames(s1, 4, "B03")
        s1 = channelnames(s1, 5, "A17")
        s1 = channelnames(s1, 6, "B17")
        @test channelnames(s1) == ["A01", "A05", "A11", "B03", "A17", "B17"]

    end

    @testset "Triggers" begin

        dats, evtTab, trigs, statusChan = readBDF(fname)
        sampRate = readBDFHeader(fname)["sampRate"][1]

        @test trigs == create_channel(
            evtTab,
            dats,
            sampRate,
            code = "code",
            index = "idx",
            duration = "dur",
        )
        @test trigs ==
              trigger_channel(read_EEG(fname, valid_triggers = collect(-1000:10000)))

        # Convert between events and channels

        dats, evtTab, trigs, statusChan = readBDF(fname)
        events = create_events(trigs, sampRate)
        channel = create_channel(events, dats, sampRate)

        @test size(channel) == size(trigs)
        @test channel == trigs

    end

end


text:
#JuliaRL/src/display.jl

using ReinforcementLearning
using PyCall
using Random
using ReinforcementLearningEnvironments
using Flux
using Logging

include("./shared.jl")

"""
    display_agent(save_dir, duration::Int = 1000)

Runs the agent saved at `save_dir`, rendering the environment,
up to `duration` steps.
"""
function display_agent(save_dir, duration::Int = 5000)
    env = LunarLander()
    agent = RLCore.load(save_dir, Agent)
    Flux.testmode!(agent)

    stop_condition = StopAfterStep(duration)
    disp_hook = DoEveryNStep(1) do t, agent, env
        env.env.pyenv.render()
    end

    disp_hook = DoEveryNStep(1) do t, agent, env
        env.env.pyenv.render()
    end

    print_hook = DoEveryNEpisode(1) do t, agent, env
        println("- Ep N $t")
    end

    hook = ComposedHook(disp_hook, print_hook)

    run(agent, env, stop_condition, hook)
end

max_stars_repo_name: FeLusiani/JuliaRL

max_issues_repo_name: FeLusiani/JuliaRL

max_forks_repo_name: FeLusiani/JuliaRL


text:
module DerivableFunctions

using Reexport
using ReverseDiff, Zygote, FiniteDiff, FiniteDifferences

@reexport using DerivableFunctionsBase

# import such that they are available downstream as if defined here
import DerivableFunctionsBase: MaximalNumberOfArguments, KillAfter, Builder
import DerivableFunctionsBase: _GetArgLength, _GetArgLengthOutOfPlace, _GetArgLengthInPlace
import DerivableFunctionsBase: GetSymbolicDerivative, SymbolicPassthrough

import DerivableFunctionsBase: _GetDeriv, _GetGrad, _GetJac, _GetHess, _GetMatrixJac, _GetDoubleJac
import DerivableFunctionsBase: _GetDerivPass, _GetGradPass, _GetJacPass, _GetHessPass, _GetDoubleJacPass, _GetMatrixJacPass
import DerivableFunctionsBase: _GetGrad!, _GetJac!, _GetHess!, _GetMatrixJac!
import DerivableFunctionsBase: _GetGradPass!, _GetJacPass!, _GetHessPass!, _GetMatrixJacPass!

import DerivableFunctionsBase: suff
suff(x::ReverseDiff.TrackedReal) = typeof(x)

## Add new backends to the output of diff_backends()
import DerivableFunctionsBase: AddedBackEnds
AddedBackEnds(::Val{true}) = [:ReverseDiff, :Zygote, :FiniteDifferences, :FiniteDiff]

include("DifferentiationOperators.jl")

end


text:
module MusicProcessing

using DSP
using FFTW
using FixedPointNumbers
using Requires
using Unitful
using IntervalSets
using Statistics
using Unitful: ns, ms, µs, s, Hz, kHz, MHz, GHz, THz
using LinearAlgebra:mul!

# types used for fixed-point 16-bit and 32-bit encoding
const PCM16Sample = Fixed{Int16, 15}
const PCM32Sample = Fixed{Int32, 31}

const Seconds = Unitful.Quantity{Int64, Unitful.𝐓,Unitful.FreeUnits{(s,),Unitful.𝐓,nothing}}
const Hertz = Unitful.Quantity{Int64,Unitful.𝐓^-1,Unitful.FreeUnits{(Hz,),Unitful.𝐓^-1,nothing}}
export Hz, kHz, s, ..
export PCM16Sample, PCM32Sample

include("util.jl")
include("SampleBuf.jl")

include("complex.jl")

include("audio.jl")
include("TFR.jl")

include("mel.jl")
include("constantq.jl")
include("chroma.jl")

export AbstractSampleBuf, SampleBuf, SpectrumBuf

function __init__()
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("display.jl")
end

end # module


text:
struct All
end
struct Matched
end

function smove_forward!(G, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
    nhd = neighbours(G, i)
    for i in nhd
        t[i], x[i] = t′, x[i] + θ[i]*(t′ - t[i])
    end
    t, x, θ
end
function smove_forward!(i::Int, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
    t[i], x[i] = t′, x[i] + θ[i]*(t′ - t[i])
    return t, x, θ
end
neighbours(::Nothing, i) = Int[]
smove_forward!(::Nothing, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag}) = t, x, θ
smove_forward!(::All, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag}) = smove_forward!(t, x, θ, t′, Z)
smove_forward!(::Nothing, i, t, x, θ, t′, Z::Union{Boomerang, FactBoomerang}) = t, x, θ
smove_forward!(::All, i, t, x, θ, t′, Z::Union{Boomerang, FactBoomerang}) = smove_forward!(t, x, θ, t′, Z)

function smove_forward!(t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
    for i in eachindex(x)
        t[i], x[i] = t′, x[i] + θ[i]*(t′ - t[i])
    end
    t, x, θ
end
function smove_forward!(G, i, t, x, θ, t′, B::Union{Boomerang, FactBoomerang})
    nhd = neighbours(G, i)
    for i in nhd
        τ = t′ - t[i]
        s, c = sincos(τ)
        t[i], x[i], θ[i] = t′, (x[i] - B.μ[i])*c + θ[i]*s + B.μ[i],
                    -(x[i] - B.μ[i])*s + θ[i]*c
    end
    t, x, θ
end

function smove_forward!(t, x, θ, t′, B::Union{Boomerang, FactBoomerang})
    for i in eachindex(x)
        τ = t′ - t[i]
        s, c = sincos(τ)
        t[i], x[i], θ[i] = t′, (x[i] - B.μ[i])*c + θ[i]*s + B.μ[i],
                    -(x[i] - B.μ[i])*s + θ[i]*c
    end
    t, x, θ
end

function event(i, t::Vector, x, θ, Z::Union{ZigZag,FactBoomerang,JointFlow})
    t[i], i, x[i], θ[i]
end

"""
    SelfMoving()

Indicates as `args[1]` that `∇ϕ` depends only on few coeffients
and takes responsibility to call `smove_forward!`.

Replaced by `ExtendedForm`.
"""
const SelfMoving = ExtendedForm
export ExtendedForm, SelfMoving

∇ϕ_(∇ϕ, t, x, θ, i, t′, Z, args...) = ∇ϕ(x, i, args...)
∇ϕ_(∇ϕ, t, x, θ, i, t′, Z, S::SelfMoving, args...) = ∇ϕ(t, x, θ, i, t′, Z, args...)
sλ(∇ϕi, i, x, θ, Z::Union{ZigZag,FactBoomerang,JointFlow}) = λ(∇ϕi, i, x, θ, Z)
sλ̄(abc, Δt) = pos(abc[1] + abc[2]*Δt)

"""
    spdmp_inner!(rng, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, (acc, num),
F::Union{ZigZag,FactBoomerang}, args...; factor=1.5, adapt=false, adaptscale=false)

[Outdated]
Inner loop of the factorised samplers: the factorised Boomerang algorithm and
the Zig-Zag sampler. Given a dependency graph `G`, gradient `∇ϕ`,
current position `x`, velocity `θ`, Queue of events `Q`, time `t`, tuning parameter `c`,
terms of the affine bounds `a`,`b` and time when the upper bounds were computed `t_old`

The sampler 1) extracts from the queue the first event time. 2) moves deterministically
according to its dynamics until event time. 3) Evaluates whether the event
time is a accepted reflection or refreshment time or shadow time. 4) If it is a
reflection time, the velocity reflects according its reflection rule, if it is a
refreshment time, the sampler updates the velocity from its prior distribution (Gaussian).
In both cases, updates `Q` according to the dependency graph `G`. The sampler proceeds
until the next accepted reflection time or refreshment time. `(num, acc)`
incrementally counts how many event times occour and how many of those are real reflection times.
"""
function spdmp_inner!(rng, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, (acc, num),
     F::Union{ZigZag,FactBoomerang}, args...; factor=1.5, adapt=false, adaptscale=false)
    n = length(x)
    while true
        i, t′ = peek(Q)
        refresh = i > n
        if refresh
            i = rand(1:n)
        end
        t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
        if refresh
            i = rand(1:n)
            t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
            if adaptscale && F isa ZigZag
                adapt_γ = 0.01; adapt_t0 = 15.; adapt_κ = 0.75
                pre = log(2.0) - sqrt(1.0 + t′)/(adapt_γ*(1.0 + t′ + adapt_t0))*log((1+acc[i])/(1.0 + 0.3*t′))
                η = (1 + t′)^(-adapt_κ)
                F.σ[i] = exp(η*pre + (1-η)*log(F.σ[i]))
                θ[i] = F.σ[i]*sign(θ[i])
            else
                if adaptscale
                    effi = (1 + 2*F.ρ/(1 - F.ρ))
                    τ = effi/(t[i]*F.λref)
                    if τ < 0.2
                        F.σ[i] = F.σ[i]*exp(((0.3t[i]/acc[i] > 1.66) - (0.3t[i]/acc[i] < 0.6))*0.03*min(1.0, sqrt(τ/F.λref)))
                    end
                end
                if F isa ZigZag && eltype(θ) <: Number
                    θ[i] = F.σ[i]*rand(rng, (-1,1))
                else
                    θ[i] = F.ρ*θ[i] + F.ρ̄*F.σ[i]*randn(rng, eltype(θ))
                end
            end

            #renew refreshment
            Q[(n + 1)] = t′ + waiting_time_ref(F)
            #update reflections
            for j in neighbours(G1, i)
                b[j] = ab(G1, j, x, θ, c, F)
                t_old[j] = t[j]
                Q[j] = t[j] + poisson_time(b[j], rand(rng))
            end
        else
            # G neighbours have moved (and all others are not our responsibility)

            ∇ϕi = ∇ϕ_(∇ϕ, t, x, θ, i, t′, F, args...)
            l, lb = sλ(∇ϕi, i, x, θ, F), sλ̄(b[i], t[i] - t_old[i])
            num += 1
            if rand(rng)*lb < l
                acc[i] += 1
                if l >= lb
                    !adapt && error("Tuning parameter `c` too small.")
                    #acc .= 0
                    #num = 0
                    adapt!(c, i, factor)
                end
                t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
                θ = reflect!(i, ∇ϕi, x, θ, F)
                for j in neighbours(G1, i)
                    b[j] = ab(G1, j, x, θ, c, F)
                    t_old[j] = t[j]
                    Q[j] = t[j] + poisson_time(b[j], rand(rng))
                end
            else
                b[i] = ab(G1, i, x, θ, c, F)
                t_old[i] = t[i]
                Q[i] = t[i] + poisson_time(b[i], rand(rng))
                continue
            end
        end
        return event(i, t, x, θ, F), t, x, θ, t′, (acc, num), c, b, t_old
    end
end

"""
    spdmp(∇ϕ, t0, x0, θ0, T, c, [G,] F::Union{ZigZag,FactBoomerang}, args...;
        factor=1.5, adapt=false)
        = Ξ, (t, x, θ), (acc, num), c

Version of spdmp which assumes that `i` only depends on coordinates
`x[j] for j in neighbours(G, i)`.

It returns a `FactTrace` (see [`Trace`](@ref)) object `Ξ`, which can be collected
into pairs `t => x` of times and locations and discretized with `discretize`.
Also returns the `num`ber of total and `acc`epted Poisson events and updated bounds
`c` (in case of `adapt==true` the bounds are multiplied by `factor` if they turn
out to be too small.) The final time, location and momentum at `T` can be obtained
with `smove_forward!(t, x, θ, T, F)`.
"""
function spdmp(∇ϕ, t0, x0, θ0, T, c, G, F::Union{ZigZag,FactBoomerang,JointFlow}, args...;
        factor=1.8, adapt=false, adaptscale=false, progress=false, progress_stops = 20, seed=Seed())
    n = length(x0)

    rng = Rng(seed)
    t′ = t0
    t = fill(t′, size(θ0)...)
    t_old = copy(t)
    G1 = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]
    if G == Matched()
        G = G1
    end
    if G == All()
        G2 = nothing
    else
        @assert all(a.second ⊇ b.second for (a,b) in zip(G, G1))
        G2 = [i => setdiff(union((G1[j].second for j in G1[i].second)...), G[i].second) for i in eachindex(G1)]
    end
    x, θ = copy(x0), copy(θ0)
    num = 0
    acc = zeros(Int, length(θ))
    Q = SPriorityQueue{Int,Float64}()
    b = [ab(G1, i, x, θ, c, F) for i in eachindex(θ)]
    for i in eachindex(θ)
        enqueue!(Q, i => poisson_time(b[i], rand(rng)))
    end
    if hasrefresh(F)
        enqueue!(Q, (n + 1) => waiting_time_ref(rng, F))
    end
    if progress
        prg = Progress(progress_stops, 1)
    else
        prg = missing
    end
    stops = ismissing(prg) ? 0 : max(prg.n - 1, 0) # allow one stop for cleanup
    tstop = T/stops
    Ξ = Trace(t0, x0, θ0, F)
    while t′ < T
        ev, t, x, θ, t′, (acc, num), c,  b, t_old = spdmp_inner!(rng, G, G1, G2, ∇ϕ, t, x, θ, Q,
                    c, b, t_old, (acc, num), F, args...; factor=factor, adapt=adapt, adaptscale=adaptscale)
        push!(Ξ, ev)
        if t′ > tstop
            tstop += T/stops
            next!(prg)
        end
   #     @show t′, T
    end
    ismissing(prg) || ProgressMeter.finish!(prg)
    #t, x, θ = smove_forward!(t, x, θ, T, F)
    Ξ, (t, x, θ), (acc, num), c
end

spdmp(∇ϕ, t0, x0, θ0, T, c, F::Union{ZigZag,FactBoomerang,JointFlow}, args...; kargs...) = spdmp(∇ϕ, t0, x0, θ0, T, c, Matched(), F, args...; kargs...)

"""
    pdmp(∇ϕ, t0, x0, θ0, T, c, F::Union{ZigZag, FactBoomerang}, args..., args) = Ξ, (t, x, θ), (acc, num), c

Outer loop of the factorised samplers, the Factorised Boomerang algorithm
and the Zig-Zag sampler. Inputs are a function `∇ϕ` giving `i`th element of gradient
of negative log target density `∇ϕ(x, i, args...)`, starting time and position `t0, x0`,
velocities `θ0`, and tuning vector `c` for rejection bounds and final clock `T`.

The process moves to time `T` with invariant mesure μ(dx) ∝ exp(-ϕ(x))dx and outputs
a collection of reflection points which, together with the initial triple `t`, `x`
`θ` are sufficient for reconstructuing continuously the continuous path.
It returns a `FactTrace` (see [`Trace`](@ref)) object `Ξ`, which can be collected
into pairs `t => x` of times and locations and discretized with `discretize`.
Also returns the `num`ber of total and `acc`epted Poisson events and updated bounds
`c` (in case of `adapt==true` the bounds are multiplied by `factor` if they turn
out to be too small.)

This version does not assume that `∇ϕ` has sparse conditional dependencies,
see [`spdmp`](@ref).
"""
pdmp(∇ϕ, t0, x0, θ0, T, c, F::Union{ZigZag,FactBoomerang,JointFlow}, args...; kargs...) = spdmp(∇ϕ, t0, x0, θ0, T, c, All(), F, args...; kargs...)


text:
using SparseArrays, PyPlot

function orig_xy(xbar,ybar, A, B, H)
  y = ybar * H
  x = xbar * ((1.0-y)*B/2.0+y*(B/2.0+A))
  return  x, y
end

function assemblePoisson(n, f, Apara, B, H)
    h = 1.0 / n #on a grid of 1x1, figure out index spacing
    N = (n+1)^2 #size of all elements, including boundary conditions
    x = h * (0:n) #x coords of the index points
    y = x #equate 2 axis

    umap = reshape(1:N, n+1, n+1)     # Index mapping from 2D grid to vector
    A = Tuple{Int64,Int64,Float64}[]  # Array of matrix elements (row,col,value)
    b = zeros(N)

    # Main loop, insert stencil in matrix for each node point
    for j = 1:n+1
        for i = 1:n+1
            row = umap[i,j]
            if i == 1
                # Dirichlet boundary condition, u = g
                push!(A, (row, row, -1.0))
                push!(A, (row, umap[i+1,j], 1.0))
                b[row] = 0.0
            elseif j == n+1
                push!(A, (row, row, 1.0))
                push!(A, (row, umap[i,j-1], -1.0))
                b[row] = 0.0
            elseif j == 1 || i == n+1
                # Dirichlet boundary condition, u = g
                push!(A, (row, row, 1.0))
                b[row] = 0.0
            else
                # Interior nodes, 5-point stencil
                push!(A, (row, row, 4.0))
                push!(A, (row, umap[i+1,j], -1.0))
                push!(A, (row, umap[i-1,j], -1.0))
                push!(A, (row, umap[i,j+1], -1.0))
                push!(A, (row, umap[i,j-1], -1.0))
                b[row] = 0.5*h^2
            end
        end
    end

    # Create CSC sparse matrix from matrix elements
    A = sparse((x->x[1]).(A), (x->x[2]).(A), (x->x[3]).(A), N, N)

    return A, b, x, y
end

function channelflow(L, B, H, n)
    Apara =  sqrt((L-B)^2/4.0-H^2)
    f(x,y)=0.0
    A, b, x, y = assemblePoisson(n, f, Apara, B, H)

    # Solve + reshape for plotting
    u = reshape(A \ b, n+1, n+1)

    h = 1.0 / n #on a grid of 1x1, figure out index spacing
    Qcap = 0.0
    X=similar(u)
    Y=similar(u)
    for j = 1:n+1
        for i = 1:n+1
           X[i,j],Y[i,j]=orig_xy(x[i],y[j], Apara, B, H)
        end
    end
    for j = 1:n+1
        for i = 1:n+1
           x1,y1=orig_xy(x[i]+h/2,y[j]+h/2, Apara, B, H)
           x0,y0=orig_xy(x[i]-h/2,y[j]-h/2, Apara, B, H)
           Qcap += u[i,j]*abs(x1-x0)*abs(y1-y0)
        end
    end

    # Plotting
    clf()
    contour(X, Y, u, 10, colors="k")
    contourf(X, Y, u, 10)
    axis("equal")
    colorbar()

  return Qcap, x, y, u
end

Q, x, y, u = channelflow(3.0, 0.5, 1.0, 20)


text:
module MarkovBounds

using MultivariateMoments, JuMP, Reexport
@reexport using Catalyst, MomentOpt, SemialgebraicSets, DynamicPolynomials

const MM = MultivariateMoments
const MO = MomentOpt

import Base: show
import LinearAlgebra: qr, nullspace, diag, Diagonal
import Parameters: @unpack

export MarkovProcess, JumpProcess, ReactionProcess, DiffusionProcess,
       JumpDiffusionProcess, ControlProcess,
       Bounds,
       ExitProbability, TerminalSetProbability, LagrangeMayer,
       extended_inf_generator, inf_generator,
       stationary_gmp, transient_gmp,
       stationary_mean, transient_mean,
       stationary_variance, transient_variance,
       stationary_probability_mass,
       optimal_control

include("processes.jl")
include("gmp.jl")
include("utils.jl")
include("distributed.jl")
end


text:
using ReadCoverage
using Test

@testset "percentile_list" begin
    x = [101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450]
    y = [103, 105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 125, 127, 129, 131, 133, 135, 137, 139, 141, 143, 145, 147, 149, 188, 203, 205, 207, 209, 211, 213, 215, 217, 219, 221, 223, 225, 227, 229, 231, 233, 235, 237, 239, 241, 243, 245, 247, 249, 276, 302, 304, 306, 308, 310, 312, 314, 316, 318, 320, 322, 324, 326, 328, 330, 332, 334, 336, 338, 340, 342, 344, 346, 348, 363, 402, 404, 406, 408, 410, 412, 414, 416, 418, 420, 422, 424, 426, 428, 430, 432, 434, 436, 438, 440, 442, 444, 446, 448, 450]

    @test ReadCoverage.percentile_list(x) == y
    @test collect(1:100) == ReadCoverage.percentile_list(collect(0:100))
end

# path_bam = "examples/RamDA_72h_A09.uniq.q40.bam"
# path_bed12 = "examples/gencode.vM15.primary_assembly.annotation.protein_coding.head.bed"
# out_prefix = "examples/out/RamDA_72h_A09.uniq.q40"
# chrom = "chr19"
# leftpos = 3205000
# rightpos = 3207000
# cov = readcoverage_bam(path_bam, chrom, leftpos, rightpos)
# x = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
# , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
# 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
#  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
# , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
# 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
#  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

# @test x == cov


text:
# # Example 2: Optimization for a perfect entangler

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`perfect_entanglers.ipynb`](@__NBVIEWER_ROOT_URL__/examples/perfect_entanglers.ipynb).
#md #
#md #     Compare this example against the [same example using Krotov's
#md #     method](https://juliaquantumcontrol.github.io/Krotov.jl/stable/examples/perfect_entanglers/).

#md # ``\gdef\op#1{\hat{#1}}``
#md # ``\gdef\init{\text{init}}``
#md # ``\gdef\tgt{\text{tgt}}``

#nb # $
#nb # \newcommand{tr}[0]{\operatorname{tr}}
#nb # \newcommand{diag}[0]{\operatorname{diag}}
#nb # \newcommand{abs}[0]{\operatorname{abs}}
#nb # \newcommand{pop}[0]{\operatorname{pop}}
#nb # \newcommand{aux}[0]{\text{aux}}
#nb # \newcommand{opt}[0]{\text{opt}}
#nb # \newcommand{tgt}[0]{\text{tgt}}
#nb # \newcommand{init}[0]{\text{init}}
#nb # \newcommand{lab}[0]{\text{lab}}
#nb # \newcommand{rwa}[0]{\text{rwa}}
#nb # \newcommand{bra}[1]{\langle#1\vert}
#nb # \newcommand{ket}[1]{\vert#1\rangle}
#nb # \newcommand{Bra}[1]{\left\langle#1\right\vert}
#nb # \newcommand{Ket}[1]{\left\vert#1\right\rangle}
#nb # \newcommand{Braket}[2]{\left\langle #1\vphantom{#2}\mid{#2}\vphantom{#1}\right\rangle}
#nb # \newcommand{op}[1]{\hat{#1}}
#nb # \newcommand{Op}[1]{\hat{#1}}
#nb # \newcommand{dd}[0]{\,\text{d}}
#nb # \newcommand{Liouville}[0]{\mathcal{L}}
#nb # \newcommand{DynMap}[0]{\mathcal{E}}
#nb # \newcommand{identity}[0]{\mathbf{1}}
#nb # \newcommand{Norm}[1]{\lVert#1\rVert}
#nb # \newcommand{Abs}[1]{\left\vert#1\right\vert}
#nb # \newcommand{avg}[1]{\langle#1\rangle}
#nb # \newcommand{Avg}[1]{\left\langle#1\right\rangle}
#nb # \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2}
#nb # \newcommand{Re}[0]{\operatorname{Re}}
#nb # \newcommand{Im}[0]{\operatorname{Im}}
#nb # $

using DrWatson
@quickactivate "GRAPETests"
#jl using Test; println("")

# This example illustrates the optimization towards a perfectly entangling
# two-qubit gate for a system of two transmon qubits with a shared transmission
# line. It uses both the indirect perfect entanglers functional shown in
# Goerz et. al., Phys. Rev. A 91, 062307 (2015) and a direct maximization of
# the gate concurrence and thus demonstrates the optimization for non-analytic
# functions via the calculation of gradients with automatic differentiation.

# ## Hamiltonian and guess pulses

# We will write the Hamiltonian in units of GHz (angular frequency; the factor
# 2π is implicit) and ns:

const GHz = 2π
const MHz = 0.001GHz
const ns = 1.0
const μs = 1000ns;

# The Hamiltonian and parameter are taken from Goerz et. al., Phys. Rev. A 91,
# 062307 (2015)., cf. Table 1 in that Reference.

⊗ = kron
const 𝕚 = 1im
const N = 6  # levels per transmon

using LinearAlgebra
using SparseArrays

function hamiltonian(;
    Ωre,
    Ωim,
    N=N,  # levels per transmon
    ω₁=4.380GHz,
    ω₂=4.614GHz,
    ωd=4.498GHz,
    α₁=-210MHz,
    α₂=-215MHz,
    J=-3MHz,
    λ=1.03,
    use_sparse=:auto
)
    𝟙 = SparseMatrixCSC{ComplexF64,Int64}(sparse(I, N, N))
    b̂₁ = spdiagm(1 => complex.(sqrt.(collect(1:N-1)))) ⊗ 𝟙
    b̂₂ = 𝟙 ⊗ spdiagm(1 => complex.(sqrt.(collect(1:N-1))))
    b̂₁⁺ = sparse(b̂₁')
    b̂₂⁺ = sparse(b̂₂')
    n̂₁ = sparse(b̂₁' * b̂₁)
    n̂₂ = sparse(b̂₂' * b̂₂)
    n̂₁² = sparse(n̂₁ * n̂₁)
    n̂₂² = sparse(n̂₂ * n̂₂)
    b̂₁⁺_b̂₂ = sparse(b̂₁' * b̂₂)
    b̂₁_b̂₂⁺ = sparse(b̂₁ * b̂₂')

    ω̃₁ = ω₁ - ωd
    ω̃₂ = ω₂ - ωd

    Ĥ₀ = sparse(
        (ω̃₁ - α₁ / 2) * n̂₁ +
        (α₁ / 2) * n̂₁² +
        (ω̃₂ - α₂ / 2) * n̂₂ +
        (α₂ / 2) * n̂₂² +
        J * (b̂₁⁺_b̂₂ + b̂₁_b̂₂⁺)
    )

    Ĥ₁re = (1 / 2) * (b̂₁ + b̂₁⁺ + λ * b̂₂ + λ * b̂₂⁺)
    Ĥ₁im = (𝕚 / 2) * (b̂₁⁺ - b̂₁ + λ * b̂₂⁺ - λ * b̂₂)

    if ((N < 5) && (use_sparse ≢ true)) || use_sparse ≡ false
        H = (Array(Ĥ₀), (Array(Ĥ₁re), Ωre), (Array(Ĥ₁im), Ωim))
    else
        H = (Ĥ₀, (Ĥ₁re, Ωre), (Ĥ₁im, Ωim))
    end
    return H

end;

# We choose a pulse duration of 400 ns. The guess pulse amplitude is 35 MHz,
# with a 15 ns switch-on/-off time. The Hamiltonian is written in a rotating
# frame, so in general, the control field is allowed to be complex-valued. We
# separate this into two control fields, one for the real part and one for the
# imaginary part. Initially, the imaginary part is zero, corresponding to a
# field exactly at the frequency of the rotating frame.

using QuantumControl.Shapes: flattop

function guess_pulses(; T=400ns, E₀=35MHz, dt=0.1ns, t_rise=15ns)

    tlist = collect(range(0, T, step=dt))
    Ωre = t -> E₀ * flattop(t, T=T, t_rise=t_rise)
    Ωim = t -> 0.0

    return tlist, Ωre, Ωim

end

tlist, Ωre_guess, Ωim_guess = guess_pulses();

# We can visualize this:

using Plots
Plots.default(
    linewidth               = 3,
    size                    = (550, 300),
    legend                  = :right,
    foreground_color_legend = nothing,
    background_color_legend = RGBA(1, 1, 1, 0.8),
)

function plot_complex_pulse(tlist, Ω; time_unit=:ns, ampl_unit=:MHz, kwargs...)

    ax1 = plot(
        tlist ./ eval(time_unit),
        abs.(Ω) ./ eval(ampl_unit);
        label="|Ω|",
        xlabel="time ($time_unit)",
        ylabel="amplitude ($ampl_unit)",
        kwargs...
    )

    ax2 = plot(
        tlist ./ eval(time_unit),
        angle.(Ω) ./ π;
        label="ϕ(Ω)",
        xlabel="time ($time_unit)",
        ylabel="phase (π)"
    )

    plot(ax1, ax2, layout=(2, 1))

end

plot_complex_pulse(tlist, Ωre_guess.(tlist) + 𝕚 * Ωim_guess.(tlist))

# ## Logical basis for two-qubit gates

# For simplicity, we will be define the qubits in the *bare* basis, i.e.
# ignoring the static coupling $J$.

function ket(i::Int64; N=N)
    Ψ = zeros(ComplexF64, N)
    Ψ[i+1] = 1
    return Ψ
end

function ket(indices::Int64...; N=N)
    Ψ = ket(indices[1]; N=N)
    for i in indices[2:end]
        Ψ = Ψ ⊗ ket(i; N=N)
    end
    return Ψ
end

function ket(label::AbstractString; N=N)
    indices = [parse(Int64, digit) for digit in label]
    return ket(indices...; N=N)
end

basis = [ket("00"), ket("01"), ket("10"), ket("11")];

# ## Defining the optimization problem

# We define the optimization with one objective for each of the four basis
# states:

using QuantumControl

H = hamiltonian(Ωre=Ωre_guess, Ωim=Ωim_guess);

objectives = [MinimalObjective(; initial_state=Ψ, generator=H) for Ψ ∈ basis];

# Note the use of `MinimalObjective` which does not have a `target_state`. This
# is because we will be optimizing for an arbitrary perfect entangler, not for
# a specific quantum gate.

# The optimization is steered by the perfect entanglers distance measure
# $D_{PE}$, that is, the geometric distance of the quantum gate obtained from
# propagating the four basis states to the polyhedron of perfect entanglers in
# the Weyl chamber. Since the logical subspace defining the qubit is embedded
# in the larger Hilbert space of the transmon, there may be loss of population
# from the logical subspace. To counter this possibility in the optimization,
# we add a unitarity measure  to $D_PE$. The two terms are added with equal
# weight.

using QuantumControl.WeylChamber: D_PE, gate_concurrence, unitarity
using QuantumControl.Functionals: gate_functional

J_T_PE = gate_functional(D_PE; unitarity_weight=0.5);

# The `gate_functional` routines used above converts the function `D_PE` that
# receives the gate $Û$ as a 4×4 matrix into a functional of the correct from
# for the `QuantumControl.optimize` routine, which is a function of the
# propagated states.

# We can check that for the guess pulse, we are not implementing a perfect
# entangler

using QuantumControl: propagate_objectives

guess_states = propagate_objectives(objectives, tlist; use_threads=true);

U_guess = [basis[i] ⋅ guess_states[j] for i = 1:4, j = 1:4]

gate_concurrence(U_guess)

#jl @test gate_concurrence(U_guess) < 0.9

# We find that the guess pulse produces a gate in the `W0*` region of the Weyl
# chamber:

using QuantumControl.WeylChamber: weyl_chamber_region
weyl_chamber_region(U_guess)

#jl @test weyl_chamber_region(U_guess) == "W0*"

# That is, the region of the Weyl chamber containing controlled-phase gates with
# a phase $> π$ (Weyl chamber coordinates $c₁ > π/2$, $c₂ < π/4$).

# This in fact allows use to use the perfect entangler functional without
# modification: if the guess pulse were in the "W1" region of the Weyl chamber,
# (close to SWAP), we would have to flip its sign, or we would optimize towards
# the local equivalence class of the SWAP gate instead of towards the perfect
# of perfect entanglers. In principle, we could use a modified functional that
# takes the absolute square of the `D_PE` term, by using
#
# ```
# J_T_PE = gate_functional(D_PE; unitarity_weight=0.5, absolute_square=true)
# ```
#
# This would specifically optimize for the *surface* of the perfect
# entanglers functional.

# The guess pulse loses about 10% of population from the logical subspace:

1 - unitarity(U_guess)

#jl @test round(1 - unitarity(U_guess), digits=1) ≈ 0.1

# We can also evaluate the geometric distance to the polyhedron of perfect
# entanglers in the Weyl chamber:

D_PE(U_guess)

# Together with the unitarity measure, this is the initial value of the
# optimization functional:

0.5 * D_PE(U_guess) + 0.5 * (1 - unitarity(U_guess))
#-
J_T_PE(guess_states, objectives)

#jl @test 0.4 < J_T_PE(guess_states, objectives) < 0.5
#jl @test 0.5 * D_PE(U_guess) + 0.5 * (1-unitarity(U_guess)) ≈ J_T_PE(guess_states, objectives) atol=1e-15

# ## Optimization

# Now, we formulate the full control problem

problem = ControlProblem(
    objectives=objectives,
    pulse_options=IdDict(
        Ωre_guess => Dict(
            :lambda_a => 10.0,
            :update_shape => t -> flattop(t, T=400ns, t_rise=15ns, func=:blackman),
        ),
        Ωim_guess => Dict(
            :lambda_a => 10.0,
            :update_shape => t -> flattop(t, T=400ns, t_rise=15ns, func=:blackman),
        ),
    ),
    tlist=tlist,
    iter_stop=100,
    J_T=J_T_PE,
    check_convergence=res -> begin
        (
            (res.J_T > res.J_T_prev) &&
            (res.converged = true) &&
            (res.message = "Loss of monotonic convergence")
        )
        (
            (res.J_T <= 1e-3) &&
            (res.converged = true) &&
            (res.message = "Found a perfect entangler")
        )
    end,
    use_threads=true,
);

# Note that we have not specified a routine to calculate the gradient.
# In  this case, our implementation of GRAPE will use automatic
# differentiation, or actually "semi-automatic differentiation",
# to calculate the gradient directly from the functional.
#
# By semi-automatic differentiation, we mean that when we consider the
# functional as a function of a set of states ``\{|ϕ_k(T)⟩\}``, we find
#
# ```math
# (∇J_T)_{ln} ≡ \frac{∂J_T(\{|ϕ_k(T)⟩\})}{∂ϵ_{ln}}
#  = 2\Re\sum_k
#     \frac{∂J_T}{∂|ϕ_k(T)⟩}
#     \frac{∂|ϕ_k(T)⟩}{∂ϵ_{ln}}
# ```
#
# We can see that ``\frac{∂J_T}{∂|ϕ_k(T)⟩}`` is very closely related to
# the definition ``|χ_k⟩ = -∂J_T/∂⟨ϕ_k|`` (up to the sign) that is a core part
# of [Krotov's method](https://github.com/JuliaQuantumControl/Krotov.jl) (an
# alternative optimization method to GRAPE). Our GRAPE
# implementation knows about this, too. It will detect whether `J_T` is a
# direct function of the forward-propagated states ``|ϕ_k(T)⟩``. It does this
# by checking whether the objectives define a target state. If they do, this
# indicates that `J_T` is a direct function of overlaps ``τ_k``. If they do not
# (our case, see the use `MinimalObjective` above), `J_T` is considered a
# direct function of ``\{|ϕ_k(T)⟩\}``. GRAPE will then look for a
# parameter `chi` that implements the above definition, or use automatic
# differentiation to generate one, which is what happens in this example.
#
# Once ``\frac{∂J_T}{∂|ϕ_k(T)⟩}`` has been calculated, GRAPE is able to
# analytically determine the full gradient ``(∇J_T)_{ln}`` without any further
# use of automatic differentiation.

# With this, we can easily find a solution to the control problem:

opt_result, file =
    @optimize_or_load(datadir(), problem; method=:GRAPE, filename="PE_OCT.jld2");
#-
opt_result

# ## Optimization result

# We extract the optimized control field from the optimization result and plot
# it

Ω_opt = opt_result.optimized_controls[1] + 𝕚 * opt_result.optimized_controls[2]

plot_complex_pulse(tlist, Ω_opt)

# We then propagate the optimized control field to analyze the resulting
# quantum gate:

opt_states = propagate_objectives(
    objectives,
    tlist;
    use_threads=true,
    controls_map=IdDict(
        Ωre_guess => opt_result.optimized_controls[1],
        Ωim_guess => opt_result.optimized_controls[2]
    )
);

U_opt = [basis[i] ⋅ opt_states[j] for i = 1:4, j = 1:4];

# We find that we have achieved a perfect entangler:

gate_concurrence(U_opt)
#jl @test round(gate_concurrence(U_opt), digits=3) ≈ 1.0

# Moreover, we have reduced the population loss to less than 2%

1 - unitarity(U_opt)

#jl @test 1 - unitarity(U_opt) < 0.02

# ## Direct maximization of the gate concurrence

# In the previous optimizations, we have optimized for a perfect entangler
# indirectly via a geometric function in the Weyl chamber. The entire reason
# that perfect entangler functional was formulated is because calculating the
# gate concurrence directly involves the eigenvalues of the unitary, see
# Kraus, Cirac, Phys. Rev. A 63, 062309 (2001) and
# Childs et al., PRA 68, 052311 (2003), which are inherently non-analytic.

# However, since we are able to obtain gradient from automatic differentiation,
# this is no longer an insurmountable obstacle

# We can define a functional for a given gate `U` that combines the gate
# concurrence and (as above) a unitarity measure to penalize loss of population
# from the logical subspace:

J_T_C = U -> 0.5 * (1 - gate_concurrence(U)) + 0.5 * (1 - unitarity(U));

# In the optimization, we will convert this functional to one that takes the
# propagated states as arguments (via the `gate_functional` routine).
# We can help the automatic differentiation making the same conversion "by
# hand". As discussed above, without further intervention, out GRAPE
# implementation will use automatic differentiation to calculate a state
# ``|χ_k⟩ = -\frac{∂}{∂⟨ϕ_k|} J_T``. Here, we can treat `J_T` as depending on
# the gate `U` and apply the chain rule:
#
# ```math
# \begin{split}
#     |χ_k⟩
#     &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} \frac{∂ U_{ik}}{∂⟨ϕ_k|} \\
#     &= - \frac{1}{2} \sum_i (∇_U J_T)_{ik} |Ψ_i⟩
# \end{split}
# ```
#
# The automatic differentiation now only has to calculate ``(∇_U J_T)_{ik}``,
# which is a much smaller-dimensional problem. The above chain rule is
# implemented like this:

using QuantumControl.Functionals: make_gate_chi

chi_C = make_gate_chi(J_T_C, objectives);

# which we can pass explicitly to the optimization.

# Running the optimization, we again are able to find a perfect entangler.

opt_result_direct, file = @optimize_or_load(
    datadir(),
    problem;
    method=:GRAPE,
    J_T=gate_functional(J_T_C),
    chi=chi_C,
    filename="PE_OCT_direct.jld2"
);
#-
opt_result_direct
#-
opt_states_direct = propagate_objectives(
    objectives,
    tlist;
    use_threads=true,
    controls_map=IdDict(
        Ωre_guess => opt_result_direct.optimized_controls[1],
        Ωim_guess => opt_result_direct.optimized_controls[2]
    )
);

U_opt_direct = [basis[i] ⋅ opt_states_direct[j] for i = 1:4, j = 1:4];
#-
gate_concurrence(U_opt_direct)
#jl @test round(gate_concurrence(U_opt_direct), digits=3) ≈ 1.0
#-
1 - unitarity(U_opt_direct)
#jl @test round(1 - unitarity(U_opt_direct), digits=3) ≈ 0.002


text:
# Julia wrapper for header: /opt/ffmpeg/include/libavutil/mem.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_malloc,
    av_malloc_array,
    av_realloc,
    av_realloc_f,
    av_reallocp,
    av_realloc_array,
    av_reallocp_array,
    av_free,
    av_mallocz,
    av_calloc,
    av_mallocz_array,
    av_strdup,
    av_memdup,
    av_freep,
    av_dynarray_add,
    av_dynarray_add_nofree,
    av_dynarray2_add,
    av_size_mult,
    av_max_alloc,
    av_memcpy_backptr,
    av_fast_realloc,
    av_fast_malloc

function av_malloc(size::Csize_t)
    ccall((:av_malloc,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_malloc_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_malloc_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_realloc(ptr,size::Csize_t)
    ccall((:av_realloc,libavutil),Ptr{Void},(Ptr{Void},Csize_t),ptr,size)
end

function av_realloc_f(ptr,nelem::Csize_t,elsize::Csize_t)
    ccall((:av_realloc_f,libavutil),Ptr{Void},(Ptr{Void},Csize_t,Csize_t),ptr,nelem,elsize)
end

function av_reallocp(ptr,size::Csize_t)
    ccall((:av_reallocp,libavutil),Cint,(Ptr{Void},Csize_t),ptr,size)
end

function av_realloc_array(ptr,nmemb::Csize_t,size::Csize_t)
    ccall((:av_realloc_array,libavutil),Ptr{Void},(Ptr{Void},Csize_t,Csize_t),ptr,nmemb,size)
end

function av_reallocp_array(ptr,nmemb::Csize_t,size::Csize_t)
    ccall((:av_reallocp_array,libavutil),Cint,(Ptr{Void},Csize_t,Csize_t),ptr,nmemb,size)
end

function av_free(ptr)
    ccall((:av_free,libavutil),Void,(Ptr{Void},),ptr)
end

function av_mallocz(size::Csize_t)
    ccall((:av_mallocz,libavutil),Ptr{Void},(Csize_t,),size)
end

function av_calloc(nmemb::Csize_t,size::Csize_t)
    ccall((:av_calloc,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_mallocz_array(nmemb::Csize_t,size::Csize_t)
    ccall((:av_mallocz_array,libavutil),Ptr{Void},(Csize_t,Csize_t),nmemb,size)
end

function av_strdup(s)
    ccall((:av_strdup,libavutil),Ptr{Uint8},(Ptr{Uint8},),s)
end

function av_memdup(p,size::Csize_t)
    ccall((:av_memdup,libavutil),Ptr{Void},(Ptr{Void},Csize_t),p,size)
end

function av_freep(ptr)
    ccall((:av_freep,libavutil),Void,(Ptr{Void},),ptr)
end

function av_dynarray_add(tab_ptr,nb_ptr,elem)
    ccall((:av_dynarray_add,libavutil),Void,(Ptr{Void},Ptr{Cint},Ptr{Void}),tab_ptr,nb_ptr,elem)
end

function av_dynarray_add_nofree(tab_ptr,nb_ptr,elem)
    ccall((:av_dynarray_add_nofree,libavutil),Cint,(Ptr{Void},Ptr{Cint},Ptr{Void}),tab_ptr,nb_ptr,elem)
end

function av_dynarray2_add(tab_ptr,nb_ptr,elem_size::Csize_t,elem_data)
    ccall((:av_dynarray2_add,libavutil),Ptr{Void},(Ptr{Ptr{Void}},Ptr{Cint},Csize_t,Ptr{Uint8}),tab_ptr,nb_ptr,elem_size,elem_data)
end

function av_size_mult(a::Csize_t,b::Csize_t,r)
    ccall((:av_size_mult,libavutil),Cint,(Csize_t,Csize_t,Ptr{Csize_t}),a,b,r)
end

function av_max_alloc(max::Csize_t)
    ccall((:av_max_alloc,libavutil),Void,(Csize_t,),max)
end

function av_memcpy_backptr(dst,back::Integer,cnt::Integer)
    ccall((:av_memcpy_backptr,libavutil),Void,(Ptr{Uint8},Cint,Cint),dst,back,cnt)
end

function av_fast_realloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_realloc,libavutil),Ptr{Void},(Ptr{Void},Ptr{Uint32},Csize_t),ptr,size,min_size)
end

function av_fast_malloc(ptr,size,min_size::Csize_t)
    ccall((:av_fast_malloc,libavutil),Void,(Ptr{Void},Ptr{Uint32},Csize_t),ptr,size,min_size)
end


text:
using Test
using EliminateGraphs

@testset "graphlib" begin
    # random graph
    eg = rand_eg(10, 0.5)
    @test check_validity(eg)
    @test eg isa EliminateGraph
    @test nv(eg) == nv0(eg) == 10

    # empty and full graph
    g = empty_eg(10)
    @test check_validity(eg)
    @test nv(g) == 10
    @test ne(g) == 0
    g = K_eg(4)
    @test check_validity(eg)
    @test ne(g) == 6
    g = K_eg(3,3)
    @test check_validity(eg)
    @test nv(g) == 6
    @test ne(g) == 9
    g =  disconnected_cliques_eg(3,3)
    @test check_validity(eg)
    @test nv(g) == 6
    @test ne(g) == 6

    # chain and ring
    g = chain_eg(10)
    @test check_validity(eg)
    @test nv(g) == 10
    @test ne(g) == 9
    g = ring_eg(10)
    @test check_validity(eg)
    @test nv(g) == 10
    @test ne(g) == 10

    # constant graphs
    @test check_validity(petersen_graph)
    @test petersen_graph |> ne == 15
    @test petersen_graph |> nv == 10
    @test check_validity(c60_graph)
    @test c60_graph |> ne == 90
    @test c60_graph |> nv == 60
end


text:
# This file was generated, do not modify it. # hide
train, test = partition(collect(eachindex(y_wind)), 0.7, shuffle=true, rng=5);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function build_reaction_table_from_vff_file(path_to_vff_file::String)::Some

    try

        # initialize -
        kegg_reaction_array = Array{KEGGReaction,1}()

        # get file buffer (array of strings) -
        vff_file_buffer = read_file_from_path(path_to_vff_file)

        # process -
        for reaction_line in vff_file_buffer

            # build KEGG reaction -
            reaction_object = KEGGReaction()

            # populate the data fields, store the reaction objects in an array

            # skip comments and empty lines -
            if (occursin("//", reaction_line) == false &&
                isempty(reaction_line) == false)

                # split -
                reaction_components = string.(split(reaction_line, ","))

                # ecnumber -
                ec_number = reaction_components[1]
                if (isnothing(ec_number) == true)
                    reaction_object.ec_number = missing
                else
                    reaction_object.ec_number = ec_number
                end

                # @show reaction_components

                # reaction number -
                rn_number = reaction_components[2]
                isnothing(rn_number) ? reaction_object.reaction_number = missing : reaction_object.reaction_number = rn_number

                # enzyme_name -
                enzyme_name = reaction_components[3]
                isnothing(enzyme_name) ? reaction_object.enzyme_name = missing : reaction_object.enzyme_name = enzyme_name

                # reaction_markup -
                reaction_markup = reaction_components[4]
                isnothing(reaction_markup) ? reaction_object.reaction_markup = missing : reaction_object.reaction_markup = reaction_markup

                # skip the rest if we have no reaction ...
                if (isnothing(reaction_markup) == false)

                    # ok, for the next fields, we need to do a little parsing ... split around <=>
                    reaction_phrases = string.(split(reaction_markup, "<=>") .|> lstrip .|> rstrip)

                    # reactant phrase -
                    left_phrase = reaction_phrases[1]
                    isnothing(left_phrase) ? reaction_object.reaction_forward = missing : reaction_object.reaction_forward = left_phrase

                    # product phrase -
                    right_phrase = reaction_phrases[2]
                    isnothing(right_phrase) ? reaction_object.reaction_reverse = missing : reaction_object.reaction_reverse = right_phrase

                    # build the stoichiometric_dictionary -
                    stoichiometric_dictionary = extract_stoichiometric_dictionary(reaction_markup)
                    reaction_object.stoichiometric_dictionary = stoichiometric_dictionary

                    # grab -
                    push!(kegg_reaction_array, reaction_object)
                end
            end
        end

        # return -
        return Some(kegg_reaction_array)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end


text:
# TODO There is https://github.com/cortner/SphericalHarmonics.jl,
#      which should be used, once it's ready.

"""
Returns the (l,m) real spherical harmonic Y_lm(r). Consistent with
https://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics
"""
function ylm_real(l::Integer, m::Integer, rvec::AbstractVector{T}) where T
    @assert 0 ≤ l
    @assert -l ≤ m ≤ l
    @assert length(rvec) == 3
    x, y, z = rvec
    r = norm(rvec)

    if l == 0  # s
        (m ==  0) && return sqrt(1 / 4T(π))
    end

    # Catch cases of numerically very small r
    if r <= 10 * eps(eltype(rvec))
        return zero(T)
    end

    if l == 1  # p
        (m == -1) && return sqrt(3 / 4T(π)) * y / r
        (m ==  0) && return sqrt(3 / 4T(π)) * z / r
        (m ==  1) && return sqrt(3 / 4T(π)) * x / r
    end

    if l == 2  # d
        (m == -2) && return sqrt(15 /  4T(π)) * (x / r) * (y / r)
        (m == -1) && return sqrt(15 /  4T(π)) * (y / r) * (z / r)
        (m ==  0) && return sqrt( 5 / 16T(π)) * (2z^2 - x^2 - y^2) / r^2
        (m ==  1) && return sqrt(15 /  4T(π)) * (x / r) * (z / r)
        (m ==  2) && return sqrt(15 /  4T(π)) * (x^2 - y^2) / r^2
    end

    if l == 3  # f
        (m == -3) && return sqrt( 35 / 32T(π)) * (3x^2 - y^2) * y / r^3
        (m == -2) && return sqrt(105 /  4T(π)) * x * y * z / r^3
        (m == -1) && return sqrt( 21 / 32T(π)) * y * (4z^2 - x^2 - y^2) / r^3
        (m ==  0) && return sqrt(  7 / 16T(π)) * z * (2z^2 - 3x^2 - 3y^2) / r^3
        (m ==  1) && return sqrt( 21 / 32T(π)) * x * (4z^2 - x^2 - y^2) / r^3
        (m ==  2) && return sqrt(105 / 16T(π)) * (x^2 - y^2) * z / r^3
        (m ==  3) && return sqrt( 35 / 32T(π)) * (x^2 - 3y^2) * x / r^3
    end

    error("The case l = $l and m = $m is not implemented")
end


text:
using DataFrames
using Levenshtein
using BioSymbols
using BioSequences
using IterTools

using Rifraf

"""indices of reads with right length and error cutoff"""
function valid_read_indices(error_range, sequences, phreds, mean_errors; length_cutoff=LENGTH_CUTOFF)
    (min_error, max_error) = error_range
    median_length = median(map(length, sequences))
    return [i for i in 1:length(phreds) if ((min_error <= mean_errors[i] <= max_error) &&
                                            (abs(length(sequences[i]) - median_length) < length_cutoff))]
end

"""random sample of `nseqs` reads"""
function sample_reads(nseqs::Int,
                      cand_indices::Vector{Int},
                      sequences::Vector{DNASeq},
                      phreds::Vector{Vector{Rifraf.Phred}},
                      names::Vector{String},
                      mean_errors::Vector{Float64})
    indices = rand(cand_indices, nseqs)
    sampled_seqs = DNASeq[sequences[i] for i in indices]
    sampled_phreds = Vector{Rifraf.Phred}[phreds[i] for i in indices]
    sampled_names = String[names[i] for i in indices]
    if any([length(s) != length(p) for (s, p) in zip(sampled_seqs, sampled_phreds)])
        for (ss, sp) in zip(sampled_seqs, sampled_phreds)
            println("$(length(ss)), $(length(sp))")
            println("")
        end
    end
    return sampled_seqs, sampled_phreds, sampled_names, indices
end

"""alignment moves"""
function get_moves(seq, reference, align_errors)
    align_scores = Scores(align_errors)
    align_bandwidth = 100
    pseq = Rifraf.RifrafSequence(DNASeq(seq), fill(-1.0, length(seq)),
                                 align_bandwidth, align_scores)
    moves = Rifraf.align_moves(DNASeq(reference), pseq; trim=true)
    moves
end

"""get [start, stop] indices to trim read ends by aligning to reference"""
function trim_ends_indices(seq, reference)
    align_errors = ErrorModel(1e5, 1e-3, 1e-3, 0.0, 0.0)
    moves = get_moves(seq, reference, align_errors)
    x = 1
    while moves[x] == Rifraf.TRACE_INSERT
        x += 1
    end
    y = length(moves)
    while moves[y] == Rifraf.TRACE_INSERT
        y -= 1
    end
    n_end = length(moves) - y
    y = length(seq) - n_end
    return x, y
end

"""align to reference and trim, to remove primers"""
function trim_reads(seqs, phreds, reference)
    refstring = String(reference)
    final_seqs = DNASeq[]
    final_phreds = Vector{Rifraf.Phred}[]
    for (seq, phred) in zip(seqs, phreds)
        if length(seq) != length(phred)
            error("mismatched lengths")
        end
        seqstring = String(seq)

        # take reverse complement, if necessary
        d1 = Levenshtein.levenshtein(seqstring, refstring)
        d2 = Levenshtein.levenshtein(String(reverse_complement(seq)),
                                     refstring)
        if d1 > d2
            seq = reverse_complement(seq)
            phred = reverse(phred)
        end

        # trim
        x, y = trim_ends_indices(seq, reference)
        push!(final_seqs, seq[x:y])
        push!(final_phreds, phred[x:y])
    end
    return final_seqs, final_phreds
end

function shorten_seq_index(seq, reference, target_len; codon_moves=false)
    align_errors = if codon_moves
        ErrorModel(10, 1e-5, 1e-5, 1.0, 1.0)
    else
        ErrorModel(10., 1., 1., 0.0, 0.0)
    end
    moves = get_moves(seq, reference, align_errors)
    indices = Rifraf.moves_to_indices(moves, length(reference), length(seq))
    stop = indices[target_len]
    if stop % 3 > 0
        stop += (3 - (stop % 3))
        stop = min(length(reference), stop)
    end
    return stop
end

function insert_base_helper(template, seq, phreds, index, base, align_errors)
    moves = get_moves(seq, template, align_errors)
    indices = Rifraf.moves_to_indices(moves, length(template), length(seq))
    read_index = indices[index]

    error_rate = mean(Rifraf.phred_to_p(phreds))
    final_seq, final_phreds = if rand(Bernoulli(error_rate)) == 1
        # small chance of sequencing error causing a deletion in the read
        # chance of missing base depends on sequence error rate
        seq, phreds
    else
        final_seq = DNASeq(seq[1:(read_index-1)], DNASeq([base]), seq[read_index:end])
        # uniform between neighbor phreds.
        phred1 = phreds[max(read_index-1, 1)]
        phred2 = phreds[min(read_index, length(template))]
        new_phred = rand(min(phred1,phred2):max(phred1,phred2))
        final_phreds = vcat(phreds[1:(read_index-1)], [new_phred], phreds[read_index:end])
        final_seq, final_phreds
    end
    return final_seq, final_phreds
end

function delete_base_helper(template, seq, phreds, index, align_errors)
    moves = get_moves(seq, template, align_errors)
    indices = Rifraf.moves_to_indices(moves, length(template), length(seq))
    read_index = indices[index]
    final_seq = copy(seq)
    final_phreds = copy(phreds)
    deleteat!(final_seq, read_index)
    deleteat!(final_phreds, read_index)
    return final_seq, final_phreds
end

"""

Get a random position in `sequence`. If `in_hompolymer`, it must be in
a homopolymer region at least `length` bp long. Otherwise it must be
in a non-homopolymer region.

"""
function sample_position(sequence, in_homopolymer; hplength=4)
    # run-length encoding
    groups = IterTools.groupby(x -> x, sequence)
    lengths = map(length, groups)

    # the index of each group in the overall sequence
    indices = [i + 1 for i in cumsum(lengths)[1:(end-1)]]
    insert!(indices, 1, 1)

    pairs = zip(lengths, indices)

    pairs = if in_homopolymer
        # filter out illegal positions
        [p for p in pairs if p[1] >= hplength]
    else
        [p for p in pairs if p[1] < hplength]
    end
    len, start = rand(pairs)
    return start + rand(1:len) - 1
end

function insert_base(template, seqs, phreds, in_homopolymer)
    # sample an existing base
    base_index = sample_position(template, in_homopolymer)
    # equal probability of insertion on either side of that base
    index = rand([base_index, base_index + 1])
    base = if in_homopolymer
        # match homopolymer bases
        template[base_index]
    else
        # do not match neighbors
        rand([DNA_A, DNA_C, DNA_G, DNA_T])
    end

    final_template = DNASeq(DNASeq(template[1:(index-1)]), DNASeq([base]), DNASeq(template[index:end]))
    align_errors = ErrorModel(10., 1., 1., 0.0, 0.0)
    final_seqs = DNASeq[]
    final_phreds = Vector{Rifraf.Phred}[]
    for (seq, phred) in zip(seqs, phreds)
        seq, phred = insert_base_helper(template, seq, phred, index, base, align_errors)
        push!(final_seqs, seq)
        push!(final_phreds, phred)
    end
    return (final_template, final_seqs, final_phreds, index)
end

function delete_base(template, seqs, phreds, in_homopolymer)
    index = sample_position(template, in_homopolymer)
    final_template = DNASeq(DNASeq(template[1:(index-1)]), DNASeq(template[(index+1):end]))
    align_errors = ErrorModel(10., 1., 1., 0.0, 0.0)

    final_seqs = DNASeq[]
    final_phreds = Vector{Rifraf.Phred}[]
    for (seq, phred) in zip(seqs, phreds)
        seq, phred = delete_base_helper(template, seq, phred, index, align_errors)
        push!(final_seqs, seq)
        push!(final_phreds, phred)
    end
    return (final_template, final_seqs, final_phreds, index)
end


text:
"""
Testing transmission of various types of chromosomes under diploidy.

Under diploidy, there are several kinds of chromosomes such as autosome and sex-chormosomes.
This file currently provides unit tests for Autosome, XChromosome, YChromosome, and
Mitochondrion. The main purpose of types tested here is to handle sex-specific transmission
of genetic materials including recombination. Actual storage of genertic data is delegated
to ChromosomalStorage type.

In order to make holding organisms in a population simpler, chromosomes are not parametrized
by sex. This design incurs wasted storage space in sex-chromosomes and organelle's chromosomes.
"""

facts("Autosome are non sex-specific.") do
    nelems = 10
    mother = Autosome{ChromosomalStorage{Int}}(nelems)
    father = Autosome{ChromosomalStorage{Int}}(nelems)

end

facts("""
All organisms have two autosomes, and the pattern of transmissioin does not depend on sex.
""") do
    nelems = 10

    context("A new chromosome can be initialized without parental chromosomes.") do
        Autosome(chr1, chr2)

    end
    context("A new chromosome can be created with one maternal and one paternal chromosomes.") do

    end
end

facts("""
    A female contains two X-chromosome. A female offspring inherits each of X-chromosomes from
    its mother and father.
    """) do
    context("A new chromosome can be initialized without parental chromosomes.") do
    end

    context("A new chromosome can be created from one maternal and one paternal chromosomes.") do
    end
end

facts("A male only contains maternally inherited X-chromosome.") do
    context("A new chromosome can be initialized without parental chromosomes.") do
    end

    context("""
        A new chromosome can be created from one maternal and one paternal chromosomes.
        However, the paternal chromosome will be ignored.
        """) do
    end
end

facts("A female does not contain Y-chromosome.") do
    context("""
        A new chromosome can be initialized without parental chromosomes. However,
        the chromosome is intentionally left uninitialized.
    """) do
    end

    context("""
        A new chromosomes can be created from one maternal and one paternal chromosomes.
        However, both chromosomes will be ignored.
    """) do
    end
end

facts("A male only contains paternally inheritec Y-chromosome.") do
    context("A new chromosome can be initialized without parental chromosomes.") do
    end

    context("""
        A new chromosome can be created from one maternal and one paternal chromosomes.
        However, the maternal chromosome will be ignored.
    """) do
    end
end

facts("A female contains maternally inherited micohondrion.") do
    context("A new chromosome can be initialized without parental chromosomes.") do
    end

    context("""
        A new chromosome can be created from one maternal and one paternal chromosomes.
        However, the paternal chromosome will be ignored.
    """) do
    end
end

facts("""
    Although male contains maternally inherited mitochondria, a male parent never passes it
    to any offspring.
    """) do
    context("""
        A new chromosome can be initialized without parental chromosomes. However, the
        chromosome is intentionally left uninitialized.
    """) do
    end

    context("""
        A new chromosomes can be created from one maternal and one paternal chromosomes.
        However, both chromosomes will be ignored.
    """) do
    end
end


text:
include("synonym.jl")
include("location.jl")
include("exon.jl")


text:
# using DelimitedFiles
# using Makie
# data = readdlm("../LinneasData190920.csv", ';')
#
# data[isnan.(data)] .= circshift(data, (-1,0))[isnan.(data)]
# data[isnan.(data)] .= circshift(data, (1,0))[isnan.(data)]
# data[isnan.(data)] .= circshift(data, (-2,0))[isnan.(data)]
# data[isnan.(data)] .= circshift(data, (2,0))[isnan.(data)]
# data[isnan.(data)] .= circshift(data, (3,0))[isnan.(data)]
# any(isnan.(data))
#
# #data = replace(data, NaN=>missing)
# #μ = mapslices(mean∘skipmissing, data, dims=1)
# #sigma = mapslices(std∘skipmissing, data, dims=1)
# #surface(0..1, 0..5, data)

# NOTE there is some loss of generality with the way recompute_ODEs is defined
# recompute_ODEs is the same across all samples, in previous implementation
# it could be different, to remedy this, adjust the definition of GibbsDefn in
# workspace.jl
"""
    mcmc(setups::Vector{<:MCMCSetup})

Gibbs sampler alternately imputing unobserved parts of the paths and updating
unknown coordinates of the parameter vector. Version suitable for multiple
trajectory samples. `setups` defines all variables required for the
initialisation of the Markov Chain
"""
function mcmc(setups::Vector{<:MCMCSetup})
    num_mcmc_steps, K = setups[1].num_mcmc_steps, length(setups)
    tu = Workspace(setups[1])
    ws, ll, θ = [tu.workspace], [tu.ll], [tu.θ]
    for k in 2:K
        tu = Workspace(setups[k])
        push!(ws, tu.workspace); push!(ll, tu.ll); push!(θ, tu.θ)
    end
    gibbs = GibbsDefn(setups[1])

    for i in 1:num_mcmc_steps
        verbose = act(Verbose(), ws[1], i)
        act(SavePath(), ws[1], i) && for k in 1:K save_path!(ws[k]) end
        for k in 1:K next_set_of_blocks(ws[k]) end
        for k in 1:K
            ll[k], acc = impute!(ws[k], ll[k], verbose, i)
            update!(ws[k].accpt_tracker, Imputation(), acc)
        end

        if act(ParamUpdate(), ws[1], i)
            for j in 1:length(gibbs)
                ll, acc, θ = update_param!(gibbs[j], θ, ws, ll, verbose, i)

                for k in 1:K
                    update!(ws[k].accpt_tracker, ParamUpdate(), j, acc)
                    update!(ws[k].θ_chain, θ[k])
                end
                verbose && print("\n")
            end
            P˟ = clone(ws[1].P[1].Target, θ[1])
            verbose && println(prod("$v=$x " for (v, x) in zip(param_names(P˟), orig_params(P˟))))
            verbose && print("------------------------------------------------",
                             "------\n")
        end
    end
    ws
end

# NOTE for the Mixed Effect models there will need to be some indicator matrix
# passed arround to indicate which indices are common and how they should be
# sampled. For the repeated.jl all θs that are passed around are the same, so
# for now I discard all the remaining θs [⋆]
function conjugate_draw(θ, XX::Vector{<:Vector}, PT, prior, updtIdx)
    μ = mustart(updtIdx)
    𝓦 = μ*μ'
    ϑ = SVector(thetaex(updtIdx, θ[1])) #NOTE [⋆] hence θ[1]
    for k in 1:length(XX)
        μ, 𝓦 = _conjugate_draw(ϑ, μ, 𝓦, XX[k], PT, updtIdx)
    end
    Σ = inv(𝓦 + inv(Matrix(prior.Σ)))
    Σ = (Σ + Σ')/2 # eliminates numerical inconsistencies
    μ_post = Σ * (μ + Vector(prior.Σ\prior.μ))
    ϑ = rand(Gaussian(μ_post, Σ))
    # expand back to one parameter per path
    [move_to_proper_place(ϑ, θ_k, updtIdx) for θ_k in θ]
end

# NOTE `ws` has UpdtIdx defined only as a type, can change it later to be a
# member so that passing ws[k].updtIdx instead of UpdtIdx() is possible
# NOTE same applies to OS (observation scheme)
#
# no blocking
function update_param!(pu::ParamUpdtDefn{ConjugateUpdt,UpdtIdx}, θ,
                       ws::Vector{<:Workspace{OS,NoBlocking}}, ll,
                       verbose=false, it=NaN) where {UpdtIdx,OS}
    K = length(ws)
    # warn if targets are different?
    θᵒ = conjugate_draw(θ, [w.XX for w in ws], ws[1].P[1].Target, pu.priors[1], UpdtIdx())   # sample new parameter

    for k in 1:K
        WW, Pᵒ, P, XXᵒ, XX = ws[k].WW, ws[k].Pᵒ, ws[k].P, ws[k].XXᵒ, ws[k].XX
        m = length(WW)
        update_laws!(P, θᵒ[k]) # hardcoded: NO Blocking
        pu.recompute_ODEs && solve_back_rec!(NoBlocking(), ws[k], P) # compute (H, Hν, c)

        for i in 1:m    # compute wiener path WW that generates XX
            inv_solve!(EulerMaruyamaBounded(), XX[i], WW[i], P[i])
        end
        # compute white noise that generates starting point
        y = XX[1].yy[1]
        z = inv_start_pt(y, ws[k].x0_prior, P[1])

        llᵒ = logpdf(ws[k].x0_prior, y)
        llᵒ += path_log_likhd(OS(), XX, P, 1:m, ws[k].fpt; skipFPT=true)
        llᵒ += lobslikelihood(P[1], y)
        print_info(verbose, it, value(ll[k]), value(llᵒ))
        verbose && k < K && print("\n")
        ll[k] = llᵒ
        set!(ws[k].z, z)
    end
    return ll, true, θᵒ
end

function update_param!(pu::ParamUpdtDefn{MetropolisHastingsUpdt,UpdtIdx}, θ,
                      ws::Vector{<:Workspace{OS,NoBlocking}}, ll, verbose=false,
                      it=NaN) where {UpdtIdx,OS}
    K = length(ws)
    #NOTE for now, let's sample parameter corresponding to a single path, this
    # will need to be changed for Mixed effect models
    θᵒ = rand(pu.t_kernel, θ[1], UpdtIdx())               # sample new parameter
    llᵒ = copy(ll)
    llr = prior_kernel_contrib(pu.t_kernel, pu.priors, θ[1], θᵒ)
    zᵒ = [copy(w.z.val) for w in ws]
    for k in 1:K
        WW, Pᵒ, P, XXᵒ, XX = ws[k].WW, ws[k].Pᵒ, ws[k].P, ws[k].XXᵒ, ws[k].XX
        m = length(WW)
        update_laws!(Pᵒ, θᵒ)
        pu.recompute_ODEs && solve_back_rec!(NoBlocking(), ws[k], Pᵒ) # compute (H, Hν, c)

        # find white noise which for a given θᵒ gives a correct starting point
        y = XX[1].yy[1]
        zᵒ[k] = inv_start_pt(y, ws[k].x0_prior, Pᵒ[1])

        find_path_from_wiener!(XXᵒ, y, WW, Pᵒ, 1:m)

        llᵒ[k] = logpdf(ws[k].x0_prior, y)
        llᵒ[k] += path_log_likhd(OS(), XXᵒ, Pᵒ, 1:m, ws[k].fpt)
        llᵒ[k] += lobslikelihood(Pᵒ[1], y)

        print_info(verbose, it, ll[k], llᵒ[k])
        verbose && k < K && print("\n")
        llr += llᵒ[k] - ll[k]
    end

    # Accept / reject
    if accept_sample(llr, verbose)
        for k in 1:K
            m = length(ws[k].WW)
            swap!(ws[k].XX, ws[k].XXᵒ, ws[k].P, ws[k].Pᵒ, 1:m)
            set!(ws[k].z, zᵒ[k])
        end
        # expand back to a vector of parameters
        return llᵒ, true, [copy(θᵒ) for _ in θ]
    else
        return ll, false, θ
    end
end


text:
macro _const(ex)
    ex = esc(ex)
    if VERSION < v"1.8.0-DEV.1148"
        return ex
    else
        return Expr(:const, ex)
    end
end

const var"@const" = var"@_const"

const Historic = if USE_HISTORIC
    Base.require(Base.PkgId(Base.UUID(0xe6ec0b50ef98488aa141280ae2eaf113), "Historic"))
else
    nothing
end

@static if Historic === nothing
    macro record(_...)
        nothing
    end
else
    Historic.@define Debug
    using .Debug: @record
end

"""
    @yield_unsafe expression

Document that `expression` must not yield to the Julia scheduler.

This macro is used purely for documentation.  The `expression` is evaluated
as-is.
"""
macro yield_unsafe(ex)
    esc(ex)
end

chaos_takemanyto() = nothing
chaos_run_finalizers() = nothing
queue_reset_done() = nothing
queue_fallback_done() = nothing


text:
# ---
# title: 891. Sum of Subsequence Widths
# id: problem891
# author: Indigo
# date: 2021-02-01

# difficulty: Hard
# categories: Array, Math
# link: <https://leetcode.com/problems/sum-of-subsequence-widths/description/>
# hidden: true
# ---
#
# Given an array of integers `A`, consider all non-empty subsequences of `A`.
#
# For any sequence S, let the  _width_  of S be the difference between the
# maximum and minimum element of S.
#
# Return the sum of the widths of all subsequences of A.
#
# As the answer may be very large, **return the answer modulo 10^9 + 7**.
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,1,3]
#     Output: 6
#     Explanation: Subsequences are [1], [2], [3], [2,1], [2,3], [1,3], [2,1,3].
#     The corresponding widths are 0, 0, 0, 1, 1, 2, 2.
#     The sum of these widths is 6.
#
#
#
#
# **Note:**
#
#   * `1 <= A.length <= 20000`
#   * `1 <= A[i] <= 20000`
#
#
## @lc code=start
using LeetCode

function sum_subseq_widths(A::Vector{Int})
    MOD = Int(1e9 + 7)
    len = length(A)
    sort!(A)
    pow = fill(2, len)
    for i in 2:len
        pow[i] = (pow[i - 1] << 1) % MOD
    end
    res = 0
    for i in 2:(len - 1)
        res += (A[i]) * ((pow[i - 1] - pow[len - i])) % MOD
    end
    return mod1((res + (A[end] - A[1]) * (pow[len - 1] - 1) % MOD), MOD)
end
## @lc code=end


text:
function get_rexhopper(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=1.0,
    contact_foot=true,
    contact_body=true,
    limits=true,
    model=:rexhopper,
    floating=true,
    contact_type=:nonlinear,
    spring=0.0,
    damper=0.0,
    parse_damper=true,
    T=Float64)

    path = joinpath(@__DIR__, "../deps/$(String(model)).urdf")
    mech = Mechanism(path; floating, T,
        gravity,
        timestep,
        parse_damper)

    # Adding springs and dampers
    set_springs!(mech.joints, spring)
    set_dampers!(mech.joints, damper)

    # joint limits
    joints = deepcopy(mech.joints)

    if limits
        joint1 = get_joint(mech, :joint1)
        joint3 = get_joint(mech, :joint3)
        joints[joint1.id] = add_limits(mech, joint1,
            rot_limits=[SVector{1}(-0.7597), SVector{1}(1.8295)])
        joints[joint3.id] = add_limits(mech, joint3,
            rot_limits=[SVector{1}(-1.8295), SVector{1}(0.7597)])
        mech = Mechanism(Origin{T}(), [mech.bodies...], [joints...];
            gravity,
            timestep)
    end

    if contact_foot
        origin = Origin{T}()
        bodies = mech.bodies
        joints = mech.joints

        normal = [0.0; 0.0; 1.0]
        models = []

        link3 = get_body(mech, :link3)
        link2 = get_body(mech, :link2)
        foot_radius = 0.0203
        ankle_radius = 0.025
        base_radius = 0.14
        p = [0.1685; 0.0025; -0.0055]
        o = foot_radius
        push!(models, contact_constraint(link3, normal;
            friction_coefficient,
            contact_origin=p,
            contact_radius=o,
            contact_type,
            name=:foot))
        p = [-0.10; -0.002; 0.01]
        o = ankle_radius
        push!(models, contact_constraint(link3, normal;
            friction_coefficient,
            contact_origin=p,
            contact_radius=o,
            contact_type,
            name=:ankle3))
        p = [0.24; 0.007; 0.005]
        push!(models, contact_constraint(link2, normal;
            friction_coefficient,
            contact_origin=p,
            contact_radius=o,
            contact_type,
            name=:ankle2))
        base_link = get_body(mech, :base_link)
        pl = [0.0; +0.075; 0.03]
        pr = [0.0; -0.075; 0.03]
        o = base_radius
        push!(models, contact_constraint(base_link, normal;
            friction_coefficient,
            contact_origin=pl,
            contact_radius=o,
            contact_type,
            name=:torso_left))
        push!(models, contact_constraint(base_link, normal;
            friction_coefficient,
            contact_origin=pr,
            contact_radius=o,
            contact_type,
            name=:torso_right))

        set_minimal_coordinates!(mech, get_joint(mech, :floating_base), [0,0,1.0, 0,0,0])

        mech = Mechanism(origin, bodies, joints, [models...];
            gravity,
            timestep)
    end
    return mech
end

function initialize_rexhopper!(mechanism::Mechanism{T};
        body_position=zeros(3),
        body_orientation=zeros(3),
        body_linear_velocity=zeros(3),
        body_angular_velocity=zeros(3)) where T

    zero_velocity!(mechanism)
    body_position += [0.0, 0.0, 0.3148]

    for joint in mechanism.joints
        !(joint.name in (:loop_joint, :floating_joint)) && set_minimal_coordinates!(mechanism, joint, zeros(input_dimension(joint)))
    end

    set_minimal_coordinates!(mechanism, get_joint(mechanism,
        :floating_base), [body_position; body_orientation])
    set_minimal_velocities!(mechanism, get_joint(mechanism,
        :floating_base), [body_linear_velocity; body_angular_velocity])
end


text:
import CoverageCore
coveragecore_cov_res = CoverageCore.process_folder()

import Pkg
Pkg.add("Coverage")
import Coverage
coverage_cov_res = Coverage.process_folder()

retry(() -> Coverage.Coveralls.submit(coverage_cov_res); delays = ExponentialBackOff(n = 5))()
retry(() -> Coverage.Codecov.submit(coverage_cov_res); delays = ExponentialBackOff(n = 5))()


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random
path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/ANNBN.jl"))

function ff(x)
    xx=200x .- 100
    n = length(xx)
    1 + (1/4000)*sum(abs2, (xx)) - prod(cos.((xx) ./ sqrt.(1:n)))
end
noise_perc=1;vars=100;
i_train=1000; xx_train=0.8*rand(i_train,vars).+0.1
yy_train=zeros(i_train); for i=1:i_train yy_train[i]=ff(xx_train[i,:]) end
# plot(yy_train[1:100]);yy_train.+=noise_perc.*(1/2 .-rand(i_train));plot!(yy_train[1:100])
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1
i_test=1000;xx_test=0.8*rand(i_test,vars).+0.1
yy_test=zeros(i_test); for i=1:i_test yy_test[i]=ff(xx_test[i,:]) end
# yy_test.+=noise_perc.*rand(i_test)
yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1

##### ANNBN RBF
neurons=Int64(floor(i_train/(vars+1)))
neurons=100 # this automatically generates number of neurons== number of clusters (§3.1)
inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)
cc1=100.0
a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
plot(err[ier],label="ANNBN",color=:black,linestyle=:dashdotdot,legend=:topleft,linewidth=3)
# savefig("Grienwak.pdf")

using DecisionTree, Statistics
num_trees=100;
regr_3 = RandomForestRegressor(n_trees=num_trees);
tmp_xx_train=copy(xx_train);tmp_yy_train=copy(yy_train);
tmp_xx_test=copy(xx_test);tmp_yy_test=copy(yy_test);
DecisionTree.fit!(regr_3, tmp_xx_train, tmp_yy_train[:,1]);
predl1 = DecisionTree.predict(regr_3, tmp_xx_train);
predl1_test = DecisionTree.predict(regr_3, tmp_xx_test);
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
plot!(err[ier],label="Random Forests",color=:black,linestyle=:dash,legend=:topleft,linewidth=3)

using PyCall,ScikitLearn
@sk_import ensemble: AdaBoostRegressor
sci_model = fit!(AdaBoostRegressor(n_estimators=100),xx_train,yy_train)
predl1=predict(sci_model,xx_train)
predl1_test=predict(sci_model,xx_test)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
plot!(err[ier],label="AdaBoost",color=:black,linestyle=:dashdot,legend=:topleft,linewidth=3)

using XGBoost
num_round = 10
bst = xgboost(copy(xx_train), num_round, label = copy(yy_train), eta = 1, max_depth = 7)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
plot!(err[ier],label="XGBoost",color=:black,linestyle=:dot,legend=:topleft,linewidth=3)


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

using Documenter, ImprovedGarbanzo

makedocs(
    sitename="ImprovedGarbanzo",
    modules = [ImprovedGarbanzo],
    pages=[
        "Start" => "index.md"
        "Reference" => "reference.md"
    ],
)

deploydocs(
    repo = "github.com/not7cd/improved-garbanzo.git",
)