text:
"""
    struct CachedPool{O,T}

        Use this pool to host the objects that are heavy to evaluate so that one wants to cache their status.
        The user should defines a compare

# Members
- name::Symbol : name of the pool
- object::O    : object
- current::T      : current status
- new::T       : the new status wants to assign later
- version::Int128 : the current version
- excited::Bool   : if set to excited, then the current status needs to be replaced with the new status
"""
mutable struct CachedPool{O,T}
    name::Symbol
    object::Vector{O}
    current::Vector{T}
    new::Vector{T}
    version::Vector{Int128}
    excited::Vector{Bool}

    function CachedPool(name::Symbol, objType::DataType, weightType::DataType)
        object = Vector{objType}(undef, 0)
        current = Vector{weightType}(undef, 0)
        _new = Vector{weightType}(undef, 0)
        version = Vector{Int128}(undef, 0)
        excited = Vector{Bool}(undef, 0)
        return new{objType,weightType}(name, object, current, _new, version, excited)
    end
    # function CachedPool{T}(obj::Vector{O}) where {O,T}
    #     weight = zeros(5, length(obj))
    #     return new{O,T}(obj, weight)
    # end
end

Base.length(pool::CachedPool) = length(pool.object)
Base.size(pool::CachedPool) = size(pool.object)
Base.show(io::IO, pool::CachedPool) = print(io, pool.object)
# Base.view(pool::Pool, inds...) = Base.view(pool.pool, inds...)

#index interface for Pool
Base.getindex(pool::CachedPool, i) = pool.object[i]
Base.setindex!(pool::CachedPool, v, i) = setindex!(pool.object, v, i)
Base.firstindex(pool::CachedPool) = 1
Base.lastindex(pool::CachedPool) = length(pool)

# iterator interface
function Base.iterate(pool::CachedPool)
    if length(pool) == 0
        return nothing
    else
        return (pool.object[1], 1)
    end
end

function Base.iterate(pool::CachedPool, state)
    if state >= length(pool) || length(pool) == 0
        return nothing
    else
        return (pool.object[state+1], state + 1)
    end
end

function append(pool::CachedPool, object)
    # @assert para isa eltype(pool.pool)
    # for (oi, o) in enumerate(pool.object)
    #     if o == object
    #         return oi #existing obj
    #     end
    # end

    id = length(pool.object) + 1
    push!(pool.object, object)

    return id #new momentum
end

function initialize!(pool::CachedPool{O,T}) where {O,T}
    N = length(pool)
    pool.current = zeros(T, N)
    pool.new = zeros(T, N)
    pool.version = ones(T, N)
    pool.excited = Vector{Bool}(zeros(Int, N))
end

function updateAll(pool::CachedPool, ignoreCache::Bool, eval::Function; kwargs...)
    N = length(pool.object)
    if ignoreCache
        T = eltype(pool.current)
        for (idx, o) in enumerate(pool.object)
            pool.current[idx] = T(eval(obj; kwargs...))
        end
    else
        error("not implemented!")
        # pool.version .= 1
        # pool.excited .= false
    end
end

"""
    struct LoopPool{T}

    Pool of loop basis. Each loop basis corresponds to a loop variable.
    A loop variable is a linear combination of N independent loops. The combination coefficients is what we call a loop basis.
    For example, if a loop is a momentum K, then

    varibale_i = K_1*basis[1, i] + K_2*basis[2, i] + K_3*basis[3, i] + ...

# Members
- name::Symbol : name of the pool
- dim::Int     : dimension of a loop variable (for example, the dimension of a momentum-frequency loop variable is (d+1) where d is the spatial dimension)
- N::Int       : number of independent loops (dimension of loop basis)
- basis::Matrix{T}    : Matrix of (N x Nb) that stores the loop basis, where Nb is the number of loop basis (or number of loop variables).
- current::Matrix{T}  : Matrix of (dim x Nb) that stores the loop variables, where Nb is the number of loop basis (or number of loop variables).
"""
mutable struct LoopPool{T}
    name::Symbol
    dim::Int #dimension
    N::Int #number of basis
    basis::Matrix{T}
    current::Matrix{T}

    function LoopPool(name::Symbol, dim::Int, N::Int, type::DataType = Float64)
        basis = Matrix{type}(undef, N, 0) # Nx0 matrix
        current = Matrix{type}(undef, dim, 0) # dimx0 matrix
        return new{type}(name, dim, N, basis, current)
    end
end

Base.length(pool::LoopPool) = size(pool.basis)[2]
Base.size(pool::LoopPool) = length(pool)
Base.show(io::IO, pool::LoopPool) = print(io, pool.basis)
# Base.view(pool::Pool, inds...) = Base.view(pool.pool, inds...)

#index interface for Pool
Base.getindex(pool::LoopPool, i) = pool.basis[:, i]
Base.setindex!(pool::LoopPool, v, i) = setindex!(pool.basis, v, i)
Base.firstindex(pool::LoopPool) = 1
Base.lastindex(pool::LoopPool) = length(pool)

# iterator interface
function Base.iterate(pool::LoopPool)
    if length(pool) == 0
        return nothing
    else
        return (pool.basis[:, 1], 1)
    end
end

function Base.iterate(pool::LoopPool, state)
    if state >= length(pool) || length(pool) == 0
        return nothing
    else
        return (pool.basis[:, state+1], state + 1)
    end
end

function update(pool::LoopPool, variable = rand(eltype(pool.current), pool.dim, pool.N))
    loopNum = size(pool.basis)[1]
    pool.current = view(variable, :, 1:loopNum) * pool.basis
    # B = view(pool.basis, 1:loopNum, :)
    # B = view(pool.basis, 1:loopNum, :)
    # C = pool.current[:, 1:length(pool)]
    # LinearAlgebra.mul!(pool.current, A, pool.basis)
    # LinearAlgebra.BLAS.gemm!('N', 'N', false, variable[:, 1:loopNum], pool.basis[1:loopNum, 1:length(pool)], false, pool.current[:, 1:length(pool)])
end

# current(pool::LoopPool, idx) = pool.current[:, idx]
current(pool::LoopPool, idx) = view(pool.current, :, idx)

function append(pool::LoopPool, basis::AbstractVector)
    for bi in 1:length(pool)
        if pool.basis[:, bi] ≈ basis
            return bi
        end
    end

    pool.basis = hcat(pool.basis, basis)
    pool.current = hcat(pool.current, rand(eltype(pool.current), pool.dim))
    # pool.size += 1
    # @assert pool.size <= size(pool.basis)[2] "Too many loop basis!. Increase maxSize when creates the LoopPool!"
    # pool.basis[:, pool.size] = basis
    return length(pool)
end


text:
__precompile__()
module LoewnerJohnEllipsoids

using Convex
export inner_ellipsoid,
	   outer_ellipsoid,
	   box,
	   barrier_value

type LinearConstraint{T<:Real}
	# Represents the constraint a * x <= b
	a::Vector{T}
	b::T
end

type Ellipsoid{T}
	# Represents an ellipsoid as the set
	# { x | x = C*u + d, ||u||_2 <= 1 }
	# or the set
	# { x | ||(x - d)' * P' * P * (x - d) <= 1 }
	# where P = inv(C)' * inv(C)
	C::Array{T, 2}
	d::Array{T, 1}
	P::Array{T, 2}
end

Ellipsoid{T}(C::Array{T, 2}, d::Array{T, 1}) = begin
	C_inv = inv(C)
	Ellipsoid(C, d, C_inv' * C_inv)
end

Ellipsoid{T}(C::Array{T,2}, d::Array{T,2}) = Ellipsoid(C, vec(d))

function barrier_value{T}(ell::Ellipsoid{T}, x)
	# Evaluates the ellipsoid's barrier function, which is defined to be < 0
	# inside the ellipsoid, 0 on the boundary, and > 0 outside
	return ((x - ell.d)' * ell.P * (x - ell.d))[1] - 1
end

type Polyhedron{T}
	constraints::Vector{LinearConstraint{T}}
end

function box{T}(lower::Vector{T}, upper::Vector{T})
	@assert length(lower) == length(upper)
	dimension = length(lower)
	constraints = Array{LinearConstraint{T}}(2 * dimension)
	for i = 1:dimension
		a = zeros(T, dimension)
		a[i] = 1
		constraints[i] = LinearConstraint(a, upper[i])
	end
	for i = 1:dimension
		a = zeros(T, dimension)
		a[i] = -1
		constraints[i + dimension] = LinearConstraint(a, -lower[i])
	end
	Polyhedron(constraints)
end

function inner_ellipsoid{T<:Real}(polyhedron::Polyhedron{T})
	@assert length(polyhedron.constraints) > 2
	dimension = length(polyhedron.constraints[1].a)
	C = Variable(dimension, dimension)
	d = Variable(dimension)
	problem = maximize(logdet(C))
	for constraint in polyhedron.constraints
		problem.constraints += norm(C*constraint.a) <= constraint.b - dot(constraint.a, d)
	end
	solve!(problem)
	return Ellipsoid(C.value, d.value)
end

function outer_ellipsoid{T}(points::Vector{T})
	@assert length(points) > 1
	dimension = length(points[1])
	P = Variable(dimension, dimension)
	q = Variable(dimension)
	problem = maximize(logdet(P))
	for point in points
		problem.constraints += norm(P*point - q) <= 1
	end
	solve!(problem)
	P, q = P.value, q.value
	C = inv(P)
	return Ellipsoid(C, C * q)
end

end # module

ext: jl


text:
# while waiting for this PR to get merged: https://github.com/Ismael-VC/Suppressor.jl/pull/12
# we'll just include the relevant code here

macro suppress_err(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            ORIGINAL_STDERR = STDERR
            err_rd, err_wr = redirect_stderr()
            err_reader = @async readstring(err_rd)
        end

        value = $(esc(block))

        if ccall(:jl_generating_output, Cint, ()) == 0
            redirect_stderr(ORIGINAL_STDERR)
            close(err_wr)
        end
        value
    end
end


text:
# SIAR model
# https://arxiv.org/pdf/2012.00443.pdf
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    S'(t) = -S(t) * (betaIa * Ia(t) + betaIs * Is(t) + betaH * H(t) + betaT * T(t)) * Ninv,
    E'(t) = S(t) * (betaIa * Ia(t) + betaIs * Is(t) + betaH * H(t) + betaT * T(t)) * Ninv - (alphaEIa + alphaEIs) * E(t),
    Ia'(t) = alphaEIa * E(t) - (alphaIaIs + alphaIaRu) * Ia(t) - xi * Ia(t),
    Is'(t) = alphaEIs * E(t) + alphaIaIs * Ia(t) - (alphaIsH + alphaIsT + alphaIsRu + alphaIsD) * Is(t),
    H'(t) = alphaIsH * Is(t) + xi * Ia(t) - (alphaHT + alphaHRd) * H(t),
    T'(t) = alphaIsT * Is(t) + alphaHT * H(t) - (alphaTRd + alphaTD) * T(t),
    Rd'(t) = alphaHRd * H(t) + alphaTRd * T(t),
    D'(t) = alphaIsD * Is(t) + alphaTD * T(t),
    y1(t) = H(t),
    y2(t) = T(t),
    y3(t) = Rd(t),
    y4(t) = D(t),
    y5(t) = Ninv
)

@time println(assess_identifiability(ode))


text:
# Variation of Information

"""
    varinfo(k1::Int, a1::AbstractVector{Int}, k2::Int, a2::AbstractVector{Int})
    varinfo(R::ClusteringResult, k0::Int, a0::AbstractVector{Int})
    varinfo(R1::ClusteringResult, R2::ClusteringResult)

Compute the variation of information between the two clusterings.

Each clustering is provided either as an instance of [`ClusteringResult`](@ref)
subtype or as a pair of arguments:
 - a number of clusters (`k1`, `k2`, `k0`)
 - a vector of point to cluster assignments (`a1`, `a2`, `a0`).
"""
function varinfo(k1::Int, a1::AbstractVector{Int},
                 k2::Int, a2::AbstractVector{Int})

    # check input arguments
    n = length(a1)
    length(a2) == n || throw(DimensionMismatch("Inconsistent array length."))

    # count & compute probabilities
    p1 = zeros(k1)
    p2 = zeros(k2)
    P = zeros(k1, k2)

    for i = 1:n
        @inbounds l1 = a1[i]
        @inbounds l2 = a2[i]
        p1[l1] += 1.0
        p2[l2] += 1.0
        P[l1, l2] += 1.0
    end

    for i = 1:k1
        @inbounds p1[i] /= n
    end
    for i = 1:k2
        @inbounds p2[i] /= n
    end
    for i = 1:(k1*k2)
        @inbounds P[i] /= n
    end

    # compute variation of information

    H1 = entropy(p1)
    H2 = entropy(p2)

    I = 0.0
    for j = 1:k2, i = 1:k1
        pi = p1[i]
        pj = p2[j]
        pij = P[i,j]
        if pij > 0.0
            I += pij * log(pij / (pi * pj))
        end
    end

    return H1 + H2 - I * 2.0
end

varinfo(R::ClusteringResult, k0::Int, a0::AbstractVector{Int}) =
    varinfo(nclusters(R), assignments(R), k0, a0)

varinfo(R1::ClusteringResult, R2::ClusteringResult) =
    varinfo(nclusters(R1), assignments(R1),
            nclusters(R2), assignments(R2))


text:
module GlobalSensitivityAnalysis

using DataStructures
using Distributions
using Statistics
using Sobol

export
    sample, analyze, SobolData

include("utils.jl")
include("sample_sobol.jl")
include("analyze_sobol.jl")
include("test_functions/ishigami.jl")

end # module


text:
using FileIO
using FilePathsBase
using Test

Threads.nthreads() <= 1 && @info "Threads.nthreads() = $(Threads.nthreads()), multithread tests will be disabled"

# Both FileIO and FilePathsBase export filename, but we only want the FileIO definition.
using FileIO: filename

struct MimeSaveTestType
end

@testset "FileIO" begin
    include("query.jl")
    include("loadsave.jl")
    include("error_handling.jl")
    include("test_mimesave.jl")
end


text:
using Base.Test

addprocs(3)

println("[TEST] with local Graphs.jl dictionary and arrays only (multicore)...")
include("fourdoortest.jl")
println("Success")
@test true

println("[TEST] plot functions...")
using Gadfly
# draw all beliefs
DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
@test true

print("[TEST] Ensure memory return is working properly...")
include("typeReturnMemRef.jl")
println("Success")

println("[TEST] packing converters work...")
# using fourdoortest data
topack = fg.f[4].attributes["data"]
dd = convert(FunctionNodeData{PackedOdo},topack)
upd = convert(FunctionNodeData{Odo}, dd)
@test topack.fnc.Zij[1] == upd.fnc.Zij[1]

# data structure conversion tests for protobuffing
println("Testing conversion to packed data structure and back")
dat = IncrementalInference.dlapi.getvertex(fg,1).attributes["data"]
pd = convert(PackedVariableNodeData, dat) #fg.v[1].attributes["data"]
# pd = convert(PackedVariableNodeData,fg.v[1].attributes["data"])
unpckd = convert(VariableNodeData, pd)

@test compare(IncrementalInference.dlapi.getvertex(fg,1).attributes["data"], unpckd)
@test IncrementalInference.dlapi.getvertex(fg,1).attributes["data"] == unpckd
println("Conversions and comparisons agree")

if false
  println("[TEST] with CloudGraphs data layer (multicore)...")
  include("fourdoortestcloudgraph.jl")
  println("[SUCESS]")
else
  warn("[NOT TESTING] CloudGraphs interface, which needs Neo4j installed")
end

println("[TEST] Pose3 evaluations...")
include("threeDimLinearProductTest.jl")
println("[SUCCESS]")


text:
# This file was generated, do not modify it. # hide
fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using RecipesBase

@recipe function f(R::T) where {T<:PLMClusteringResult}
    χ = isinf(R.ϵ) ? 2.0 : R.ϵ
    for (i,idxs) in enumerate(R.clusters)
        addlabel = true
        for c in idxs
            @series begin
                label --> (addlabel ? "MC$i" : "")
                linecolor --> i
                models(R)[c], χ
            end
            addlabel = false
        end
    end
    if length(size(R.complex)) > 0
        D = hcat(map(mean, models(R))...)'
        R.complex, D
    end
end

max_stars_repo_name: wildart/PLMC

max_issues_repo_name: wildart/PLMC

max_forks_repo_name: wildart/PLMC


text:
############
# bindings #
############

const CONTEXT_BINDING = Symbol("__CONTEXT__")
const METADATA_BINDING = Symbol("__METADATA__")
const CONFIG_BINDING = Symbol("__trace__")

############
# @context #
############

"""
    Cassette.@context Ctx

Defined and return a new Cassette context type with the name `Ctx`.
"""
macro context(Ctx)
    @assert isa(Ctx, Symbol) "context name must be a Symbol"
    name = Expr(:quote, Ctx)
    return esc(quote
        struct $Ctx{T<:$Cassette.Tag} <: $Cassette.Context
            @inline $Ctx(f) = new{$Cassette.tagtype(f)}()
        end
        $Cassette.@primitive function Cassette.tagtype(x) where {__CONTEXT__<:$Ctx}
            return $Cassette.tagtype(__CONTEXT__, x)
        end
        $Cassette.@execution function Core.getfield(x, name) where {__CONTEXT__<:$Ctx}
            return $Cassette._getfield(x, name)
        end
        $Cassette.@execution function Core.setfield!(x, name, y) where {__CONTEXT__<:$Ctx}
            return $Cassette._setfield!(x, name, y)
        end
        $Cassette.@execution function Core.arrayref(check, x, i) where {__CONTEXT__<:$Ctx}
            return $Cassette._arrayref(check, x, i)
        end
        $Cassette.@execution function Core.arrayset(check, x, y, i) where {__CONTEXT__<:$Ctx}
            return $Cassette._arrayset(check, x, y, i)
        end
        $Ctx
    end)
end

#########
# @pass #
#########

"""
    Cassette.@pass transform

Return a Cassette pass that applies `transform` to all overdubbed method bodies. `transform`
must be callable with the following signature:

    transform(signature::Type{Tuple{...}}, method_body::CodeInfo)::CodeInfo

Note that this macro expands to an `eval` call and thus should only be called at top-level.
Furthermore, to avoid world-age issues, `transform` should not be overloaded after it has
been registered with `@pass`.
"""
macro pass(transform)
    Pass = gensym("PassType")
    name = Expr(:quote, :($__module__.$Pass))
    line = Expr(:quote, __source__.line)
    file = Expr(:quote, __source__.file)
    return esc(quote
        struct $Pass <: $Cassette.AbstractPass end
        (::Type{$Pass})(signature, codeinfo) = $transform(signature, codeinfo)
        eval($Cassette, $Cassette.overdub_transform_call_definition($name, $line, $file))
        $Pass()
    end)
end

#####################
# @prehook/posthook #
#####################

"""
    Cassette.@prehook contextual_method_definition

Place in front of a contextual method definition to define a callback that Cassette
executes before method calls that match the contextual method definition's signature.

For example, the following code uses `@prehook` to increment a counter stored in
`__trace__.metadata` every time a method is called with one or more arguments matching a
given type:

    using Cassette: @context, @prehook

    @context CountCtx

    mutable struct Count{T}
        count::Int
    end

    @prehook function (f::Any)(input::T, ::T...) where {T,__CONTEXT__<:CountCtx,__METADATA__<:Count{T}}
        __trace__.metadata.count += 1
    end

Prehooks are generally useful for inserting side-effects that do not ultimately affect
the trace's execution path. However, prehooks are allowed to (for example) mutate input
argument state.

For details regarding the format of `contextual_method_definition`, see the contextual
dispatch documentation.
"""
macro prehook(method)
    return contextual_definition!(:($Cassette.prehook), method)
end

"""
    Cassette.@posthook contextual_method_definition

Place in front of a contextual method definition to define a callback that Cassette executes
after method calls that match the contextual method definition's signature. Note that the
contextual method definition's signature differs from those of the matched method calls in
that the first argument is the method calls' output.

For example, the following code uses `@posthook` to increment a counter stored in
`__trace__.metadata` every time a method call's output type matches the type of its
input:

    using Cassette: @context, @posthook

    @context CountCtx

    mutable struct Count{T}
        count::Int
    end

    @posthook function (f::Any)(output::T, input::T, ::T...) where {T,__CONTEXT__<:CountCtx,__METADATA__<:Count{T}}
        __trace__.metadata.count += 1
    end

Posthooks are generally useful for inserting side-effects that do not ultimately affect
the trace's execution path. However, posthooks are allowed to (for example) mutate
input/output argument state.

For details regarding the format of `contextual_method_definition`, see the contextual
dispatch documentation.
"""
macro posthook(method)
    return contextual_definition!(:($Cassette.posthook), method)
end

##############
# @execution #
##############

"""
    Cassette.@execution contextual_method_definition

Place in front of a contextual method definition to overload Cassette's primitive execution
behavior for method calls matching the contextual method definition's signature. Note that
this execution behavior does not apply to non-primitives, which, by definition, Cassette
will simply trace through.

For details regarding the format of `contextual_method_definition`, see the contextual
dispatch documentation.
"""
macro execution(method)
    return contextual_definition!(:($Cassette.execution), method)
end

################
# @isprimitive #
################

"""
    Cassette.@isprimitive contextual_method_signature

Place in front of a contextual method signature to mark matching method calls as Cassette
primitives. Cassette primitives are executed using the contextual method definitions
provided via `Cassette.@execution`.

For details regarding the format of `contextual_method_signature`, see the contextual
dispatch documentation.
"""
macro isprimitive(signature)
    body = Expr(:block)
    push!(body.args, :(return true))
    return contextual_definition!(:($Cassette.is_user_primitive), signature, body)
end

##############
# @primitive #
##############

"""
    Cassette.@primitive contextual_method_definition

A convenience macro for simultaneously applying `Cassette.@execution` and
`Cassette.@isprimitive` to the provided contextual method definition.
"""
macro primitive(method)
    @assert is_method_definition(method)
    signature = deepcopy(first(method.args))
    return esc(quote
        $Cassette.@execution $method
        $Cassette.@isprimitive $signature
    end)
end

############
# @overdub #
############

"""
    Cassette.@overdub Ctx expression

A convenience macro for overdubbing and executing `expression` within the context `Ctx`.
"""
macro overdub(ctx, ex)
  :(overdub($(esc(ctx)), () -> $(esc(ex)))())
end

########
# @Box #
########

"""
    Cassette.@Box([V, M])

Used within contextual method signatures to indicate the type of a value `V` wrapped in a
`Cassette.Box` with metadata of type `M`. `V` and `M` default to `Any` if unspecified.

For example:

    using Cassette: @context, @prehook, @Box

    @context BoxCtx

    @prehook function (f::Any)(x::@Box) where {__CONTEXT__<:BoxCtx}
        println("The value ", unbox(__trace__.context, x),
                " was wrapped in a Cassette.Box and passed to ", f)
    end

    @prehook function (f::Any)(x::@Box(Number)) where {__CONTEXT__<:BoxCtx}
        println("The numeric value ", unbox(__trace__.context, x),
                " was wrapped in a Cassette.Box and passed to ", f)
    end

    @prehook function (f::Any)(x::@Box(Number,String)) where {__CONTEXT__<:BoxCtx}
        println("The numeric value ", unbox(__trace__.context, x),
                " was wrapped in a Cassette.Box and passed to ", f,
                " and carries the message: ", meta(__trace__.context, x))
    end

Note that `Cassette.@Box` can only used within the scope of a contextual method signature.

For further details, see the contextual metadata propagation documentation.
"""
macro Box(args...)
    error("`Cassette.@Box` can only used within the scope of a contextual method signature")
end

#############
# utilities #
#############

function typify_signature(signature)
    if isa(signature, Expr)
        f, args = signature.args[1], signature.args[2:end]
        if signature.head == :where
            return Expr(:where, typify_signature(f), args...)
        elseif signature.head == :call
            if !(isa(f, Expr) && f.head == :(::))
                return Expr(:call, :(::typeof($f)), args...)
            end
            return signature
        end
    end
    error("malformed signature: $signature")
end

function contextual_definition!(f, method)
    @assert is_method_definition(method)
    signature, body = method.args
    return contextual_definition!(f, signature, body)
end

function contextual_definition!(f, signature::Expr, body::Expr)
    @assert(signature.head == :where,
            "method signature missing `where` clause; `$(CONTEXT_BINDING) <: ContextType` "*
            "must be defined in your method definition's `where` clause")

    signature = typify_signature(signature)

    # add metadata typevar (if not already present)
    has_metadata_typevar = any(signature.args) do x
        x == METADATA_BINDING || isa(x, Expr) && x.head == :(<:) && x.args[1] == METADATA_BINDING
    end
    !has_metadata_typevar && push!(signature.args, METADATA_BINDING)

    # add world age typevar
    world_binding = gensym("world")
    push!(signature.args, world_binding)

    config_arg = :($CONFIG_BINDING::$Cassette.TraceConfig{$CONTEXT_BINDING,$METADATA_BINDING,$world_binding})

    call_args = signature.args[1].args
    for i in 1:length(call_args)
        x = call_args[i]
        if Base.Meta.isexpr(x, :(::))
            xtype = last(x.args)
            if Base.Meta.isexpr(xtype, :macrocall) && xtype.args[1] == Symbol("@Box")
                box_args = xtype.args[3:end]
                if isempty(box_args)
                    U, M = :Any, :Any
                elseif length(box_args) == 1
                    U, M = first(box_args), :Any
                elseif length(box_args) == 2
                    U, M = box_args
                else
                    error("incorrect usage of `@Box`: $(xtype)")
                end
                new_xtype = :($Cassette.Box{$CONTEXT_BINDING,<:$U,<:Any,<:$M})
            else
                new_xtype = :(Union{$Cassette.Box{<:Any,<:$xtype},$xtype})
            end
            x.args[end] = new_xtype
        end
    end

    signature.args[1] = Expr(:call, f, config_arg, call_args...)

    pushfirst!(body.args, Expr(:meta, :inline))

    return esc(Expr(:function, signature, body))
end


text:
Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[]) # JuliaLang/julia/pull/28625

using CLIMA, Documenter, Literate

generated_dir = joinpath(@__DIR__, "src", "generated") # generated files directory
mkpath(generated_dir)

include("list_of_experiments.jl")        # defines a dict `experiments`
include("list_of_tutorials.jl")          # defines a dict `tutorials`
include("list_of_extending_clima.jl")    # defines a dict `extending_clima`
include("list_of_discussions.jl")        # defines a dict `apis`
include("list_of_apis.jl")               # defines a dict `discussions`

pages = Any[
    "Home" => "index.md",
    "Installation" => "Installation.md",
    "Experiments" => experiments,
    "Tutorials" => tutorials,
    "Extending CLIMA" => extending_clima,
    "APIs" => apis,
    "Theory & design philosophy" => discussions,
    # TODO: Move everything below here into one of the above sections
    "Developer docs" => Any[
        "CodingConventions.md",
        "AcceptableUnicode.md",
        "VariableList.md",
    ],
]

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

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
    # prettyurls = !("local" in ARGS),
    # canonical = "https://climate-machine.github.io/CLIMA/stable/",
)

makedocs(
    sitename = "CLIMA",
    doctest = false,
    strict = false,
    linkcheck = false,
    format = format,
    checkdocs = :exports,
    # checkdocs = :all,
    clean = true,
    modules = [Documenter, CLIMA],
    pages = pages,
)

include("clean_build_folder.jl")

deploydocs(
    repo = "github.com/climate-machine/CLIMA.git",
    target = "build",
    push_preview = true,
)


text:
module Assets

using ..Ahorn
using Cairo

const crowIdle = open(Cairo.read_from_png, Ahorn.abs"../assets/crow18.png")
const crowPeck1 = open(Cairo.read_from_png, Ahorn.abs"../assets/crow20.png")
const crowPeck2 = open(Cairo.read_from_png, Ahorn.abs"../assets/crow21.png")

const risingLava = open(Cairo.read_from_png, Ahorn.abs"../assets/rising_lava.png")
const lavaSandwich = open(Cairo.read_from_png, Ahorn.abs"../assets/lava_sandwich.png")
const speechBubble = open(Cairo.read_from_png, Ahorn.abs"../assets/core_message.png")
const speaker = open(Cairo.read_from_png, Ahorn.abs"../assets/sound_source.png")
const tentacle = open(Cairo.read_from_png, Ahorn.abs"../assets/tentacles.png")
const northernLights = open(Cairo.read_from_png, Ahorn.abs"../assets/northern_lights.png")
const summitBackgroundManager = open(Cairo.read_from_png, Ahorn.abs"../assets/summit_background_manager.png")
const cutsceneNode = open(Cairo.read_from_png, Ahorn.abs"../assets/cutscene_node.png")

const missingImage = open(Cairo.read_from_png, Ahorn.abs"../assets/missing_image.png")

# For backwards compability... Spelling words is sometimes too hard...
const lavaSanwitch = lavaSandwich

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# Scatter plots

function PlotScatter(label_id, values::Union{AbstractArray{T},Ptr{T},Ref{T}}, count::Integer,
                     xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0,
                     stride::Integer = sizeof(T)) where {T<:ImPlotData}
    ImPlot_PlotScatter(label_id, values, count, xscale, x0, offset, stride)
end

function PlotScatter(label_id, values::Union{AbstractArray{T},Ptr{T},Ref{T}}, count::Integer,
                     xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0,
                     stride::Integer = sizeof(Float64)) where {T<:Real}
    ImPlot_PlotScatter(label_id, Float64.(values), count, xscale, x0, offset, stride)
end

function PlotScatter(label_id, x::Union{AbstractArray{T},Ptr{T},Ref{T}},
                     y::Union{AbstractArray{T},Ptr{T},Ref{T}}, count::Integer,
                     offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    ImPlot_PlotScatter(label_id, x, y, count, offset, stride)
end

function PlotScatter(label_id, x::Union{AbstractArray{T},Ptr{T},Ref{T}},
                     y::Union{AbstractArray{T},Ptr{T},Ref{T}}, count::Integer,
                     offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    ImPlot_PlotScatter(label_id, Float64.(x), Float64.(y), count, offset, stride)
end

function PlotScatter(label_id::String, x::AbstractArray{T}, y::AbstractArray{T};
                  count::Integer = min(length(x), length(y)),
                  offset::Integer = 0, stride::Integer = 1) where {T<:ImPlotData}

    ImPlot_PlotScatter(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotScatter(label_id::String, x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real,T2<:Real}
    PlotScatter(label_id::String, promote(x,y)...; kwargs...)
end

function PlotScatter(label_id::String, y::AbstractArray{T};count::Integer=length(y),
                  xscale::Real = 1.0, x0::Real = 0,
                  offset::Integer=0, stride::Integer=1) where {T<:ImPlotData}
    ImPlot_PlotScatter(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotScatter(label_id::String, x::UnitRange{<:Integer}, y::AbstractArray{T};
                  xscale::Real = 1.0, x0::Real = 0,
                  ) where {T<:ImPlotData}

    count::Cint = length(x)
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(sizeof(T))
    ImPlot_PlotScatter(label_id, y, count, xscale, x0, offset, stride)
end

function PlotScatter(label_id::String, x::StepRange, y::AbstractArray{T};
                     xscale::Real = 1.0, x0::Real = 0,
                     ) where {T<:ImPlotData}

    x.stop < 1 && throw("Range out of bounds")
    count::Cint = length(x)
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(x.step * sizeof(T))
    ImPlot_PlotScatter(label_id, y, count, xscale, x0, offset, stride)
end

# xfield, yfield should be propertynames of eltype(structvec)
function PlotScatter(
    label_id::String,
    structvec::Vector{T}, xfield::Symbol, yfield::Symbol;
    count::Integer = length(structvec), offset::Integer = 0,
    stride::Integer = 1
) where T

    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = (pointer(structvec, 1) + x_offset) |> Ptr{Tx}
    y_ptr = (pointer(structvec, 1) + y_offset) |> Ptr{Ty}
    p = pointer(structvec, 1)

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own = false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own = false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    ImPlot_PlotScatter(label_id, x, y, count, offset, stride)
end

PlotScatterG(label_id, getter, data, count, offset=0) =
ImPlot_PlotScatterG(label_id, getter, data, count, offset)

line_functions = [
    (Cfloat, Cfloat) => ImPlot_PlotScatter_FloatPtrFloatPtr
    (Cdouble, Cdouble) => ImPlot_PlotScatter_doublePtrdoublePtr
    (Int32, Int32) => ImPlot_PlotScatter_S32PtrS32Ptr
    (UInt64, UInt64) => ImPlot_PlotScatter_U64PtrU64Ptr
    (UInt16, UInt16) => ImPlot_PlotScatter_U16PtrU16Ptr
    (Int64, Int64) => ImPlot_PlotScatter_S64PtrS64Ptr
    (UInt8, UInt8) => ImPlot_PlotScatter_U8PtrU8Ptr
    (Int16, Int16) => ImPlot_PlotScatter_S16PtrS16Ptr
    (UInt32, UInt32) => ImPlot_PlotScatter_U32PtrU32Ptr
    (Int8, Int8) => ImPlot_PlotScatter_S8PtrS8Ptr
]
for ((T1, T2), func) in line_functions
    @eval ImPlot_PlotScatter(label_id::String, xs::AbstractArray{<:$T1}, ys::AbstractArray{<:$T2}, count, offset, stride) = $func(label_id, xs, ys, count, offset, stride)
end

line_functions_nox = [
    UInt16 => ImPlot_PlotScatter_U16PtrInt
    Cfloat => ImPlot_PlotScatter_FloatPtrInt
    Int64 => ImPlot_PlotScatter_S64PtrInt
    UInt8 => ImPlot_PlotScatter_U8PtrInt
    Int16 => ImPlot_PlotScatter_S16PtrInt
    UInt32 => ImPlot_PlotScatter_U32PtrInt
    Int8 => ImPlot_PlotScatter_S8PtrInt
    Cdouble => ImPlot_PlotScatter_doublePtrInt
    Int32 => ImPlot_PlotScatter_S32PtrInt
    UInt64 => ImPlot_PlotScatter_U64PtrInt
]
for (T1, func) in line_functions_nox
    @eval ImPlot_PlotScatter(label_id::String, values::AbstractArray{<:$T1}, count, xscale=1.0, x0=0, offset=0, stride=$(sizeof(T1))) = $func(label_id, values, count, xscale, x0, offset, stride)
end

# ImPlot_PlotScatterG                   ImPlot_PlotScatter_S64PtrInt           ImPlot_PlotScatter_U32PtrU32Ptr
# ImPlot_PlotScatter_FloatPtrFloatPtr   ImPlot_PlotScatter_S64PtrS64Ptr        ImPlot_PlotScatter_U64PtrInt
# ImPlot_PlotScatter_FloatPtrInt        ImPlot_PlotScatter_S8PtrInt            ImPlot_PlotScatter_U64PtrU64Ptr
# ImPlot_PlotScatter_S16PtrInt          ImPlot_PlotScatter_S8PtrS8Ptr          ImPlot_PlotScatter_U8PtrInt
# ImPlot_PlotScatter_S16PtrS16Ptr       ImPlot_PlotScatter_U16PtrInt           ImPlot_PlotScatter_U8PtrU8Ptr
# ImPlot_PlotScatter_S32PtrInt          ImPlot_PlotScatter_U16PtrU16Ptr        ImPlot_PlotScatter_doublePtrInt
# ImPlot_PlotScatter_S32PtrS32Ptr       ImPlot_PlotScatter_U32PtrInt           ImPlot_PlotScatter_doublePtrdoublePtr


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

solves a simple random linear optimization problem (LP):
    min  c'x
    s.t. Ax = b, x >= 0
=#

using SparseArrays

struct LinearOptNative{T <: Real} <: ExampleInstanceNative{T}
    m::Int
    n::Int
    nz_frac::Float64
end

function build(inst::LinearOptNative{T}) where {T <: Real}
    (m, n, nz_frac) = (inst.m, inst.n, inst.nz_frac)
    @assert 0 < nz_frac <= 1

    # generate random data
    A = (nz_frac >= 1) ? rand(T, m, n) : sprand(T, m, n, nz_frac)
    A .*= 10
    b = vec(sum(A, dims = 2))
    c = rand(T, n)
    G = Diagonal(-one(T) * I, n) # TODO uniformscaling
    h = zeros(T, n)
    cones = Cones.Cone{T}[Cones.Nonnegative{T}(n)]

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
"""
Test orthonormalization of wavefunction in real space
"""
function test_ortho( dVol, Nstates, psi )
    print("\nTest norm:\n")
    print( "===========================\n")
    for is = 1:Nstates
        print("State: #$(is): $(dot( psi[:,is], psi[:,is] ) * dVol)\n")
    end
    print("\nTest ortho w.r.t state #1:\n")
    print( "===========================\n")
    for is = 2:Nstates
        print("State: #$(is): $(dot( psi[:,is], psi[:,1] ) * dVol)\n")
    end
    print("\n")
end


text:
"""
    ccm_with_summary(source,
            target,
            timeseries_lengths;
            average_measure::Symbol = :median,
            uncertainty_measure::Symbol = :quantile,
            quantiles = [0.327, 0.673],
            kwargs...)

## Algorithm

Compute the cross mapping between a `source` series and a `target` series over
different `timeseries_lengths` and return summary statistics of the results.

## Arguments
- **`source`**: The data series representing the putative source process.
- **`target`**: The data series representing the putative target process.
- **`timeseries_lengths`**: Time series length(s) for which to compute the
    cross mapping(s).

## Summary keyword arguments
- **`average_measure`**: Either `:median` or `:mean`. Default is `:median`.
- **`uncertainty_measure`**: Either `:quantile` or `:std`. Default is `:quantile`.
- **`quantiles`**: Compute uncertainty over quantile(s) if `uncertainty_measure`
    is `:quantile`. Default is `[0.327, 0.673]`, roughly corresponding to 1s for
    normally distributed data.

## Keyword arguments to `crossmap`
- **`dim`**: The dimension of the state space reconstruction (delay embedding)
    constructed from the `target` series. Default is `dim = 3`.
- **`τ`**: The embedding lag for the delay embedding constructed from `target`.
    Default is `τ = 1`.
- **`η`**: The prediction lag to use when predicting scalar values of `source`
    from the delay embedding of `target`.
    `η > 0` are forward lags (causal; `source`'s past influences `target`'s future),
    and `η < 0` are backwards lags (non-causal; `source`'s' future influences
    `target`'s past). Adjust the prediction lag if you
    want to performed lagged ccm
    [(Ye et al., 2015)](https://www.nature.com/articles/srep14750).
    Default is `η = 0`, as in
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079).
    *Note: The sign of the lag `η` is organized to conform with the conventions in
    [TransferEntropy.jl](), and is opposite to the convention used in the
    [`rEDM`](https://cran.r-project.org/web/packages/rEDM/index.html) package
    ([Ye et al., 2016](https://cran.r-project.org/web/packages/rEDM/index.html)).*
- **`libsize`**: Among how many delay embedding points should we sample time indices
    and look for nearest neighbours at each cross mapping realization (of which there
    are `n_reps`)?
- **`n_reps`**: The number of times we draw a library of `libsize` points from the
    delay embedding of `target` and try to predict `source` values. Equivalently,
    how many times do we cross map for this value of `libsize`?
    Default is `n_reps = 100`.
- **`replace`**: Sample delay embedding points with replacement? Default is `replace = true`.
- **`theiler_window`**: How many temporal neighbors of the delay embedding
    point `target_embedding(t)` to exclude when searching for neighbors to
    determine weights for predicting the scalar point `source(t + η)`.
    Default is `theiler_window = 0`.
- **`tree_type`**: The type of tree to build when looking for nearest neighbors.
    Must be a tree type from NearestNeighbors.jl. For now, this is either
    `BruteTree`, `KDTree` or `BallTree`.
- **`distance_metric`**: An instance of a `Metric` from Distances.jl. `BallTree` and `BruteTree` work with any `Metric`.
    `KDTree` only works with the axis aligned metrics `Euclidean`, `Chebyshev`,
    `Minkowski` and `Cityblock`. Default is `metric = Euclidean()` *(note the instantiation of the metric)*.
- **`correspondence_measure`**: The function that computes the correspondence
    between actual values of `source` and predicted values. Can be any
    function returning a similarity measure between two vectors of values.
    Default is `correspondence_measure = StatsBase.cor`, which returns values on ``[-1, 1]``.
    In this case, any negative values are usually filtered out (interpreted as zero coupling) and
    a value of ``1`` means perfect prediction.
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
    also proposes to use the root mean square deviation, for which a value of ``0`` would
    be perfect prediction.

## References
Sugihara, George, et al. "Detecting causality in complex ecosystems."
Science (2012): 1227079.
[http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)

Ye, Hao, et al. "Distinguishing time-delayed causal interactions using convergent cross mapping." Scientific Reports 5 (2015): 14750.
[https://www.nature.com/articles/srep14750](https://www.nature.com/articles/srep14750)

Ye, H., et al. "rEDM: Applications of empirical dynamic modeling from time series." R Package Version 0.4 7 (2016).
[https://cran.r-project.org/web/packages/rEDM/index.html](https://cran.r-project.org/web/packages/rEDM/index.html)
"""
function ccm_with_summary(source,
            target,
            timeseries_lengths;
            average_measure::Symbol = :median,
            uncertainty_measure::Symbol = :quantile,
            quantiles = [0.327, 0.673],
            kwargs...)

    if average_measure ∈ [:median, :mean]
            average = zeros(Float64, length(timeseries_lengths))
        end

        if uncertainty_measure == :quantile
            uncertainties = zeros(Float64, length(timeseries_lengths), length(quantiles))
        elseif uncertainty_measure == :std
            uncertainties = zeros(Float64, length(timeseries_lengths))
        end

        for (i, ts_len) in enumerate(timeseries_lengths)
            correlations = crossmap(source[1:ts_len], target[1:ts_len]; kwargs...)

            if average_measure == :mean
                average[i] = mean(correlations)
            elseif average_measure == :median
                average[i] = median(correlations)
            end

            if uncertainty_measure == :quantile
                for (j, quant) in enumerate(quantiles)
                    uncertainties[i, j] = quantile(correlations, quant)
                end
            elseif uncertainty_measure == :std
                uncertainties[i] = std(correlations)
            end
        end

        if average_measure == :none
            return uncertainties
        elseif uncertainty_measure == :none
            return average
        else
            return average, uncertainties
        end
end

"""
    ccm(source,
            target,
            timeseries_lengths;
            kwargs...) -> Vector{Vector{Float64}}

## Algorithm
Compute the cross mapping between a `source` series and a `target` series over
different `timeseries_lengths`.

## Arguments
- **`source`**: The data series representing the putative source process.
- **`target`**: The data series representing the putative target process.
- **`timeseries_lengths`**: Time series length(s) for which to compute the
    cross mapping(s).

## Keyword arguments to `crossmap`
- **`dim`**: The dimension of the state space reconstruction (delay embedding)
    constructed from the `target` series. Default is `dim = 3`.
- **`τ`**: The embedding lag for the delay embedding constructed from `target`.
    Default is `τ = 1`.
- **`η`**: The prediction lag to use when predicting scalar values of `source`
    fromthe delay embedding of `target`.
    `η > 0` are forward lags (causal; `source`'s past influences `target`'s future),
    and `η < 0` are backwards lags (non-causal; `source`'s' future influences
    `target`'s past). Adjust the prediction lag if you
    want to performed lagged ccm
    [(Ye et al., 2015)](https://www.nature.com/articles/srep14750).
    Default is `η = 0`, as in
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079).
    *Note: The sign of the lag `η` is organized to conform with the conventions in
    [TransferEntropy.jl](), and is opposite to the convention used in the
    [`rEDM`](https://cran.r-project.org/web/packages/rEDM/index.html) package
    ([Ye et al., 2016](https://cran.r-project.org/web/packages/rEDM/index.html)).*
- **`libsize`**: Among how many delay embedding points should we sample time indices
    and look for nearest neighbours at each cross mapping realization (of which there
    are `n_reps`)?
- **`n_reps`**: The number of times we draw a library of `libsize` points from the
    delay embedding of `target` and try to predict `source` values. Equivalently,
    how many times do we cross map for this value of `libsize`?
    Default is `n_reps = 100`.
- **`replace`**: Sample delay embedding points with replacement? Default is `replace = true`.
- **`theiler_window`**: How many temporal neighbors of the delay embedding
    point `target_embedding(t)` to exclude when searching for neighbors to
    determine weights for predicting the scalar point `source(t + η)`.
    Default is `theiler_window = 0`.
- **`tree_type`**: The type of tree to build when looking for nearest neighbors.
    Must be a tree type from NearestNeighbors.jl. For now, this is either
    `BruteTree`, `KDTree` or `BallTree`.
- **`distance_metric`**: An instance of a `Metric` from Distances.jl. `BallTree` and `BruteTree` work with any `Metric`.
    `KDTree` only works with the axis aligned metrics `Euclidean`, `Chebyshev`,
    `Minkowski` and `Cityblock`. Default is `metric = Euclidean()` *(note the instantiation of the metric)*.
- **`correspondence_measure`**: The function that computes the correspondence
    between actual values of `source` and predicted values. Can be any
    function returning a similarity measure between two vectors of values.
    Default is `correspondence_measure = StatsBase.cor`, which returns values on ``[-1, 1]``.
    In this case, any negative values are usually filtered out (interpreted as zero coupling) and
    a value of ``1`` means perfect prediction.
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
    also proposes to use the root mean square deviation, for which a value of ``0`` would
    be perfect prediction.

## References
Sugihara, George, et al. "Detecting causality in complex ecosystems."
Science (2012): 1227079.
[http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)

Ye, Hao, et al. "Distinguishing time-delayed causal interactions using convergent cross mapping." Scientific Reports 5 (2015): 14750.
[https://www.nature.com/articles/srep14750](https://www.nature.com/articles/srep14750)

Ye, H., et al. "rEDM: Applications of empirical dynamic modeling from time series." R Package Version 0.4 7 (2016).
[https://cran.r-project.org/web/packages/rEDM/index.html](https://cran.r-project.org/web/packages/rEDM/index.html)
"""
function ccm(source,
            target,
            timeseries_lengths;
            kwargs...)

    [crossmap(source[1:ts_len], target[1:ts_len]; kwargs...)
            for ts_len in timeseries_lengths]
end

"""
    convergentcrossmap(source,
            target,
            timeseries_lengths;
            summarise::Bool = true,
            average_measure::Symbol = :median,
            uncertainty_measure::Symbol = :quantile,
            quantiles = [0.327, 0.673],
            kwargs...)

## Algorithm
Compute the cross mapping between a `source` series and a `target` series over
different `timeseries_lengths`. If `summarise = true`, then call `ccm_with_summary`.
If `summarise = false`, then call `ccm` (returns raw crossmap skills).

## Arguments
- **`source`**: The data series representing the putative source process.
- **`target`**: The data series representing the putative target process.
- **`timeseries_lengths`**: Time series length(s) for which to compute the
    cross mapping(s).

## Summary keyword arguments
- **`summarise`**: Should cross map skills be summarised for each time series length?
    Default is `summarise = true`.
- **`average_measure`**: Either `:median` or `:mean`. Default is `:median`.
- **`uncertainty_measure`**: Either `:quantile` or `:std`. Default is `:quantile`.
- **`quantiles`**: Compute uncertainty over quantile(s) if `uncertainty_measure`
    is `:quantile`. Default is `[0.327, 0.673]`, roughly corresponding to 1s for
    normally distributed data.

## Keyword arguments to `crossmap`
- **`dim`**: The dimension of the state space reconstruction (delay embedding)
    constructed from the `target` series. Default is `dim = 3`.
- **`τ`**: The embedding lag for the delay embedding constructed from `target`.
    Default is `τ = 1`.
- **`η`**: The prediction lag to use when predicting scalar values of `source`
    fromthe delay embedding of `target`.
    `η > 0` are forward lags (causal; `source`'s past influences `target`'s future),
    and `η < 0` are backwards lags (non-causal; `source`'s' future influences
    `target`'s past). Adjust the prediction lag if you
    want to performed lagged ccm
    [(Ye et al., 2015)](https://www.nature.com/articles/srep14750).
    Default is `η = 0`, as in
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079).
    *Note: The sign of the lag `η` is organized to conform with the conventions in
    [TransferEntropy.jl](), and is opposite to the convention used in the
    [`rEDM`](https://cran.r-project.org/web/packages/rEDM/index.html) package
    ([Ye et al., 2016](https://cran.r-project.org/web/packages/rEDM/index.html)).*
- **`libsize`**: Among how many delay embedding points should we sample time indices
    and look for nearest neighbours at each cross mapping realization (of which there
    are `n_reps`)?
- **`n_reps`**: The number of times we draw a library of `libsize` points from the
    delay embedding of `target` and try to predict `source` values. Equivalently,
    how many times do we cross map for this value of `libsize`?
    Default is `n_reps = 100`.
- **`replace`**: Sample delay embedding points with replacement? Default is `replace = true`.
- **`theiler_window`**: How many temporal neighbors of the delay embedding
    point `target_embedding(t)` to exclude when searching for neighbors to
    determine weights for predicting the scalar point `source(t + η)`.
    Default is `theiler_window = 0`.
- **`tree_type`**: The type of tree to build when looking for nearest neighbors.
    Must be a tree type from NearestNeighbors.jl. For now, this is either
    `BruteTree`, `KDTree` or `BallTree`.
- **`distance_metric`**: An instance of a `Metric` from Distances.jl. `BallTree` and `BruteTree` work with any `Metric`.
    `KDTree` only works with the axis aligned metrics `Euclidean`, `Chebyshev`,
    `Minkowski` and `Cityblock`. Default is `metric = Euclidean()` *(note the instantiation of the metric)*.
- **`correspondence_measure`**: The function that computes the correspondence
    between actual values of `source` and predicted values. Can be any
    function returning a similarity measure between two vectors of values.
    Default is `correspondence_measure = StatsBase.cor`, which returns values on ``[-1, 1]``.
    In this case, any negative values are usually filtered out (interpreted as zero coupling) and
    a value of ``1`` means perfect prediction.
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
    also proposes to use the root mean square deviation, for which a value of ``0`` would
    be perfect prediction.

## References
Sugihara, George, et al. "Detecting causality in complex ecosystems."
Science (2012): 1227079.
[http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)

Ye, Hao, et al. "Distinguishing time-delayed causal interactions using convergent cross mapping." Scientific Reports 5 (2015): 14750.
[https://www.nature.com/articles/srep14750](https://www.nature.com/articles/srep14750)

Ye, H., et al. "rEDM: Applications of empirical dynamic modeling from time series." R Package Version 0.4 7 (2016).
[https://cran.r-project.org/web/packages/rEDM/index.html](https://cran.r-project.org/web/packages/rEDM/index.html)
"""
function convergentcrossmap(source,
            target,
            timeseries_lengths;
            summarise::Bool = true,
            average_measure::Symbol = :median,
            uncertainty_measure::Symbol = :quantile,
            quantiles = [0.327, 0.673],
            kwargs...)

    validate_average_measure(average_measure)
    validate_uncertainty_measure(uncertainty_measure)
    validate_output_selection(average_measure, uncertainty_measure, summarise)

    if summarise
        ccm_with_summary(source,
            target,
            timeseries_lengths;
            average_measure = average_measure,
            uncertainty_measure = uncertainty_measure,
            quantiles = quantiles,
            kwargs...)
    else
        ccm(source,
            target,
            timeseries_lengths;
            kwargs...)
    end

end

export
ccm_with_summary,
ccm,
convergentcrossmap


text:
module OneClassActiveLearning

using Reexport
using SVDD

using MLLabelUtils, MLKernels, MLDataUtils
using ROCAnalysis
using ValueHistories
using JSON, Unmarshal
using DataStructures
using Printf
using DelimitedFiles
using Statistics
using LinearAlgebra
using Dates
using Pkg
using JuMP
import StatsBase: countmap

using Formatting
using Memento

include("QueryStrategies/QueryStrategies.jl")
@reexport using .QueryStrategies

include("oracle.jl")
include("data_util.jl")
include("evaluate.jl")
include("result.jl")
include("serialize.jl")
include("run.jl")
include("log_util.jl")

const LOGGER = getlogger(@__MODULE__)

function __init__()
    Memento.register(LOGGER)
end

export
    load_data,
    SplitStrategy,
    FullSplitStrat, UnlabeledSplitStrat, UnlabeledAndLabeledInlierSplitStrat,
    LabeledSplitStrat, LabeledInlierSplitStrat, LabeledOutlierSplitStrat,
    DataSplits,
    get_train, get_test, get_query, calc_mask,
    get_splits_and_init_pools, get_initial_pools,
    Oracle, PoolOracle,

    ConfusionMatrix,
    cohens_kappa,
    matthews_corr,
    roc_auc,
    get_positive,
    get_negative,
    tpr,
    fpr,
    tnr,
    sensitivity,
    specificity,
    recall,
    get_n,
    f1_score,

    DataStats, Result
end


text:
module BTrDB

using HTTP, JSON
import Base: values

export
    collections, streams, # metadata queries
    refresh, stream_from_uuid, create, obliterate, settags, setannotations, # stream management functions
    nearest, earliest, latest, values, windows, aligned_windows, # data retrieval functions
    insert, delete, flush  # data management functions

###############################################################################
# Includes
###############################################################################

include("exceptions.jl")
include("types.jl")
include("utility.jl")

###############################################################################
# Module Variables & Constants
###############################################################################

APIKEY = ENV["BTRDB_API_KEY"]
ENDPOINTS = ENV["BTRDB_ENDPOINTS"]
BASEURL = "https://" * join([string.(split(ENDPOINTS, ":")[1]), "v5"], "/")

###############################################################################
# Status Related
###############################################################################

# intentionally not exported
function info()
    url = join([BASEURL, "info"], "/")
    response = apicall(url, "{}")
    data = JSON.parse(response)

    if data["stat"] != nothing
        if haskey(data["stat"], "code")
            throw(BTrDBException(data["stat"]["message"]))
        end
        throw(BTrDBException("unknown error occurred"))
    end

    return data
end

###############################################################################
# Query for Stream or Collections
###############################################################################

"""
    collections()

Returns an array of all collection strings
"""
function collections()
    return collections("")
end

"""
    collections(starts_with::String)

Returns an array of all collection strings that start with a given String
"""
function collections(starts_with::String)
    label = "collections"
    url = join([BASEURL, "listcollections"], "/")
    payload = "{\"prefix\": \"$starts_with\"}"
    payload = Dict(
        "prefix" => starts_with
    )

    response = apicall(url, JSON.json(payload))
    return parse_api_results(response, label)
end

"""
    streams(collection::String)

Returns an Array of `Stream` objects found with the supplied collection.
"""
function streams(collection::String)

    label = "results"
    url = join([BASEURL, "lookupstreams"], "/")

    payload = Dict(
        "collection" => collection,
        "isCollectionPrefix" => true,
        "tags" => [],
        "annotations" => [],
    )

    response = apicall(url, JSON.json(payload))

    return [Stream(s) for s in parse_api_results(response, label)]
end

###############################################################################
# Stream Management Functions
###############################################################################

"""
    create(uuid::String, collection::String, tags::Dict{String, String}, annotations::Dict{String, Any})

Creates a new stream on the server.
"""
function create(uuid::String, collection::String, tags::Dict{String, String}, annotations::Dict{String, Any})

    url = join([BASEURL, "create"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid),
        "collection" => collection,
        "tags" => dict2array(tags),
        "annotations" => dict2array(annotations)
    )

    response = apicall(url, JSON.json(payload))
    return refresh(uuid)
end

"""
    refresh(uuid::String)

Returns a new `Stream` object with the latest metadata from the server.
"""
function refresh(uuid::String)
    url = join([BASEURL, "streaminfo"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid),
        "omitVersion" => false,
        "omitDescriptor" => false
    )
    response = apicall(url, JSON.json(payload))
    data = JSON.parse(response)

    if data["stat"] != nothing
        if haskey(data["stat"], "msg")
            throw(BTrDBException(data["stat"]["msg"]))
        end
        println(data["stat"])
        throw(BTrDBException("unknown error occurred"))
    end

    s = Stream(data["descriptor"])
    s.version = parse(Int64, data["versionMajor"])
    s.propertyVersion = parse(Int64, data["descriptor"]["propertyVersion"])
    return s
end

"""
    refresh(stream::Stream)

Returns a new `Stream` object with the latest metadata from the server.
"""
refresh(stream::Stream) = refresh(stream.uuid)

"""
    stream_from_uuid(uuid::String)

Returns a new `Stream` object with the latest metadata from the server.
Alias for `refresh(uuid::String)`
"""
stream_from_uuid(uuid::String) = refresh(uuid)

"""
    obliterate(uuid::String)

Deletes a `Stream` with a given UUID from the server permanently.
"""
function obliterate(uuid::String)
    url = join([BASEURL, "obliterate"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid)
    )
    response = apicall(url, JSON.json(payload))
    data = JSON.parse(response)

    if data["stat"] != nothing
        if haskey(data["stat"], "code")
            throw(BTrDBException(data["stat"]["message"]))
        end
        throw(BTrDBException("unknown error occurred"))
    end
end

"""
    obliterate(stream::Stream)

Deletes a `Stream` from the server permanently.
"""
obliterate(stream::Stream) = obliterate(stream.uuid)

"""
    flush(uuid::String)

Forces the database to commit buffered data.  This function is generally
not needed by end users.
"""
function flush(uuid::String)
    url = join([BASEURL, "flush"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid)
    )
    response = apicall(url, JSON.json(payload))
    data = JSON.parse(response)

    if data["stat"] != nothing
        if haskey(data["stat"], "code")
            throw(BTrDBException(data["stat"]["message"]))
        end
        throw(BTrDBException("unknown error occurred"))
    end
end

"""
    flush(stream::Stream)

Forces the database to commit buffered data.  This function is generally
not needed by end users.
"""
flush(stream::Stream) = flush(stream.uuid)

"""
    settags(stream::Stream, tags::Dict{String, String})

Overwrites existing tags.

**Notes**

Tags are used internally by the system and are not generally meant to be
manipulated by end users.

If you would like to add your own custom metadata, please use annotations.
"""
function settags(stream::Stream, tags::Dict{String, String})
    url = join([BASEURL, "setstreamtags"], "/")
    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "expectedPropertyVersion" => stream.propertyVersion,
        "tags" => dict2array(tags),
        "collection" => stream.collection
    )

    response = apicall(url, JSON.json(payload))
    result = JSON.parse(response)
    checkstat(result)

    return refresh(stream)
end

"""
    setannotations(stream::Stream, annotations::Dict{String, Any})

Overwrites existing annotations.

**Notes**

Stream annotations are provided so that you may store Stream related custom
metadata.

You cannot delete annotations at this time though this will change in a
future update.
"""
function setannotations(stream::Stream, annotations::Dict{String, Any})
    url = join([BASEURL, "setstreamannotations"], "/")
    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "expectedPropertyVersion" => stream.propertyVersion,
        "changes" => dict2array(annotations),
    )

    response = apicall(url, JSON.json(payload))
    result = JSON.parse(response)
    checkstat(result)

    return refresh(stream)
end

###############################################################################
# Stream Data Retrieval Functions
###############################################################################

"""
    nearest(stream::Stream, time::Int64, version::Int=0, backward::Bool=false)

Finds the closest point in the stream to a specified time.

Return the point nearest to the specified `time` in nanoseconds since
Epoch in the stream with `version` while specifying whether to search
forward or backward in time. If `backward` is false, the returned point
will be >= `time`. If backward is true, the returned point will be <
`time`. The version of the stream used to satisfy the query is returned.

**Arguments**

* stream : The `Stream` object to search
* time : The `Int64` time (in nanoseconds since Epoch) to search near
* version : `Int64` version of the stream to use in search
* backward : `boolean` (true) to search backwards from time, else false for forward

"""
function nearest(stream::Stream, time::Int64, version::Int=0, backward::Bool=false)
    url = join([BASEURL, "nearest"], "/")

    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "time" => string(time),
        "backward" => backward,
        "versionMajor" => stream.version
    )
    response = apicall(url, JSON.json(payload))
    result = JSON.parse(response)
    checkstat(result)

    points = RawPoint(result["value"])
    return points
end

"""
    earliest(stream::Stream, version::Int=0)

Returns the first point of data in the stream.

**Arguments**

* stream : The `Stream` object to search
* version : `Int64` version of the stream to use in search
"""
earliest(stream::Stream, version::Int=0) = nearest(stream, MINIMUM_TIME, version, false)

"""
    latest(stream::Stream, version::Int=0)

Returns last point of data in the stream.

**Arguments**

* stream : The `Stream` object to search
* version : `Int64` version of the stream to use in search
"""
latest(stream::Stream, version::Int=0) = nearest(stream, MAXIMUM_TIME, version, true)

"""
    values(stream::Stream, start::Int64, stop::Int64, version::Int=0)

Read raw values from BTrDB between time [start, stop) in nanoseconds.

RawValues queries BTrDB for the raw time series data points between
`start` and `end` time, both in nanoseconds since the Epoch for the
specified stream `version`.

**Arguments**

* start : The start time (`Int64`) in nanoseconds for the range to be queried (inclusive)
* stop : The end time (`Int64`) in nanoseconds for the range to be queried (exclusive)
* version : `Int64` version of the stream to use when querying data

**Notes**

Note that the raw data points are the original values at the sensor's
native sampling rate (assuming the time series represents measurements
from a sensor). This is the lowest level of data with the finest time
granularity. In the tree data structure of BTrDB, this data is stored in
the vector nodes.
"""
function values(stream::Stream, start::Int64, stop::Int64, version::Int=0)
    url = join([BASEURL, "rawvalues"], "/")
    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "start" => string(start),
        "end" => string(stop),
        "versionMajor" => string(version)
    )
    response = apicall(url, JSON.json(payload))
    points = RawPoint.(parse_api_results(response, "values"))

    return points
end

"""
    windows(stream::Stream, start::Int64, stop::Int64, width::Int, depth::Int, version::Int=0)

Read arbitrarily-sized windows of data from BTrDB.  StatPoint objects
will be returned representing the data for each window.

**Arguments**

* start : The start time (`Int64`) in nanoseconds for the range to be queried (inclusive)
* stop : The end time (`Int64`) in nanoseconds for the range to be queried (exclusive)
* version : `Int64` version of the stream to use when querying data
* width :  the number (`Int`) of nanoseconds in each window, subject to the depth parameter.
* depth : The precision of the window duration as a power of 2 in nanoseconds. E.g 30 would make the window duration accurate to roughly 1 second

**Notes**

Windows returns arbitrary precision windows from BTrDB. It is slower
than AlignedWindows, but still significantly faster than RawValues. Each
returned window will be `width` nanoseconds long. `start` is inclusive,
but `end` is exclusive (e.g if end < start+width you will get no
results). That is, results will be returned for all windows that start
at a time less than the end timestamp. If (`end` - `start`) is not a
multiple of width, then end will be decreased to the greatest value less
than end such that (end - start) is a multiple of `width` (i.e., we set
end = start + width * floordiv(end - start, width). The `depth`
parameter is an optimization that can be used to speed up queries on
fast queries. Each window will be accurate to 2^depth nanoseconds. If
depth is zero, the results are accurate to the nanosecond. On a dense
stream for large windows, this accuracy may not be required. For example
for a window of a day, +- one second may be appropriate, so a depth of
30 can be specified. This is much faster to execute on the database
side.
"""
function windows(stream::Stream, start::Int64, stop::Int64, width::Int, depth::Int, version::Int=0)
    url = join([BASEURL, "windows"], "/")
    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "start" => string(start),
        "end" => string(stop),
        "width" => string(width),
        "depth" => depth,
        "versionMajor" => version
    )
    response = apicall(url, JSON.json(payload))
    points = StatPoint.(parse_api_results(response, "values"))

    return points
end

"""
    aligned_windows(stream::Stream, start::Int64, stop::Int64, pointwidth::Int, version::Int=0)

Read statistical aggregates of windows of data from BTrDB.

Query BTrDB for aggregates (or roll ups or windows) of the time series
with `version` between time `start` (inclusive) and `end` (exclusive) in
nanoseconds. Each point returned is a statistical aggregate of all the
raw data within a window of width 2**`pointwidth` nanoseconds. These
statistical aggregates currently include the mean, minimum, and maximum
of the data and the count of data points composing the window.

Note that `start` is inclusive, but `end` is exclusive. That is, results
will be returned for all windows that start in the interval [start, end).
If end < start+2^pointwidth you will not get any results. If start and
end are not powers of two, the bottom pointwidth bits will be cleared.
Each window will contain statistical summaries of the window.
Statistical points with count == 0 will be omitted.

**Arguments**

* start : The start time (`Int64`) in nanoseconds for the range to be queried (inclusive)
* stop : The end time (`Int64`) in nanoseconds for the range to be queried (exclusive)
* pointwidth : the number of ns between data points (2**pointwidth)
* version : `Int64` version of the stream to use when querying data

**Notes**

As the window-width is a power-of-two, it aligns with BTrDB internal
tree data structure and is faster to execute than `windows()`.
"""
function aligned_windows(stream::Stream, start::Int64, stop::Int64, pointwidth::Int, version::Int=0)
    url = join([BASEURL, "alignedwindows"], "/")
    payload = Dict(
        "uuid" => encodeUUID(stream.uuid),
        "start" => string(start),
        "end" => string(stop),
        "pointWidth" => string(pointwidth),
        "versionMajor" => version
    )
    response = apicall(url, JSON.json(payload))
    points = StatPoint.(parse_api_results(response, "values"))

    return points

end

###############################################################################
# Stream Data Management Functions
###############################################################################

"""
    insert(uuid::String, data::Array{Pair{Int64, Float64},1})

Insert new data in the form (time, value) into the series.

Inserts a list of new (time, value) tuples into the series. The tuples
in the list need not be sorted by time. If the arrays are larger than
appropriate, this function will automatically chunk the inserts. As a
consequence, the insert is not necessarily atomic, but can be used with
a very large array.

**Arguments**

* uuid : a `String` for the stream's unique identifier
* data : an `Array` of `Pair` objects representings points where the first pair item is the time (`Int64`) and the second item is the value (`Float64`)
"""
function insert(uuid::String, data::Array{Pair{Int64, Float64},1})
    objs = [Dict("time" => ii[1], "value" => ii[2]) for ii in data]

    url = join([BASEURL, "insert"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid),
        "sync" => true,
        "values" => objs
    )
    response = apicall(url, JSON.json(payload))
    result = JSON.parse(response)
    checkstat(result)

    return result["versionMajor"]
end

"""
    insert(stream::Stream, data::Array{Pair{Int64, Float64},1})

Insert new data in the form (time, value) into the series.

Inserts a list of new (time, value) tuples into the series. The tuples
in the list need not be sorted by time. If the arrays are larger than
appropriate, this function will automatically chunk the inserts. As a
consequence, the insert is not necessarily atomic, but can be used with
a very large array.

**Arguments**

* stream : the `Stream` to use for data insertion
* data : an `Array` of `Pair` objects representings points where the first pair item is the time (`Int64`) and the second item is the value (`Float64`)
"""
insert(stream::Stream, data::Array{Pair{Int64, Float64},1}) = insert(stream.uuid, data)

"""
    delete(uuid::String, start::Int64, stop::Int64)

"Delete" all points between `start` (inclusive) and `end` (exclusive),
both in nanoseconds. As BTrDB has persistent multiversioning, the
deleted points will still exist as part of an older version of the
stream.

**Arguments**

* uuid : the `String` UUID of the stream to use for data deletion
* start : the start time (`Int64`) in nanoseconds for the range to be deleted (inclusive)
* stop : the end time (`Int64`) in nanoseconds for the range to be deleted (exclusive)
"""
function delete(uuid::String, start::Int64, stop::Int64)

    url = join([BASEURL, "delete"], "/")
    payload = Dict(
        "uuid" => encodeUUID(uuid),
        "start" => string(start),
        "end" => string(stop)
    )
    response = apicall(url, JSON.json(payload))
    result = JSON.parse(response)
    checkstat(result)

    return result["versionMajor"]
end

"""
    delete(stream::Stream, start::Int64, stop::Int64)

"Delete" all points between `start` (inclusive) and `end` (exclusive),
both in nanoseconds. As BTrDB has persistent multiversioning, the
deleted points will still exist as part of an older version of the
stream.

**Arguments**

* stream : the `Stream` to use for data deletion
* start : the start time (`Int64`) in nanoseconds for the range to be deleted (inclusive)
* stop : the end time (`Int64`) in nanoseconds for the range to be deleted (exclusive)
"""
delete(stream::Stream, start::Int64, stop::Int64) = delete(stream.uuid, start, stop)

end # module


text:
using FileIO: load
import LibSndFile
using SampledSignals
using FFTW
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--within"
            help = "video/audio as input to be aligned with --with"
            required = true
        "--find-offset-of"
            help = "another video/audio file"
            required = true
    end

    return parse_args(s)
end

function main()
    # Parse the program parameters
    parsed_args = parse_commandline()
    input_filepath = parsed_args["within"]
    find_offest_input_filepath = parsed_args["find-offset-of"]

    # Convert the input files to ogg, samplerate of 16KHz
    convert_to_audio(input_filepath)
    convert_to_audio(find_offest_input_filepath)

    # Load the audio files
    audio_1 = load(audio_filepath(input_filepath))
    audio_2 = load(audio_filepath(find_offest_input_filepath))

    # Pad the signals
    audio_max_size = 2 * max(length(audio_2), length(audio_1))

    audio_1_padded = padded(audio_1, audio_max_size)
    audio_2_padded = padded(audio_2, audio_max_size)
    audio_2_padded_reversed = reverse(audio_2_padded, dims=1)

    # Perform the convolutions
    convolution = irfft(rfft(audio_1_padded) .* rfft(audio_2_padded_reversed), length(audio_1_padded))

    # Find the lag in seconds
    peak_bin, peak_value = argmax(convolution[:, 1]), maximum(convolution[:, 1])
    lag_seconds = peak_bin / audio_1.samplerate
    # The FFT performs a ciruclar convolution, so the peak_value can be on the second half of the padded audio signal
    # which would be longer than the audio file itself. Ensure the offset is lower than the length of the audio file.
    convolution_length_seconds =  length(convolution[:, 1]) / audio_1.samplerate
    lag_seconds = min(lag_seconds, convolution_length_seconds - lag_seconds)

    # Improve the lag
    p = interpolate(convolution[:, 1], peak_bin)
    lag_seconds_improved = lag_seconds + p / audio_1.samplerate

    println(lag_seconds_improved)
end

function convert_command(filepath)
    output_filepath = audio_filepath(filepath)
    return `ffmpeg -i $filepath -vn -ar 16000 -ac 1 $output_filepath -n`
end

function audio_filepath(filepath)
    # Keep the actual filename, but not the folders
    parts = rsplit(filepath, "/"; limit=2)
    return "$(tempdir())/$(parts[length(parts)]).ogg"
end

function convert_to_audio(filepath)
    if (!isfile(audio_filepath(filepath)))
        run(convert_command(filepath))
    end
end

function padded(audio, expected_size)
    pad_size = expected_size - length(audio)
    audio_pad = SampledSignals.SampleBuf(zeros(Float32, pad_size, 1), audio.samplerate)
    audio_padded = vcat(audio, audio_pad)

    return audio_padded
end

function interpolate(signal, peak_bin)
    # Quadratic Interpolation of Spectral Peaks
    # https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
    α = signal[peak_bin - 1]
    β = signal[peak_bin]
    γ = signal[peak_bin + 1]

    p = 1 / 2 * (α - γ) / (α - 2 * β + γ)
return p
end

main()


text:
module TestDistributions

# using Revise
using Test
using MLJBase
using CategoricalArrays
import Distributions
import Random.seed!
seed!(1234)

## UNIVARIATE NOMINAL

v = collect("asdfghjklzxc")
d = UnivariateFinite(v, [0.09, 0.02, 0.1, 0.1,
                          0.1, 0.1, 0.1, 0.11,
                          0.01, 0.1, 0.07, 0.1])
@test pdf(d, 's') ≈ 0.02
@test mode(d) == 'k'
@test rand(d, 5) == ['a', 'z', 'a', 'k', 'z']
@test Set(levels(d)) == Set(v)

v = collect("abcd")
d = UnivariateFinite(v, [0.2, 0.3, 0.1, 0.4])
sample = rand(d, 10^4)
freq_given_level = Distributions.countmap(sample)
pairs  = collect(freq_given_level)
sort!(pairs, by=pair->pair[2], alg=QuickSort)
sorted_levels = first.(pairs)
@test sorted_levels == ['c', 'a', 'b', 'd']

# test unseen values in pool get zero probability:
v = levels!(categorical(collect("abcd")), collect("abcdf"))
d = UnivariateFinite(v, [0.2, 0.3, 0.1, 0.4])
@test d.prob_given_level['f'] == 0
vp = broadcast(identity, v)
d = UnivariateFinite(vp, [0.2, 0.3, 0.1, 0.4])
@test d.prob_given_level['f'] == 0
vpp = levels!(categorical(["x", "y"]), ["x", "y", "z"])
d = UnivariateFinite(vpp, [0.2, 0.8])
@test d.prob_given_level["z"] == 0

v = categorical(['a', 'b', 'a', 'b', 'c', 'b', 'a', 'a'])
d = Distributions.fit(UnivariateFinite, v)
@test pdf(d, 'a') ≈ 0.5
@test pdf(d, 'b') ≈ 0.375
@test pdf(d, 'c') ≈ 0.125

# to check fitting to categorical returns zero prob for missing
# levels, we add and drop new level to a categorical version of v:
w = levels!(v, ['a', 'b', 'c', 'f'])
e = Distributions.fit(UnivariateFinite, w)
@test e.prob_given_level['f'] == 0

# arithmetic
d1 = UnivariateFinite(['a', 'b'], [0.2, 0.8])
d2 = UnivariateFinite(['b', 'c'], [0.3, 0.7])
d = average([d1, d2])
@test d.prob_given_level['a'] ≈ 0.1 && d.prob_given_level['b'] ≈ 0.55 &&  d.prob_given_level['c'] ≈ 0.35

end # module

true


text:
@testset "parser" begin

    @testset "paste event" begin

        str = "This is a good parser"
        @test T.parse_sequence(str) === T.PasteEvent(str)

    end

    @testset "parse_single_char_sequence" begin

        # @info "Enter"
        enter_event = T.KeyPressedEvent(T.ENTER)
        @test T.parse_sequence("\n") === enter_event
        @test T.parse_sequence("\r") === enter_event

        # @info "Esc"
        @test T.parse_sequence("\e") === T.KeyPressedEvent(T.ESC)

        # @info "Teb"
        @test T.parse_sequence("\t") === T.KeyPressedEvent(T.TAB)

        # @info "Backspace"
        @test T.parse_sequence("\x7F") === T.KeyPressedEvent(T.BACKSPACE)

        # @info "Null"
        @test T.parse_sequence("\0") === T.KeyPressedEvent(T.NULL)

        # @info "Ctrl 'a' - 'z', '1' - '4'"
        for i=0x01:0x1A
            if i == 'i'-'a'+0x01
                # char = '\t'
                e = T.KeyPressedEvent(T.TAB)
            elseif i == 'j'-'a'+0x01
                # char = '\n'
                e = T.KeyPressedEvent(T.ENTER)
            elseif i == 'm'-'a'+0x01
                # char = '\r'
                e = T.KeyPressedEvent(T.ENTER)
            else
                e = T.KeyPressedEvent(Char(i-0x01+'a'), T.CTRL)
            end

            @test T.parse_sequence(string(Char(i))) === e
        end
        for i=0x1C:0x1F
            char = Char(i-0x1C+'4')
            @test T.parse_sequence(string(Char(i))) === T.KeyPressedEvent(char, T.CTRL)
        end

        # @info "Char"
        @test T.parse_sequence("p") === T.KeyPressedEvent('p')
        @test T.parse_sequence("😁") === T.KeyPressedEvent('😁') # TODO: This doesn't seem right

    end

    @testset "parse_esc_leaded_sequence" begin

        @testset "F1 - F4" begin

            @test T.parse_sequence("\eOP") === T.KeyPressedEvent(T.F1)
            @test T.parse_sequence("\eOQ") === T.KeyPressedEvent(T.F2)
            @test T.parse_sequence("\eOR") === T.KeyPressedEvent(T.F3)
            @test T.parse_sequence("\eOS") === T.KeyPressedEvent(T.F4)

            @testset "F1 - F4 fallback" begin

                str = "\eOz"
                @test T.parse_sequence(str) === T.PasteEvent(str)

            end

        end

        @testset "parse_csi" begin

            @testset "parse_csi fallback" begin

                str = "\e["
                @test T.parse_sequence(str) === T.PasteEvent(str)

            end

            function create_adjoint(adjoint_code::Int)
                if adjoint_code == 2
                    return T.SHIFT
                elseif adjoint_code == 3
                    return T.CTRL
                elseif adjoint_code == 4
                    return T.SHIFT_ALT
                elseif adjoint_code == 5
                    return T.ALT
                elseif adjoint_code == 6
                    return T.SHIFT_CTRL
                elseif adjoint_code == 7
                    return T.CTRL_ALT
                elseif adjoint_code == 8
                    return T.SHIFT_CTRL_ALT
                else
                    return -1
                end
            end

            @testset "vt sequence" begin

                @testset "vt sequence fallback" begin

                    str = "\e[;~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[a;~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[;b~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[a;b~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[1;9~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)

                    str = "\e[~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[a~"
                    @test T.parse_sequence(str) === T.PasteEvent(str)

                end

                function test_vt(; adjoint_code=-1, adjoint=T.NO_ADJOINT)
                    if adjoint_code != -1
                        adjoint_str = ";$(adjoint_code)"
                    else
                        adjoint_str = ""
                    end

                    @test T.parse_sequence("\e[1$(adjoint_str)~") === T.KeyPressedEvent(T.HOME, adjoint)
                    @test T.parse_sequence("\e[2$(adjoint_str)~") === T.KeyPressedEvent(T.INSERT, adjoint)
                    @test T.parse_sequence("\e[3$(adjoint_str)~") === T.KeyPressedEvent(T.DELETE, adjoint)
                    @test T.parse_sequence("\e[4$(adjoint_str)~") === T.KeyPressedEvent(T.END, adjoint)
                    @test T.parse_sequence("\e[5$(adjoint_str)~") === T.KeyPressedEvent(T.PAGEUP, adjoint)
                    @test T.parse_sequence("\e[6$(adjoint_str)~") === T.KeyPressedEvent(T.PAGEDOWN, adjoint)
                    @test T.parse_sequence("\e[7$(adjoint_str)~") === T.KeyPressedEvent(T.HOME, adjoint)
                    @test T.parse_sequence("\e[8$(adjoint_str)~") === T.KeyPressedEvent(T.END, adjoint)
                    @test T.parse_sequence("\e[9$(adjoint_str)~") === T.PasteEvent("\e[9$(adjoint_str)~")
                    @test T.parse_sequence("\e[10$(adjoint_str)~") === T.PasteEvent("\e[10$(adjoint_str)~") # F0
                    @test T.parse_sequence("\e[11$(adjoint_str)~") === T.KeyPressedEvent(T.F1, adjoint)
                    @test T.parse_sequence("\e[12$(adjoint_str)~") === T.KeyPressedEvent(T.F2, adjoint)
                    @test T.parse_sequence("\e[13$(adjoint_str)~") === T.KeyPressedEvent(T.F3, adjoint)
                    @test T.parse_sequence("\e[14$(adjoint_str)~") === T.KeyPressedEvent(T.F4, adjoint)
                    @test T.parse_sequence("\e[15$(adjoint_str)~") === T.KeyPressedEvent(T.F5, adjoint)
                    @test T.parse_sequence("\e[16$(adjoint_str)~") === T.PasteEvent("\e[16$(adjoint_str)~")
                    @test T.parse_sequence("\e[17$(adjoint_str)~") === T.KeyPressedEvent(T.F6, adjoint)
                    @test T.parse_sequence("\e[18$(adjoint_str)~") === T.KeyPressedEvent(T.F7, adjoint)
                    @test T.parse_sequence("\e[19$(adjoint_str)~") === T.KeyPressedEvent(T.F8, adjoint)
                    @test T.parse_sequence("\e[20$(adjoint_str)~") === T.KeyPressedEvent(T.F9, adjoint)
                    @test T.parse_sequence("\e[21$(adjoint_str)~") === T.KeyPressedEvent(T.F10, adjoint)
                    @test T.parse_sequence("\e[22$(adjoint_str)~") === T.PasteEvent("\e[22$(adjoint_str)~")
                    @test T.parse_sequence("\e[23$(adjoint_str)~") === T.KeyPressedEvent(T.F11, adjoint)
                    @test T.parse_sequence("\e[24$(adjoint_str)~") === T.KeyPressedEvent(T.F12, adjoint)
                    @test T.parse_sequence("\e[25$(adjoint_str)~") === T.PasteEvent("\e[25$(adjoint_str)~") # F13
                    @test T.parse_sequence("\e[26$(adjoint_str)~") === T.PasteEvent("\e[26$(adjoint_str)~") # F14
                    @test T.parse_sequence("\e[27$(adjoint_str)~") === T.PasteEvent("\e[27$(adjoint_str)~")
                    @test T.parse_sequence("\e[28$(adjoint_str)~") === T.PasteEvent("\e[28$(adjoint_str)~") # F15
                    @test T.parse_sequence("\e[29$(adjoint_str)~") === T.PasteEvent("\e[29$(adjoint_str)~") # F16
                    @test T.parse_sequence("\e[30$(adjoint_str)~") === T.PasteEvent("\e[30$(adjoint_str)~")
                    @test T.parse_sequence("\e[31$(adjoint_str)~") === T.PasteEvent("\e[31$(adjoint_str)~") # F17
                    @test T.parse_sequence("\e[32$(adjoint_str)~") === T.PasteEvent("\e[32$(adjoint_str)~") # F18
                    @test T.parse_sequence("\e[33$(adjoint_str)~") === T.PasteEvent("\e[33$(adjoint_str)~") # F19
                    @test T.parse_sequence("\e[34$(adjoint_str)~") === T.PasteEvent("\e[34$(adjoint_str)~") # F20
                    @test T.parse_sequence("\e[35$(adjoint_str)~") === T.PasteEvent("\e[35$(adjoint_str)~")
                end

                @testset "with AdjointKey" begin

                    # +------------------------------------------------------------------+
                    # | the form of the sequence: "\e[<code>;<adjoint_code>~" and code=1:35 |
                    # +------------------------------------------------------------------+
                    for i in 2:8
                        test_vt(adjoint_code=i, adjoint=create_adjoint(i))
                    end

                end

                @testset "without AdjointKey" begin

                    # +------------------------------------------------------+
                    # | the form of the sequence: "\e[<code>~" and code=1:35 |
                    # +------------------------------------------------------+
                    test_vt()

                end

            end

            @testset "xterm sequence" begin

                @testset "vt sequence fallback" begin

                    str = "\e[;"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[a;8A"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[1;8AA"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[1;aa"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[1;9A"
                    @test T.parse_sequence(str) === T.PasteEvent(str)

                    str = "\e["
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[99"
                    @test T.parse_sequence(str) === T.PasteEvent(str)
                    str = "\e[a"
                    @test T.parse_sequence(str) === T.PasteEvent(str)

                end

                function test_xterm(; adjoint_code=-1, adjoint=T.NO_ADJOINT)
                    if adjoint_code != -1
                        adjoint_str = "1;$(adjoint_code)"
                    else
                        adjoint_str = ""
                    end

                    @test T.parse_sequence("\e[$(adjoint_str)A") === T.KeyPressedEvent(T.UP, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)B") === T.KeyPressedEvent(T.DOWN, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)C") === T.KeyPressedEvent(T.RIGHT, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)D") === T.KeyPressedEvent(T.LEFT, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)E") === T.PasteEvent("\e[$(adjoint_str)E")
                    @test T.parse_sequence("\e[$(adjoint_str)F") === T.KeyPressedEvent(T.END, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)G") === T.KeyPressedEvent('5', adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)H") === T.KeyPressedEvent(T.HOME, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)I") === T.PasteEvent("\e[$(adjoint_str)I")
                    @test T.parse_sequence("\e[$(adjoint_str)J") === T.PasteEvent("\e[$(adjoint_str)J")
                    @test T.parse_sequence("\e[$(adjoint_str)K") === T.PasteEvent("\e[$(adjoint_str)K")
                    @test T.parse_sequence("\e[$(adjoint_str)L") === T.PasteEvent("\e[$(adjoint_str)L")
                    @test T.parse_sequence("\e[$(adjoint_str)M") === T.PasteEvent("\e[$(adjoint_str)M")
                    @test T.parse_sequence("\e[$(adjoint_str)N") === T.PasteEvent("\e[$(adjoint_str)N")
                    @test T.parse_sequence("\e[$(adjoint_str)O") === T.PasteEvent("\e[$(adjoint_str)O")
                    if adjoint_code != -1
                    @test T.parse_sequence("\e[$(adjoint_str)P") === T.KeyPressedEvent(T.F1, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)Q") === T.KeyPressedEvent(T.F2, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)R") === T.KeyPressedEvent(T.F3, adjoint)
                    @test T.parse_sequence("\e[$(adjoint_str)S") === T.KeyPressedEvent(T.F4, adjoint)
                    end
                    @test T.parse_sequence("\e[$(adjoint_str)T") === T.PasteEvent("\e[$(adjoint_str)T")
                    @test T.parse_sequence("\e[$(adjoint_str)U") === T.PasteEvent("\e[$(adjoint_str)U")
                    @test T.parse_sequence("\e[$(adjoint_str)V") === T.PasteEvent("\e[$(adjoint_str)V")
                    @test T.parse_sequence("\e[$(adjoint_str)W") === T.PasteEvent("\e[$(adjoint_str)W")
                    @test T.parse_sequence("\e[$(adjoint_str)X") === T.PasteEvent("\e[$(adjoint_str)X")
                    @test T.parse_sequence("\e[$(adjoint_str)Y") === T.PasteEvent("\e[$(adjoint_str)Y")
                    @test T.parse_sequence("\e[$(adjoint_str)Z") === T.KeyPressedEvent(T.BACKTAB, adjoint)
                end

                @testset "with AdjointKey" begin

                    # +-----------------------------------------------------------------+
                    # | the form of the sequence: "\e[1;<adjoint_code><code>" and code=A:Z |
                    # +-----------------------------------------------------------------+
                    for i in 2:8
                        test_xterm(adjoint_code=i, adjoint=create_adjoint(i))
                    end

                end

                @testset "without AdjointKey" begin

                    # +----------------------------------------------------+
                    # | the form of the sequence: "\e[<code>" and code=A:Z |
                    # +----------------------------------------------------+
                    test_xterm()

                end

            end

        end

        @testset "ALT + Char" begin

            @test T.parse_sequence("\eA") === T.KeyPressedEvent('A', T.ALT)

        end

        @testset "parse_esc_leaded_sequence fallback" begin

            str = "\efallback"
            @test T.parse_sequence(str) === T.PasteEvent(str)

        end

    end

end