text:
using Markdown # hide

# ## All-KNN
# By Eric S. Téllez
#

ENV["GKSwstype"] = "nul" # hide

using SimilaritySearch
# Computing the $k$ nearest neighbors of a dataset (all vs all) is a useful task to take knowledge
# of a given dataset. This is a core task for some clustering algorithms and non-linear dimensional reduction
# algorithms.
#
# Given a metric database $(X, dist)$  and a relatively small $k$ value, the goal is
# to compute $\{ knn(x) \mid x \in X \}$ taking into account that each $x_i \in X$, and therefore, $x_i$ should be
# removed from the $i$-th $knn$ result set.
#
# Solving `allknn` fast and accuratelly is the goal of this example.

const dim = 16
k = 15
db = MatrixDatabase(randn(Float32, dim, 10^5))
dist = SqL2Distance()
G = SearchGraph(; dist, db)
opt = OptimizeParameters(kind=MinRecall(0.9))
itime = @elapsed (index!(G); optimize!(G, opt))
nothing # hide
# Now we can solve all $k$nn

allknntime = @elapsed knns, dists = allknn(G, k; parallel=true)
nothing # hide
# ## Differences between `allknn(G, k)` and `searchbatch(G, X, k)`
#
# We can solve similarly with `searchbatch` but self-references should be removed later, and more important,
# `allknn` use special pivoting/boosting strategies that yields to faster searches.

searchtime = @elapsed sknns, sdists = searchbatch(G, db, k; parallel=true)
nothing # hide

# ## About the cost of construction + `allknn` instead of a brute force computation.
# `allknn` for `ExhaustiveSearch` doesn't perform any optimization but removes self references.
etime = @elapsed eknns, edists = allknn(ExhaustiveSearch(; db, dist), k; parallel=true)
nothing # hide

# ## Comparing solution times
#
# indexing, allknn, and searchbatch times
itime, allknntime, searchtime

# full cost `allknn`
itime + allknntime
# full cost `searchbatch`
itime + searchtime
# exhaustive `allknn`
etime

# ## Quality
# macro recall of `allknn`
macrorecall(eknns, knns)
# macro recall of `searchbatch`
macrorecall(eknns, sknns)

# ## Final notes:
# Exhaustive search will fetch the exact solution but it has a higher cost and this could be more
# notorious as dataset's size increases.
# Also note that even when we pass `parallel=true` it runs in a single thread due to the html generation pipeline.
#

#- Using BeutifulMakie style
# ## Dependencies
using Pkg # hide
Pkg.status(["SimilaritySearch"]) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
import Base

"""Base class for a single optimization objective.

All objectives must have a field `initial_state` and a field `generator`, at
minimum.
"""
abstract type AbstractControlObjective end

"""Standard optimization objective.

```julia
Objective(;
    initial_state=<initial_state>,
    generator=<genenerator>,
    target_state=<target_state>
)
```

describes an optimization objective where the time evoluation of the given
`initial_state` under the given `generator` aims towards `target_state`. The
`generator` here is e.g. a time-dependent Hamiltonian or Liouvillian.

The most common control problems in quantum control, e.g. state-to-state
transitions or quantum gate implementations can be expressed by simultaneously
fulfilling multiple objectives of this type.

Note that the objective can only be instantiated via keyword arguments.
"""
struct Objective{ST, GT} <: AbstractControlObjective
    initial_state :: ST
    generator :: GT
    target_state :: ST
    function Objective(;initial_state::ST, generator::GT, target_state::ST) where {ST, GT}
        new{ST, GT}(initial_state, generator, target_state)
    end
end

"""Standard optimization objective with a weight.

```julia
WeightedObjective(;
    initial_state=<initial_state>,
    generator=<genenerator>,
    target_state=<target_state>,
    weight=<weight>
)
```

initializes a control objective like [`Objective`](@ref), but with an
additional `weight` parameter (a float generally between 0 and 1) that weights
the objective relative to other objectives that are part of the same control
problem.
"""
struct WeightedObjective{ST, GT} <: AbstractControlObjective
    initial_state :: ST
    generator :: GT
    target_state :: ST
    weight :: Float64
    function WeightedObjective(;initial_state::ST, generator::GT, target_state::ST, weight::Float64) where {ST, GT}
        new{ST, GT}(initial_state, generator, target_state, weight)
    end
end

"""A full control problem with multiple objectives.

```julia
ControlProblem(
   objectives=<list of objectives>,
   pulse_options=<dict of controls to pulse options>,
   tlist=<time grid>,
   kwargs...
)
```

Note that the control problem can only be instantiated via keyword arguments.

The `objectives` are a list of [`AbstractControlObjective`](@ref) instances,
each defining an initial state and a dynamical generator for the evolution of
that state. Usually, the objective will also include a target state (see
[`Objective`](@ref)) and possibly a weight (see [`WeightedObjective`](@ref)).

The `pulse_options` are a dictionary (`IdDict`) mapping controls that occur in
the `objectives` to properties specific to the control method.

The `tlist` is the time grid on which the time evolution of the initial states
of each objective should be propagated.

The remaining `kwargs` are keyword arguments that are passed directly to the
optimal control method. These typically include e.g. the optimization
functional.

The control problem is solved by finding a set of controls that simultaneously
fulfill all objectives.
"""
struct ControlProblem{OT<:AbstractControlObjective, OPT<:AbstractDict, KWT}
    objectives :: Vector{OT}
    pulse_options :: OPT # TODO: make this optional?
    tlist :: Vector{Float64}
    kwargs :: KWT
    function ControlProblem(;objectives, pulse_options, tlist, kwargs...)
        kwargs_dict = Dict(kwargs)  # make the kwargs mutable
        new{typeof(objectives[1]), typeof(pulse_options), typeof(kwargs_dict)}(
            objectives, pulse_options, tlist, kwargs_dict
        )
    end
end

function Base.copy(problem::ControlProblem)
    return ControlProblem(
        objectives=problem.objectives,
        pulse_options=problem.pulse_options,
        tlist=problem.tlist;
        problem.kwargs...
    )
end

# adjoint for the nested-tuple dynamical generator (e.g. `(H0, (H1, ϵ))`)
function dynamical_generator_adjoint(G::Tuple)
    result = []
    for part in G
        if isa(part, Tuple)
            push!(result, (Base.adjoint(part[1]), part[2]))
        else
            push!(result, Base.adjoint(part))
        end
    end
    return Tuple(result)
end

# fallback adjoint
dynamical_generator_adjoint(G) = Base.adjoint(G)

"""
```julia
adjoint(objective)
```

Adjoint of a control objective. The adjoint objective contains the adjoint of
the dynamical generator `obj.generator`, and adjoints of the
`obj.initial_state` / `obj.target_state` if these exist and have an adjoint.
"""
function Base.adjoint(obj::Objective)
    initial_state_adj = obj.initial_state
    try initial_state_adj = Base.adjoint(obj.initial_state) catch end
    generator_adj = dynamical_generator_adjoint(obj.generator)
    target_adj = obj.target_state
    try target_adj = Base.adjoint(obj.target_state) catch end
    return Objective(
        initial_state=initial_state_adj,
        generator=generator_adj,
        target_state=target_adj
    )
end

function Base.adjoint(obj::WeightedObjective)
    initial_state_adj = obj.initial_state
    try initial_state_adj = Base.adjoint(obj.initial_state) catch end
    generator_adj = dynamical_generator_adjoint(obj.generator)
    target_adj = obj.target_state
    try target_adj = Base.adjoint(obj.target_state) catch end
    return WeightedObjective(
        initial_state=initial_state_adj,
        generator=generator_adj,
        target_state=target_adj,
        weight=obj.weight
    )
end

ext: jl


text:
using ITensors,
      Test

include("util.jl")

@testset "MPO Basics" begin
  N = 6
  sites = [Index(2,"Site") for n=1:N]
  @test length(MPO()) == 0
  O = MPO(sites)
  @test length(O) == N

  str = split(sprint(show, O), '\n')
  @test str[1] == "MPO"
  @test length(str) == length(O) + 2

  O[1] = ITensor(sites[1], prime(sites[1]))
  @test hasindex(O[1],sites[1])
  @test hasindex(O[1],prime(sites[1]))
  P = copy(O)
  @test hasindex(P[1],sites[1])
  @test hasindex(P[1],prime(sites[1]))

  @testset "orthogonalize" begin
    phi = randomMPS(sites)
    K = randomMPO(sites)
    orthogonalize!(phi, 1)
    orthogonalize!(K, 1)
    orig_inner = inner(phi, K, phi)
    orthogonalize!(phi, div(N, 2))
    orthogonalize!(K, div(N, 2))
    @test inner(phi, K, phi) ≈ orig_inner
  end

  @testset "inner <y|A|x>" begin
    phi = randomMPS(sites)
    K = randomMPO(sites)
    @test maxLinkDim(K) == 1
    psi = randomMPS(sites)
    phidag = dag(phi)
    prime!(phidag)
    phiKpsi = phidag[1]*K[1]*psi[1]
    for j = 2:N
      phiKpsi *= phidag[j]*K[j]*psi[j]
    end
    @test phiKpsi[] ≈ inner(phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi,K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
        mpo_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors2 = ITensor[ITensor() for ii in 1:N]
        mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mpo_tensors[1] = randomITensor(mpo_link_inds[1], sites[1], sites[1]')
        mps_tensors[1] = randomITensor(mps_link_inds[1], sites[1])
        mps_tensors2[1] = randomITensor(mps_link_inds[1], sites[1])
        for ii in 2:N-1
            mpo_tensors[ii] = randomITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
            mps_tensors[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
            mps_tensors2[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        end
        mpo_tensors[N] = randomITensor(mpo_link_inds[N-1], sites[N], sites[N]')
        mps_tensors[N] = randomITensor(mps_link_inds[N-1], sites[N])
        mps_tensors2[N] = randomITensor(mps_link_inds[N-1], sites[N])
        K   = MPO(N, mpo_tensors, 0, N+1)
        psi = MPS(N, mps_tensors, 0, N+1)
        phi = MPS(N, mps_tensors2, 0, N+1)
        orthogonalize!(psi, 1; maxdim=link_dim)
        orthogonalize!(K, 1; maxdim=link_dim)
        orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
        phidag = dag(phi)
        prime!(phidag)
        phiKpsi = phidag[1]*K[1]*psi[1]
        for j = 2:N
          phiKpsi *= phidag[j]*K[j]*psi[j]
        end
        @test scalar(phiKpsi) ≈ inner(phi,K,psi)
    end
  end

  @testset "inner <By|A|x>" begin
    phi = makeRandomMPS(sites)

    K = makeRandomMPO(sites,chi=2)
    J = makeRandomMPO(sites,chi=2)

    psi = makeRandomMPS(sites)
    phidag = dag(phi)
    prime!(phidag, 2)
    Jdag = dag(J)
    prime!(Jdag)
    for j ∈ eachindex(Jdag)
      swapprime!(inds(Jdag[j]),2,3)
      swapprime!(inds(Jdag[j]),1,2)
      swapprime!(inds(Jdag[j]),3,1)
    end

    phiJdagKpsi = phidag[1]*Jdag[1]*K[1]*psi[1]
    for j ∈ eachindex(psi)[2:end]
      phiJdagKpsi = phiJdagKpsi*phidag[j]*Jdag[j]*K[j]*psi[j]
    end

    @test scalar(phiJdagKpsi) ≈ inner(J,phi,K,psi)

    ## Do contraction manually.
    #O = 1.
    #for j ∈ eachindex(phi)
    #    psij = reshape(array(psi[j]),2)
    #    phij = reshape(array(phi[j]),2)
    #    Kj = reshape(array(K[j]),2,2)
    #    Jj = reshape(array(J[j]),2,2)
    #    O *= ((transpose(Jj)*phij)'*transpose(Kj)*psij)[]
    #end
    #@test O ≈ inner(J,phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(J,phi,K,badpsi)
  end

  @testset "errorMPOProd" begin
    phi = makeRandomMPS(sites)
    K = makeRandomMPO(sites,chi=2)

    psi = makeRandomMPS(sites)

    dist = sqrt(abs(1 + (inner(phi,phi) - 2*real(inner(phi,K,psi)))
                        /inner(K,psi,K,psi)))
    @test dist ≈ errorMPOProd(phi,K,psi)

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    # Apply K to phi and check that errorMPOProd is close to 0.
    Kphi = applyMPO(K,phi;method="naive", cutoff=1E-8)
    @test errorMPOProd(Kphi, K, phi) ≈ 0. atol=1e-4

    @test_throws DimensionMismatch applyMPO(K,badpsi;method="naive", cutoff=1E-8)
    @test_throws DimensionMismatch errorMPOProd(phi,K,badpsi)
  end

  @testset "applyMPO" begin
    phi = randomMPS(sites)
    K   = randomMPO(sites)
    @test maxLinkDim(K) == 1
    psi = randomMPS(sites)
    psi_out = applyMPO(K, psi,maxdim=1)
    @test inner(phi,psi_out) ≈ inner(phi,K,psi)
    @test_throws ArgumentError applyMPO(K, psi, method="fakemethod")

    badsites = [Index(2,"Site") for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch applyMPO(K,badpsi)

    # make bigger random MPO...
    for link_dim in 2:5
        mpo_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors  = ITensor[ITensor() for ii in 1:N]
        mps_tensors2 = ITensor[ITensor() for ii in 1:N]
        mpo_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mps_link_inds = [Index(link_dim, "r$ii,Link") for ii in 1:N-1]
        mpo_tensors[1] = randomITensor(mpo_link_inds[1], sites[1], sites[1]')
        mps_tensors[1] = randomITensor(mps_link_inds[1], sites[1])
        mps_tensors2[1] = randomITensor(mps_link_inds[1], sites[1])
        for ii in 2:N-1
            mpo_tensors[ii] = randomITensor(mpo_link_inds[ii], mpo_link_inds[ii-1], sites[ii], sites[ii]')
            mps_tensors[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
            mps_tensors2[ii] = randomITensor(mps_link_inds[ii], mps_link_inds[ii-1], sites[ii])
        end
        mpo_tensors[N] = randomITensor(mpo_link_inds[N-1], sites[N], sites[N]')
        mps_tensors[N] = randomITensor(mps_link_inds[N-1], sites[N])
        mps_tensors2[N] = randomITensor(mps_link_inds[N-1], sites[N])
        K   = MPO(N, mpo_tensors, 0, N+1)
        psi = MPS(N, mps_tensors, 0, N+1)
        phi = MPS(N, mps_tensors2, 0, N+1)
        orthogonalize!(psi, 1; maxdim=link_dim)
        orthogonalize!(K, 1; maxdim=link_dim)
        orthogonalize!(phi, 1; normalize=true, maxdim=link_dim)
        psi_out = applyMPO(deepcopy(K), deepcopy(psi); maxdim=10*link_dim, cutoff=0.0)
        @test inner(phi, psi_out) ≈ inner(phi, K, psi)
    end
  end
  @testset "add" begin
    shsites = spinHalfSites(N)
    K = randomMPO(shsites)
    L = randomMPO(shsites)
    M = sum(K, L)
    @test length(M) == N
    psi = randomMPS(shsites)
    k_psi = applyMPO(K, psi, maxdim=1)
    l_psi = applyMPO(L, psi, maxdim=1)
    @test inner(psi, sum(k_psi, l_psi)) ≈ inner(psi, M, psi) atol=5e-3
  end

  @testset "nmultMPO" begin
    psi = randomMPS(sites)
    K = randomMPO(sites)
    L = randomMPO(sites)
    @test maxLinkDim(K) == 1
    @test maxLinkDim(L) == 1
    KL = nmultMPO(K, L, maxdim=1)
    psi_kl_out = applyMPO(K, applyMPO(L, psi, maxdim=1), maxdim=1)
    @test inner(psi,KL,psi) ≈ inner(psi, psi_kl_out) atol=5e-3

    badsites = [Index(2,"Site") for n=1:N+1]
    badL = randomMPO(badsites)
    @test_throws DimensionMismatch nmultMPO(K,badL)
  end

  sites = spinHalfSites(N)
  O = MPO(sites,"Sz")
  @test length(O) == N # just make sure this works

  @test_throws ArgumentError randomMPO(sites, 2)
  @test_throws ErrorException linkindex(MPO(N, fill(ITensor(), N), 0, N + 1), 1)
end


text:
globalUnique = 0
using DataStructures

function compile_method_instance(method_instance::Core.MethodInstance, world)
    function postprocess(ir)
        # get rid of jfptr wrappers
        for llvmf in functions(ir)
            startswith(LLVM.name(llvmf), "jfptr_") && unsafe_delete!(ir, llvmf)
        end

        return
    end

    # set-up the compiler interface
    last_method_instance = nothing
    call_stack = Vector{Core.MethodInstance}()
    dependencies = MultiDict{Core.MethodInstance,LLVM.Function}()

    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        ir = LLVM.Module(ref)
        postprocess(ir)

        # find the function that this module defines
        llvmfs = filter(llvmf -> !isdeclaration(llvmf) &&
                                 linkage(llvmf) == LLVM.API.LLVMExternalLinkage,
                        collect(functions(ir)))

        llvmf = nothing
        if length(llvmfs) == 1
            llvmf = first(llvmfs)
        elseif length(llvmfs) > 1
            llvmfs = filter!(llvmf -> startswith(LLVM.name(llvmf), "julia_"), llvmfs)
            if length(llvmfs) == 1
                llvmf = first(llvmfs)
            end
        end

        @assert llvmf !== nothing

        insert!(dependencies, last_method_instance, llvmf)
    end
    function hook_emit_function(method_instance, code, world)
        push!(call_stack, method_instance)

        # check for recursion
        if method_instance in call_stack[1:end-1]
            @error "call_stack" call_stack
            error("recursion is currently not supported")
        end
    end
    function hook_emitted_function(method, code, world)
        @assert last(call_stack) == method
        last_method_instance = pop!(call_stack)
    end
    params = Base.CodegenParams(cached             = false,
                                track_allocations  = false,
                                code_coverage      = false,
                                static_alloc       = false,
                                prefer_specsig     = true,
                                module_activation  = hook_module_activation,
                                emit_function      = hook_emit_function,
                                emitted_function   = hook_emitted_function)

    # get the code
    ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                (Any, UInt, Bool, Bool, Base.CodegenParams),
                method_instance, world, #=wrapper=#false, #=optimize=#false, params)
    if ref == C_NULL
        throw(InternalCompilerError(job, "the Julia compiler could not generate LLVM IR"))
    end
    llvmf = LLVM.Function(ref)
    ir = LLVM.parent(llvmf)
    postprocess(ir)

    return llvmf, dependencies
end

function irgen(method_instance::Core.MethodInstance, world)
    entry, dependencies = compile_method_instance(method_instance, world)
    mod = LLVM.parent(entry)

    begin
        # we disable Julia's compilation cache not to poison it with GPU-specific code.
        # as a result, we might get multiple modules for a single method instance.
        cache = Dict{String,String}()

        for called_method_instance in keys(dependencies)
            llvmfs = dependencies[called_method_instance]

            # link the first module
            llvmf = popfirst!(llvmfs)
            llvmfn = LLVM.name(llvmf)
            link!(mod, LLVM.parent(llvmf))

            # process subsequent duplicate modules
            for dup_llvmf in llvmfs
                # don't link them, but note the called function name in a cache
                dup_llvmfn = LLVM.name(dup_llvmf)
                cache[dup_llvmfn] = llvmfn
            end
        end

        # resolve function declarations with cached entries
        for llvmf in filter(isdeclaration, collect(functions(mod)))
            llvmfn = LLVM.name(llvmf)
            if haskey(cache, llvmfn)
                def_llvmfn = cache[llvmfn]
                replace_uses!(llvmf, functions(mod)[def_llvmfn])

                @assert isempty(uses(llvmf)) job
                unsafe_delete!(LLVM.parent(llvmf), llvmf)
            end
        end
    end

    llvmfn = replace(LLVM.name(entry), r"_\d+$"=>"")
    ## append a global unique counter
    global globalUnique
    globalUnique += 1
    llvmfn *= "_$globalUnique"
    LLVM.name!(entry, llvmfn)

    # minimal required optimization
    ModulePassManager() do pm
        linkage!(entry, LLVM.API.LLVMExternalLinkage)
        internalize!(pm, [LLVM.name(entry)])

        always_inliner!(pm)
        run!(pm, mod)
    end

    return mod, entry
end

function get_methodinstance(@nospecialize(sig))
    world = typemax(UInt)
    m = ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), sig, world)
    meth = m.func::Method

    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig, meth.sig)::Core.SimpleVector
    if VERSION >= v"1.2.0-DEV.320"
        meth = Base.func_for_method_checked(meth, ti, env)
    else
        meth = Base.func_for_method_checked(meth, ti)
    end

    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)
    return linfo::Core.MethodInstance, world
end

function lookup_sig(@nospecialize(func), @nospecialize(tt))
    isa(func, Core.Builtin) && error("function is not a generic function")
    return Base.signature_type(func, tt)::Type
end

function irgen(@nospecialize(func), @nospecialize(tt))
    sig = lookup_sig(func, tt)
    irgen(sig)
end

function irgen(@nospecialize(sig))
    mi, world = get_methodinstance(sig)
    return irgen(mi, world)
end


text:
using Plots
using Random
using DifferentialEquations
using LinearAlgebra
using DynamicalSystems
using Parameters
using BasicReservoirComputing

include("utilities.jl")

@inline @inbounds function colpitts_eom(x, p, t)
    α = p[1]; Γ = p[2]; q = p[3]; η = p[4]
    du1 = α*x[2]
    du2 = -Γ*(x[1] + x[3]) - q*x[2]
    du3 = η*(x[2] + 1 - exp(-x[1]))
    return SVector{3}(du1, du2, du3)
end

@inline @inbounds function colpitts_jac(x, p, t)
    α = p[1]; Γ = p[2]; q = p[3]; η = p[4]
    J = @SMatrix [0  α  0;
                 -Γ -q -Γ;
                 η*exp(-x[1]) η 0]
    return J
end

@inline @inbounds function CL63_eom(x, p, t)
    xe,ye,ze, xt,yt,zt, xo,yo,zo = x
    sigma, rho, beta, S, k1, k2, tau, c, ce, cz = p

    dxedt = sigma*(ye - xe) - ce*(S*xt + k1)
    dyedt = rho*xe - ye - xe*ze + ce*(S*yt + k1)
    dzedt = xe*ye - beta*ze

    # Tropical atmosphere system
    dxtdt = sigma*(yt - xt) - c*(S*xo + k2) - ce*(S*xe + k1)
    dytdt = rho*xt - yt - xt*zt + c*(S*yo + k2) + ce*(S*ye + k1)
    dztdt = xt*yt - beta*zt + cz*zo

    # Tropical ocean system
    dxodt = tau*sigma*(yo - xo) - c*(xt + k2)
    dyodt = tau*rho*xo - tau*yo - tau*S*xo*zo + c*(yt + k2)
    dzodt = tau*S*xo*yo - tau*beta*zo - cz*zt
    return SVector{9}(dxedt,dyedt,dzedt,
                      dxtdt,dytdt,dztdt,
                      dxodt,dyodt,dzodt)
end

@inline @inbounds function CL63_jac(x, p, t)
    xe,ye,ze, xt,yt,zt, xo,yo,zo = x
    σ, ρ, β, S, k1, k2, τ, c, ce, cz = p
    J = @SMatrix [-σ σ 0 -ce*S 0 0 0 0 0;
                  ρ-ze -1 -xe 0 ce*S 0 0 0 0;
                  ye xe -β 0 0 0 0 0 0;
                  -ce*S 0 0 -σ σ 0 -c*S 0 0;
                  0 ce*S 0 ρ-zt -1 -xt 0 c*S 0;
                  0 0 0 yt xt -β 0 0 cz;
                  0 0 0 -c 0 0 -τ*σ τ*σ 0;
                  0 0 0 0 c 0 τ*(σ-S*zo) -τ -τ*S*xo;
                  0 0 0 0 0 -cz τ*S*yo τ*S*xo -τ*β]
    return J
end

function test_rc(system::DynamicalSystem, params::NamedTuple;
                 Δt = 0.01, train_time=100, test_time = 50, n_test = 100)
    D=length(system.u0)
    utrain, utest_arr = make_data(system; train_time = train_time,
                                  test_time = test_time, n_test = n_test)

    @unpack N, SR, ρA, α, σ, β, σb, nspin = params

    # Build and Train RC
    reservoir = rc(D, N, ρA, Δt, SR, α, σ, σb, β)
    train_RC(reservoir, utrain, spinup=nspin)

    # Forecast RC
    utrue = utest_arr[1][:, nspin:end]
    nsteps = size(utrue)[2]
    upred = forecast_RC(reservoir, nsteps, uspin = utest_arr[1][:, 1:nspin])
    t = collect(0:nsteps-1)*Δt

    p1 = plot_prediction(t, upred, utrue)

    # Test on many points
    vpt = test_RC(reservoir, utest_arr; ϵ=0.3, spinup=nspin)
    p2 = histogram(vpt, legend=false,
                   overwrite_figure=false, bins = 15)
    xlabel!(p2, "Valid Prediction Time")
    pfinal = plot(p1, p2, layout = (1, 2), show=true)
    display(pfinal)
    return reservoir
end

CL63 = (σ = 0.007, α = 1.00, SR = 0.01, β = 2.201375e-08, σb = 0.66,
        N = 2000, ρA = 0.02, nspin = 200)
L63 = (σ = 0.084,α = 0.60, SR = 0.80, β = 8.493901e-08, σb = 1.60,
       N = 800, ρA = 0.02, nspin = 200)
L9610d = (σ = 0.005, α = 0.72, SR = 0.21, β = 7.640822e-09, σb = 1.47,
          N = 2000, ρA = 0.02, nspin = 200)
L965d = (σ = 0.060, α = 0.70, SR = 0.58, β = 6.332524e-09, σb = 1.59,
         N = 800, ρA = 0.02, nspin = 200)
colpitts = (σ = 0.100, α = 1.00, SR = 1.20, β = 1.000000e-08, σb = 2.00,
            N = 800, ρA = 0.02, nspin = 200)
rossler = (σ = 0.066, α = 0.47, SR = 0.50, β = 2.101845e-09, σb = 1.23,
           N = 800, ρA = 0.02, nspin = 200)

sys_L63 = Systems.lorenz()
sys_L965d = Systems.lorenz96(5, F=8.0)
sys_L9610d = Systems.lorenz96(10, F=8.0)
sys_rossler = Systems.roessler()
sys_colpitts = ContinuousDynamicalSystem(colpitts_eom,
                                         rand(3),
                                         [5.0, 0.08, 0.7, 6.3],
                                         colpitts_jac)
sys_CL63 = ContinuousDynamicalSystem(CL63_eom,
                                     [-3.1, -3.1, 20.7, -3.1, -3.1, 20.7, -3.1, -3.1, 20.7],
                                     [10, 28, 8/3., 1, 10, -11, 0.1, 1, 0.08, 1],
                                     CL63_jac)
# change arguments to test different systems
trained_reservoir = test_rc(sys_L63, L63, test_time = 25);

max_stars_repo_name: japlatt/BasicReservoirComputing

max_issues_repo_name: japlatt/BasicReservoirComputing

max_forks_repo_name: japlatt/BasicReservoirComputing


text:
using EventSimulation
using Documenter

makedocs(sitename="EventSimulation.jl")

deploydocs(repo = "github.com/bkamins/EventSimulation.jl.git",
           target = "build", deps = nothing, make = nothing)


text:
# Second iteration of explicitly exloring the marginalization of discrete variables onto the continuous state space.
# Although this code is still excessive and messy, this is a significant feature expansion; and future versions will
# generalize the marginalization process to allow for implicit hypothesis exploration.  The messy explicit version of code is
# intended to help develop the required generalistic unit tests.  The unit tests will then be validated, frozen and used to
# confirm that future "algebraic" marginalization (implicit) versions operate correctly.
# FYI, the complexity of general multihypothesis convolutions can be deceiving, however, note that the coding
# complexity is contained for each indivual factor at a time.  Global Bayes tree inference then creates the symphony
# of non-Gaussian (multimodal) posterior beliefs from the entire factor graph.
#
# 2018/6/01 @dehann

"""
    $SIGNATURES

Return common vectors `(allmhp, certainidx,uncertnidx)` used for dealing with multihypo cases.
"""
function getHypothesesVectors(mhp::Vector{Float64})::Tuple{Vector{Int},Vector{Int},Vector{Int}}
  allmhp = 1:length(mhp)
  certainidx = allmhp[mhp .== 0.0]  # TODO remove after gwp removed
  uncertnidx = allmhp[0.0 .< mhp]
  return (allmhp,certainidx,uncertnidx)
end

"""
    $(SIGNATURES)

This function explicitly codes that marginalization of a discrete categorical selection variable for ambiguous data association situations.  Improved implementations should implicitly induce the same behaviour through summation (integration) when marginalizing any number of discrete variables.  This function populates `allelements` with particle indices associated with particular multihypothesis selection while `activehypo` simultaneously contains the hypothesis index and factor graph variables associated with that hypothesis selection.  The return value `certainidx` are the hypotheses that are not in question.

Notes:
- Issue 427, race condition during initialization since n-ary variables not resolvable without other init.

```
# `allelements` example BearingRange [:x1, 0.5:l1a, 0.5:l1b]
# sfidx = (1=:x1,2=:l1a,3=:l1b)
if solvefor :x1, then allelem = [mhidx.==:l1a; mhidx.==l1b]
if solvefor :l1a, then allelem = [mhidx.==:l1a] and ARR[solvefor][:,mhidx.==:l1b]=ARR[:l1b][:,mhidx.==:l1b]
if solvefor :l1b, then allelem = [mhidx.==:l1b] and ARR[solvefor][:,mhidx.==:l1a]=ARR[:l1a][:,mhidx.==:l1a]
if solvefor 1, then allelem = [mhidx.==2; mhidx.==3]
if solvefor 2, then allelem = [mhidx.==2] and ARR[solvefor][:,mhidx.==3]=ARR[3][:,mhidx.==3]
if solvefor 3, then allelem = [mhidx.==3] and ARR[solvefor][:,mhidx.==2]=ARR[2][:,mhidx.==2]

# `activehypo` in example mh=[1.0;0.5;0.5]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3
sfidx=3, mhidx=2:  3 should take a value from 2
sfidx=3, mhidx=3:  ah = [1;3]

# `activehypo` in example mh=[1.0;0.33;0.33;0.34]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=1, mhidx=4:  ah = [1;4]

sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3
sfidx=2, mhidx=4:  2 should take a value from 4

sfidx=3, mhidx=2:  3 should take a value from 2
sfidx=3, mhidx=3:  ah = [1;3]
sfidx=3, mhidx=4:  3 should take a value from 4

sfidx=4, mhidx=2:  4 should take a value from 2
sfidx=4, mhidx=3:  4 should take a value from 3
sfidx=4, mhidx=4:  ah = [1;4]
```
"""
function assembleHypothesesElements!(mh::Categorical,
                                     maxlen::Int,
                                     sfidx::Int,
                                     lenXi::Int,
                                     isinit::Vector{Bool}=ones(Bool, lenXi)  )
  #
  allelements = []
  activehypo = []
  mhidx = Int[]

  allidx = 1:maxlen
  # @show sfidx
  allmhp, certainidx, uncertnidx = getHypothesesVectors(mh.p)
  # allmhp = 1:length(mh.p)
  # certainidx = allmhp[mh.p .== 0.0]  # TODO remove after gwp removed
  # uncertnidx = allmhp[0.0 .< mh.p]

  # select only hypotheses that can be used (ie variables have been initialized)
  @assert !(sum(isinit) == 0 && sfidx == certainidx) # cannot init from nothing for any hypothesis

  mhh = if sum(isinit) < lenXi - 1
    @assert isLeastOneHypoAvailable(sfidx, certainidx, uncertnidx, isinit)
    @info "not all hypotheses initialized, but at least one available -- see #427"
    mhp = deepcopy(mh.p)
    suppressmask = isinit .== false
    suppressmask[sfidx] = false
    mhp[suppressmask] .= 0.0
    mhp ./= sum(mhp)
    Categorical(mhp)
  else
    mh
  end

  # prep mmultihypothesis selection values
  mhidx = rand(mhh, maxlen) # selection of which hypothesis is correct

  pidx = 0
  sfincer = sfidx in certainidx
  for pval in mh.p
    pidx += 1
    pidxincer = pidx in certainidx # ??
    # permutation vectors for later computation
    iterarr = allidx[mhidx .== pidx]
    iterah = Int[]
    if !pidxincer && sfincer # 1e-15 <= pval && mh.p[sfidx] < 1e-10  # proxy for sfidx in certainidx
      # solve for one of the certain variables containing uncertain hypotheses in others
      iterah = sort(union(certainidx, pidx)) # sort([sfidx;pidx])
      # DONE -- supports n-ary factors in multihypo mode
    elseif pidxincer && !sfincer || sfidx == pidx # pval < 1e-15 && mh.p[sfidx] >= 1e-10
      # solve for one of the uncertain variables
      iterah = sort(union(certainidx, sfidx)) # sort([sfidx;pidx])
      # EXPERIMENTAL -- support more than binary factors in multihypo mode
    elseif pidxincer && sfincer # pval < 1e-15 && mh.p[sfidx] < 1e-10
      iterarr = Int[]
      iterah = Int[] # may be moot anyway, but double check first
    elseif !pidxincer && !sfincer # pval >= 1e-15 && mh.p[sfidx] >= 1e-10
      iterah = uncertnidx #allmhp[mh.p .> 1e-15]
    else
      error("Unknown hypothesis case, got sfidx=$(sfidx) with mh.p=$(mh.p), pidx=$(pidx)")
    end
    push!(allelements, iterarr)
    push!(activehypo, (pidx,iterah))
  end

  # if mhidx == sfidx
  # ah = sort(union([sfidx;], certainidx))

  return certainidx, allelements, activehypo, mhidx
end
function assembleHypothesesElements!(mh::Nothing,
                                     maxlen::Int,
                                     sfidx::Int,
                                     lenXi::Int,
                                     isinit::Vector{Bool}=ones(Bool, lenXi)  )
  #
  allelements = []
  activehypo = []
  mhidx = Int[]

  allidx = 1:maxlen
  certainidx = 1:lenXi
  doneall = false
  for i in certainidx
    if !doneall
      push!(allelements, allidx)
      push!(activehypo, (i,certainidx))
      doneall = true
    else
      push!(allelements, Int[])
      push!(activehypo, (i,Int[]))
    end
  end
  return certainidx, allelements, activehypo, mhidx # certainidx = allhp
end

#


text:
module Dynesty

# Write your package code here.
using PyCall
using StatsBase

export NestedSampler,
       DynamicNestedSampler,
       runplot,
       traceplot,
       cornerpoints,
       cornerplot,
       boundplot,
       cornerbound,
       sample, resample_equal

const dynesty = PyNULL()

"""
    `dyplot`
Object that holds the plotting submodule of the `dynesty` package
"""
const dyplot = PyNULL()

"""
    `NestedSampler`
Julia interface to the `dynesty` NestedSampler class. Note that we do
not pass the loglikelihood or prior transform function here. Instead this
is passed to the `sample`` call.

# Example
```julia
# define a distribution
d = MvNormal(ones(10))

loglikelihood(x) = logpdf(d, x)

prior_transform(p) = -10.0 .+ 20.0.*p

smplr = NestedSampler(10)

res = sample(loglikelihood, prior_transform, smplr; dlogz=0.5)
```
"""
Base.@kwdef struct NestedSampler
    ndim::Int
    nlive::Int = 500
    bound::String = "multi"
    sample::String = "auto"
    periodic = nothing
    reflective = nothing
    update_interval::Float64 = 1.5
    first_update = nothing
    gradient = nothing
    walks::Int = 25
    facc::Float64 = 0.5
    slices::Int = 5
    fmove::Float64 = 0.9
    max_move::Int = 100
end

function NestedSampler(ndim::Int; kwargs...)
    return NestedSampler(;ndim=ndim, kwargs...)
end

"""
    `DyamicNestedSampler`
Julia interface to the `dynesty` DynamicNestedSampler class. Note that we do
not pass the loglikelihood or prior transform function here. Instead this
is passed to the `sample`` call.

# Example
```julia
# define a distribution
d = MvNormal(ones(10))

loglikelihood(x) = logpdf(d, x)

prior_transform(p) = -10.0 .+ 20.0.*p

smplr = NestedSampler(10)

# sample using dynamic nested sampling with 500 initial live points
res = sample(loglikelihood, prior_transform, smplr; dlogz=0.5, nlive_init=500)
```
"""
Base.@kwdef struct DynamicNestedSampler
    ndim::Int
    bound::String = "multi"
    sample::String = "auto"
    periodic = nothing
    reflective = nothing
    update_interval::Float64 = 1.5
    first_update = nothing
    gradient = nothing
    walks::Int = 25
    facc::Float64 = 0.5
    slices::Int = 5
    fmove::Float64 = 0.9
    max_move::Int = 100
end

function DynamicNestedSampler(ndim::Int; kwargs...)
    return DynamicNestedSampler(;ndim=ndim, kwargs...)
end

"""
    DynestyOutput
A objects that holds the dynesty output as well as the sampler.
This object can be passed to the plotting functions to produce the usual
dynesty diagnostic plots.

# Notes
We have also implemented a limited number of Julia's `Dictionary` interface so you
can access the output in the usual Dynesty manner

# Example
```julia
# define a distribution
d = MvNormal(ones(10))

loglikelihood(x) = logpdf(d, x)

# We only look at finite region of parameter space
prior_transform(p) = -5.0 .+ 10.0.*p

smplr = NestedSampler(10)

# sample using dynamic nested sampling with 500 initial live points
res = sample(loglikelihood, prior_transform, smplr; dlogz=0.5, nlive_init=500)

# fetch the samples
res[:samples]

# print the keys
keys(res)
```

"""
struct DynestyOutput{D,S}
    dict::D
    sampler::S
end

Base.getindex(d::DynestyOutput, key) = getindex(d.dict, key)
Base.keys(d::DynestyOutput) = keys(d.dict)
Base.get(f::Function, d::DynestyOutput, key) = get(f, d.dict, key)
Base.values(d::DynestyOutput) = values(d.dict)

"""
    `sample(loglikelihood, prior_transform, s::NestedSampler; kwargs...)`

Runs dynesty's NestedSampler algorithm with the specified loglikelihood and prior_transform.
The loglikelihood and prior_transform are functions. For the specific relevant kwargs see the
dynesty documentation at [https://dynesty.readthedocs.io/]
"""
function StatsBase.sample(loglikelihood, prior_transform, s::NestedSampler; kwargs...)
    dysampler = dynesty.NestedSampler(loglikelihood,
                                      prior_transform,
                                      s.ndim,;
                                      nlive = s.nlive,
                                      bound = s.bound,
                                      sample = s.sample,
                                      periodic = s.periodic,
                                      reflective = s.reflective,
                                      update_interval = s.update_interval,
                                      first_update = s.first_update,
                                      gradient = s.gradient,
                                      walks = s.walks,
                                      facc = s.facc,
                                      slices = s.slices,
                                      fmove = s.fmove,
                                      max_move = s.max_move
                                    )
    dysampler.run_nested(;kwargs...)
    return DynestyOutput(dysampler.results, dysampler)
end

"""
    `sample(loglikelihood, prior_transform, s::DynamicNestedSampler; kwargs...)`

Runs dynesty's DynamicNestedSampler algorithm with the specified loglikelihood and prior_transform.
The loglikelihood and prior_transform are functions. For the specific relevant kwargs see the
dynesty documentation at [https://dynesty.readthedocs.io/]
"""
function StatsBase.sample(loglikelihood, prior_transform, s::DynamicNestedSampler; kwargs...)
    dysampler = dynesty.DynamicNestedSampler(loglikelihood,
                                      prior_transform,
                                      s.ndim,;
                                      bound = s.bound,
                                      sample = s.sample,
                                      periodic = s.periodic,
                                      reflective = s.reflective,
                                      update_interval = s.update_interval,
                                      first_update = s.first_update,
                                      gradient = s.gradient,
                                      walks = s.walks,
                                      facc = s.facc,
                                      slices = s.slices,
                                      fmove = s.fmove,
                                      max_move = s.max_move
                                    )
    dysampler.run_nested(;kwargs...)
    return DynestyOutput(dysampler.results, dysampler)
end

"""
    resample_equal(res::DynestyOutput, nsamples::Int)
Resample the `dynesty` nested sampling run so that the samples have equal weighting.
This uses the `StatsBase` algorithm under the hood.

The results are a vector of vectors where the inner vector corresponds to the samples.
"""
function resample_equal(res::DynestyOutput, nsamples::Int)
    samples, weights = res["samples"].T, exp.(res["logwt"].T .- res["logz"][end])
    println(size(eachcol(samples)))
    sample(collect(eachcol(samples)), Weights(weights), nsamples)
end

"""
    merge(args::DynestyOutput...; print_progres=true)
Runs dynesty's merge_runs to combine multiple separate dynesty runs.
"""
function Base.merge(args::DynestyOutput...; print_progress=true)
    runs = getproperty.([args...], :dict)
    return dynesty.utils.merge_runs(runs; print_progress)
end

# internal ess estimator
function ess(res::DynestyOutput)
    weights = exp.(res["logwt"] .- res["logz"][end])
    return inv(sum(abs2, weights))
end

# Make the dynesty plots with a bit of fun metaprogramming
PLOTS = [:runplot, :traceplot, :cornerpoints, :cornerplot, :boundplot, :cornerbound]
for P in PLOTS
    @eval begin
        """
            $($P)(d::DynestyOutput; kwargs...)
        Produces the $($P) plot from the dynesty.plotting module.
        For a list of possible kwargs see the dynesty documentation at
        [dynesty.readthedocs.io]
        """
        function ($P)(d::DynestyOutput, args...; kwargs...)
            f = getproperty(dyplot, Symbol($P))
            f(d.dict, args...; kwargs...)
        end
    end
end

function __init__()
    copy!(dynesty, pyimport_conda("dynesty","dynesty"))
    copy!(dyplot, pyimport("dynesty.plotting"))
    # Define a hack to get merge to work without doing silly dict conversion
    # Not needed for new versions
    #py"""
    #import dynesty
    #def merge(x, print_progress=True):
    #    xres = [dynesty.results.Results(r) for r in x]
    #    mruns = dynesty.utils.merge_runs(xres, print_progress)
    #    return mruns
    #"""
end

end


text:
module AlbedoProperties

using Reexport
@reexport using Dates, Statistics, DrWatson, ClimateBase

using NCDatasets: NCDataset, dimnames, NCDatasets
export NCDataset, dimnames

Time = ClimateBase.Time

using DimensionalData
export DimensionalData
AbDimArray = DimensionalData.AbstractDimArray
export AbDimArray
include("core.jl")
include("statistics.jl")
include("residuals.jl")
include("albedo.jl")

end

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
using Plots

function hailLength(x::Int)
    n = 0
    while x != 1
        if x % 2 == 0
            x = Int(x/2)
        else
            x = 3x +1
        end
        n += 1
    end
    return n
end

lengths = [hailLength(x0) for x0 in 2:10^7]

histogram(lengths, bins=1000, normed=:true,
    fill=(:blue, true), la=0, legend=:none,
    xlims=(0, 500), ylims=(0, 0.012),
    xlabel="Length", ylabel="Frequency")


text:
using Test
using ContentTrees
using BSDiff
using Tar

registry = expanduser("~/.julia/registries/General")
new_tarball = open(Tar.rewrite, `git -C $registry archive 'master@{now}'`)
old_tarball = open(Tar.rewrite,	`git -C $registry archive 'master@{1 week ago}'`)
patch_file = bsdiff(old_tarball, new_tarball)

old_root = tempname()
old_hash = extract_tree(old_tarball, old_root)
new_root = tempname()
patch_tree(old_root, new_root, patch_file)


text:
using CubicEoS
# using CubicEoSDatabase: Data
using LinearAlgebra

function vt_flash_helmoltz_diff_saturation_tries(; mix, nmol, volume, RT)
    issinglephase, vt_stab_tries = vt_stability(mix, nmol, volume, RT)
    @assert !issinglephase

    function get_state(conc, sat)
        return [(conc * (volume * sat) ./ nmol)..., sat]
    end

    grad = fill(NaN, ncomponents(mix) + 1)
    for sat₁ in 0.01:0.01:0.99
        ΔA_vec = fill(NaN, length(vt_stab_tries))
        error_vec = fill(false, length(vt_stab_tries))
        for (j, vt_try) in enumerate(vt_stab_tries)
            state = get_state(vt_try.concentration, sat₁)
            ΔA, error = NaN, false
            try
                ΔA, grad = helmholtz_diff!(state, grad)
            catch e
                @warn "helmholtz_diff! error" sat₁ e
                ΔA, error = NaN, true
            end
            ΔA_vec[j] = ΔA
            error_vec[j] = error
        end
        print(sat₁)
        print('\t', join(error_vec, '\t'))
        print('\t', join(ΔA_vec, '\t'))
        println()
    end
    return nothing
end

function vt_flash_helmoltz_diff_saturation(
    ;
    mix,
    nmol,
    volume,
    RT,
    conc₁,
    helmholtz_diff!::Function,
    satinit,
)
    state = fill(NaN, ncomponents(mix) + 1)
    grad = similar(state)

    function upd_state(conc₁, sat₁)
        N₁ = conc₁ * (volume * sat₁)
        state[1:end-1] .= N₁ ./ nmol
        state[end] = sat₁
        return nothing
    end

    sat₁ = satinit
    scale = 1
    scale_step = 0.9
    upd_state(conc₁, sat₁)
    for step in 1:50
        dir = 1 * [(conc₁ * (volume * sat₁) ./ nmol)..., 1]
        dir ./= norm(dir, 2)
        scalemax = constrain_step(state, dir) / 0.9
        satmax = state[end] + scalemax * dir[end]

        sat₁ = satinit * scale
        upd_state(conc₁, sat₁)
        ΔA, error = NaN, false
        try
            ΔA, grad = helmholtz_diff!(state, grad)
        catch e
            @warn "helmholtz_diff! error" sat₁ e
            ΔA, error = NaN, true
        end
        scale *= scale_step

        print(join([sat₁, satmax], '\t'))
        print('\t', error)
        print('\t', join([ΔA, norm(grad)], '\t'))
        println()
    end
    # upd_state(conc₁, 0.8)
    # ΔA, grad = helmholtz_diff!(state, grad)
    return nothing
end

mix = load(
    BrusilovskyEoSMixture;
    names=("methane", "n-pentane"),
    # component_dbs=(Data.martinez(), Data.brusilovsky_comp()),
    # mix_eos_db=Data.brusilovsky_mix()
)
# dump(mix)
# # exit()
# [0.547413, 0.452587] mole fractions from (Mikyska, 2013, Example 2)
# V = 1e-6 from old vtflash
# (5e3 * V) is (5 kmol m⁻³ * V) just an ordinary value
# T = 370 Kelvins at 5 kmol m⁻³ should correspond to two-phase
# T = 440 Kelvins at 5 kmol m⁻³ should correspond to single-phase

V = 1e-6
N = (5000 * V) .* [0.547413, 0.452587]
RT = 370 * CubicEoS.GAS_CONSTANT_SI

# closures
constrain_step, helmholtz_diff_grad!, helmholtz_diff! = CubicEoS.vt_flash_closures(mix, N, V, RT)

# vtstability
issinglephase, vt_stab_tries = vt_stability(mix, N, V, RT)
conc₁ = CubicEoS.__vt_flash_init_conc_choose(vt_stab_tries)
println(stderr, "Is single phase: ", issinglephase)
println(stderr, "Concentration from vt_stability ", conc₁)
# for vttry in vt_stab_tries
#     dump(vttry)
# end

# initial state debug
# state = fill(NaN, ncomponents(mix) + 1)
# found_initial = CubicEoS.vt_flash_initial_state!(
#     state, N, V, conc₁, helmholtz_diff!, constrain_step;
#     sat₁max=0.1,
#     steps=100,
#     step_scale=0.8,
#     helmholtz_thresh=-1e-5,
# )
# println(stderr, "vt_flash_initial_state", " state is found: ", found_initial)
# println(helmholtz_diff!(state, similar(state)))

# ΔA from saturation and initial concentration
# vt_flash_helmoltz_diff_saturation(
#     ;
#     mix=mix,
#     nmol=N,
#     volume=V,
#     RT=RT,
#     conc₁=conc₁,
#     helmholtz_diff! =helmholtz_diff!,
#     satinit=1,
# )

# vt_flash_helmoltz_diff_saturation_tries(mix=mix, nmol=N, RT=RT, volume=V)

# ∇P = fill(NaN, ncomponents(mix) + 1)
# CubicEoS.__vt_flash_pressure_gradient!(∇P, mix, N, V, RT)
# dump(∇P)

# ℍ = fill(NaN, (ncomponents(mix) + 1, ncomponents(mix) + 1))
# # state = fill(NaN, size(ℍ)[1])
# CubicEoS.__vt_flash_hessian!(ℍ, state, mix, N, V, RT)
# dump(ℍ)

state = CubicEoS.vt_flash_newton(mix, N, V, RT)
dump(state)

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using TextWrap

text = """Reformat the single paragraph in 'text' to fit in lines of no more
    than 'width' columns, and return a new string containing the entire
    wrapped paragraph.  As with wrap(), tabs are expanded and other
    whitespace characters converted to space.  See TextWrapper class for
    available keyword args to customize wrapping behaviour."""

println("# Wrapped at 80 chars")
print_wrapped(text, width=80)
println("\n\n# Wrapped at 70 chars")
print_wrapped(text, width=70)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Quaternions: argq
using LinearAlgebra

let # test accessors
    c = rand(4)
    @test c[1] == real(quat(c...))
    @test c[2:4] == imag(quat(c...))
end

# creating random examples
sample(QT::Type{Quaternion{T}}) where {T <: Integer} = QT(rand(-100:100, 4)..., false)
sample(QT::Type{Quaternion{T}}) where {T <: AbstractFloat} = QT(rand(Bool) ? quatrand() : nquatrand())
sample(CT::Type{Complex{T}}) where {T <: Integer} = CT(rand(-100:100, 2)...)
sample(CT::Type{Complex{T}}) where {T <: AbstractFloat} = CT(randn(2)...)
sample(T, n) = T[sample(T) for _ in 1:n]

# test algebraic properties of quaternions
for _ in 1:10, T in (Float32, Float64, Int32, Int64)
    q, q1, q2, q3 = sample(Quaternion{T}, 4)

    # skewfield
    test_group(q1, q2, q3, +, zero(q), -)
    test_group(q1, q2, q3, *, one(q), inv)
    test_multiplicative(q1, q2, *, norm)

    # complex embedding
    c1, c2 = sample(Complex{T}, 2)
    test_multiplicative(c1, c2, *, Quaternion)
    test_multiplicative(c1, c2, +, Quaternion)
end

let # test rotations
    qx = qrotation([1, 0, 0], pi / 4)
    @test qx * qx ≈ qrotation([1, 0, 0], pi / 2)
    @test qx^2 ≈ qrotation([1, 0, 0], pi / 2)
    theta = pi / 8
    qx = qrotation([1, 0, 0], theta)
    c = cos(theta); s = sin(theta)
    Rx = [1 0 0; 0 c -s; 0 s c]
    @test rotationmatrix(qx) ≈ Rx
    theta = pi / 6
    qy = qrotation([0, 1, 0], theta)
    c = cos(theta); s = sin(theta)
    Ry = [c 0 s; 0 1 0; -s 0 c]
    @test rotationmatrix(qy) ≈ Ry
    theta = 4pi / 3
    qz = qrotation([0, 0, 1], theta)
    c = cos(theta); s = sin(theta)
    Rz = [c -s 0; s c 0; 0 0 1]
    @test rotationmatrix(qz) ≈ Rz

    @test rotationmatrix(qx * qy * qz) ≈ Rx * Ry * Rz
    @test rotationmatrix(qy * qx * qz) ≈ Ry * Rx * Rz
    @test rotationmatrix(qz * qx * qy) ≈ Rz * Rx * Ry

    a, b = qrotation([0, 0, 1], deg2rad(0)), qrotation([0, 0, 1], deg2rad(180))
    @test slerp(a, b, 0.0) ≈ a
    @test slerp(a, b, 1.0) ≈ b
    @test slerp(a, b, 0.5) ≈ qrotation([0, 0, 1], deg2rad(90))

    @test angle(qrotation([1, 0, 0], 0)) ≈ 0
    @test angle(qrotation([0, 1, 0], pi / 4)) ≈ pi / 4
    @test angle(qrotation([0, 0, 1], pi / 2)) ≈ pi / 2

    let # test numerical stability of angle
        ax = randn(3)
        for θ in [1e-9, π - 1e-9]
            q = qrotation(ax, θ)
            @test angle(q) ≈ θ
        end
    end
end

for _ in 1:100
    let # test specialfunctions
        c = Complex(randn(2)...)
        q, q2 = sample(Quaternion{Float64}, 4)
        unary_funs = [exp, log, sin, cos, sqrt, inv, conj, abs2, norm]
        # since every quaternion is conjugate to a complex number,
        # one can establish correctness as follows:
        for fun in unary_funs
            @test fun(Quaternion(c)) ≈ Quaternion(fun(c))
            @test q2 * fun(q) * inv(q2) ≈ fun(q2 * q * inv(q2))
        end

        @test exp(log(q)) ≈ q
        @test exp(zero(q)) ≈ one(q)
    end

    let # test qrotation and angleaxis inverse
        ax = randn(3); ax = ax / norm(ax)
        Θ = π * rand()
        q = qrotation(ax, Θ)
        @test angle(q) ≈ Θ
        @test axis(q) ≈ ax
        @test angleaxis(q)[1] ≈ Θ
        @test angleaxis(q)[2] ≈ ax
    end

    let # test argq
        q, q2 = sample(Quaternion{Float64}, 2)
        @test q2 * argq(q) * inv(q2) ≈ argq(q2 * q * inv(q2))
        v = Quaternion(0, randn(3)...)
        @test argq(v) * norm(v) ≈ v
    end

    let # test normalize
        q = quatrand()
        @test norm(normalize(q)) ≈ 1
        @test normalize(q).norm
        @test q ≈ norm(q) * normalize(q)
        qn = nquatrand()
        @test qn.norm
        @test normalize(qn) === qn
    end

    let # test rotation <=> rotationmatrix
        q1 = nquatrand()
        q2 = qrotation(rotationmatrix(q1), q1)
        @test q1 ≈ q2
    end

    let # test slerp and linpol if q1 = 1
        q1 = quat(1, 0, 0, 0.)
        # there are numerical stability issues with slerp atm
        θ = clamp(rand() * 3.5, deg2rad(5e-1), π)
        ax = randn(3)
        q2 = qrotation(ax, θ)
        t = rand()
        slerp(q1, q2, 0.) ≈ q1
        @test slerp(q1, q2, 0.) ≈ q1
        @test slerp(q1, q2, 1.) ≈ q2
        @test slerp(q1, q2, t) ≈ qrotation(ax, t * θ)
        @test norm(slerp(q1, q2, t)) ≈ 1
        @test slerp(q1, q2, 0.5) ≈ qrotation(ax, 0.5 * θ)
        @test linpol(q1, q2, 0.5) ≈ qrotation(ax, 0.5 * θ)

    end
    let # test conjugation invariance
        q, q1, q2 = sample(Quaternion{Float64}, 3)
        ⊗(s, t) = s * t * inv(s)
        t = rand()
        @test q ⊗ slerp(q1, q2, t) ≈ slerp(q ⊗ q1, q ⊗ q2, t)
        @test q ⊗ linpol(q1, q2, t) ≈ linpol(q ⊗ q1, q ⊗ q2, t)
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_dashdraggable2

"""
    ''_dashdraggable2(;kwargs...)
    ''_dashdraggable2(children::Any;kwargs...)
    ''_dashdraggable2(children_maker::Function;kwargs...)

A DashDraggable2 component.
ExampleComponent is an example component.
It takes a property, `label`, and
displays it.
It renders an input with the property `value`
which is editable by the user.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `value` (String; optional): The value displayed in the input.
"""
function ''_dashdraggable2(; kwargs...)
        available_props = Symbol[:children, :id, :value]
        wild_props = Symbol[]
        return Component("''_dashdraggable2", "DashDraggable2", "dash_draggable", available_props, wild_props; kwargs...)
end

''_dashdraggable2(children::Any; kwargs...) = ''_dashdraggable2(;kwargs..., children = children)
''_dashdraggable2(children_maker::Function; kwargs...) = ''_dashdraggable2(children_maker(); kwargs...)

max_stars_repo_name: tcgvn/dash-draggable

max_issues_repo_name: tcgvn/dash-draggable

max_forks_repo_name: tcgvn/dash-draggable


text:
module TestLogUniform
using Test
using Distributions
import Random

@testset "LogUniform" begin
    rng = Random.MersenneTwister(0)

    @test pdf(LogUniform(1f0, 2f0), 1) isa Float32
    @test pdf(LogUniform(1, 2), 1f0) isa Float32
    @test pdf(LogUniform(1, 2), 1) isa Float64
    @test quantile(LogUniform(1, 2), 1) isa Float64
    @test quantile(LogUniform(1, 2), 1f0) isa Float32
    @testset "$f" for f in [pdf, cdf, quantile, logpdf, logcdf]
        @test @inferred(f(LogUniform(1,2), 1)) isa Float64
        @test @inferred(f(LogUniform(1,2), 1.0)) isa Float64
        @test @inferred(f(LogUniform(1.0,2), 1.0)) isa Float64
        @test @inferred(f(LogUniform(1.0f0,2), 1)) isa Float32
        @test @inferred(f(LogUniform(1.0f0,2), 1f0)) isa Float32
        @test @inferred(f(LogUniform(1,2), 1f0)) isa Float32
    end

    d = LogUniform(1,10)
    @test eltype(d) === Float64
    @test 1 <= rand(rng, d) <= 10
    @test rand(rng, d) isa eltype(d)
    @test @inferred(quantile(d, 0))   ≈ 1
    @test quantile(d, 0.5) ≈ sqrt(10) # geomean
    @test quantile(d, 1)   ≈ 10
    @test mode(d) ≈ 1
    @test !insupport(d, 0)
    @test @inferred(minimum(d)) === 1
    @test @inferred(maximum(d)) === 10
    @test partype(d) === Int
    @test truncated(d, 2, 14) === LogUniform(2,10)
    @test truncated(d, 0, 8) === LogUniform(1, 8)
    @test truncated(d; upper=8) === LogUniform(1, 8)
    @test truncated(d; lower=3) === LogUniform(3, 10)

    # numbers obtained by calling scipy.stats.loguniform
    @test @inferred(std(d)        ) ≈ 2.49399867607628
    @test @inferred(mean(d)       ) ≈ 3.908650337129266
    @test @inferred(pdf(d, 1.0001)) ≈ 0.43425105679757203
    @test @inferred(pdf(d, 5     )) ≈ 0.08685889638065035
    @test @inferred(pdf(d, 9.9999)) ≈ 0.04342988248915007
    @test @inferred(cdf(d, 1.0001)) ≈ 4.342727686266485e-05
    @test @inferred(cdf(d, 5     )) ≈ 0.6989700043360187
    @test @inferred(cdf(d, 9.9999)) ≈ 0.999995657033466
    @test @inferred(median(d)     ) ≈ 3.1622776601683795
    @test @inferred(logpdf(d, 5)  ) ≈ -2.443470357682056

    for _ in 1:10
        lo = rand(rng)
        hi = lo + 10*rand(rng)
        dist = LogUniform(lo,hi)
        q = rand(rng)
        @test cdf(dist, quantile(dist, q)) ≈ q

        u = Uniform(log(lo), log(hi))
        @test exp(quantile(u, q)) ≈ quantile(dist, q)
        @test exp(median(u)) ≈ median(dist)
        x = rand(rng, dist)
        @test cdf(u, log(x)) ≈ cdf(dist, x)

        @test @inferred(entropy(dist)) ≈  Distributions.expectation(x->-logpdf(dist,x), dist)
    end

    @test kldivergence(LogUniform(1,2), LogUniform(1,2)) ≈ 0 atol=100eps(Float64)
    @test isfinite(kldivergence(LogUniform(1,2), LogUniform(1,10)))
    @test kldivergence(LogUniform(1.1,10), LogUniform(1,2)) === Inf
    @test kldivergence(LogUniform(0.1,10), LogUniform(1,2)) === Inf
    @test kldivergence(LogUniform(0.1,1), LogUniform(1,2))  === Inf
    @test @inferred(kldivergence(LogUniform(0.1f0,1), LogUniform(1,2)))  === Inf32

    for _ in 1:10
        aq = 10*rand(rng)
        ap = aq + 10*rand(rng)
        bp = ap + 10*rand(rng)
        bq = bp + 10*rand(rng)
        p = LogUniform(ap, bp)
        q = LogUniform(aq, bq)
        @test @inferred(kldivergence(p, q)) ≈
            kldivergence(Uniform(log(ap), log(bp)), Uniform(log(aq), log(bq)))
    end
end

end#module


text:
# This file is a part of DataDeps.jl. License is MIT.
## Core path determining stuff

const standard_loadpath = joinpath.([
    Base.DEPOT_PATH; homedir(); # Common all systems

    @static if Sys.iswindows()
        vcat(get.(Ref(ENV),
           ["APPDATA", "LOCALAPPDATA",
            "ProgramData", "ALLUSERSPROFILE", # Probably the same, on all systems where both exist
            "PUBLIC", "USERPROFILE"], # Home Dirs ("USERPROFILE" is probably the same as homedir()
           [String[]])...)
    else
        ["/scratch", "/staging", # HPC common folders
         "/usr/share", "/usr/local/share"] # Unix Filestructure
    end], "datadeps")

# ensure at least something in the loadpath exists when instaleld
mkpath(first(standard_loadpath))

########################################################################################################################
## Package reletive path determining

"""
    try_determine_package_datadeps_dir(filepath)

Takes a path to a file.
If that path is in a package's folder,
Then this returns a path to the deps/data dir for that package (as a Nullable).
Which may or may not exist.
If not in a package returns null
"""
function try_determine_package_datadeps_dir(filepath)::Union{String, Nothing}
    olddir=filepath
    curdir = dirname(filepath)
    while(olddir!=curdir) # At root `dirname` returns it's input
        olddir = curdir
        curdir = dirname(curdir)
        datadeps_dir = joinpath(curdir, "deps","data")
        if isdir(datadeps_dir)
            return datadeps_dir
        end
    end
    return nothing
end

############################################

"""
    preferred_paths(calling_filepath; use_package_dir=true)

returns the datadeps load_path
plus if calling_filepath is provided and `use_package_dir=true`
and is currently inside a package directory then it also includes the path to the dataseps in that folder.
"""
function preferred_paths(rel=nothing; use_package_dir=true)
    cands = String[]
    if use_package_dir
        @assert rel != nothing
        pkg_deps_root = try_determine_package_datadeps_dir(rel)
        pkg_deps_root != nothing && push!(cands, pkg_deps_root)
    end

    append!(cands, env_list("DATADEPS_LOAD_PATH", []))
    if !env_bool("DATADEPS_NO_STANDARD_LOAD_PATH", false)
        append!(cands, standard_loadpath)
    end
    cands
end

####################################################################################################################
## Permission checking stuff
@enum AccessMode F_OK=0b0000 X_OK=0b0001 W_OK=0b0010 XW_OK=0b0011 R_OK=0b0100 XR_OK=0b0101 WX_OK=0b0110 XWR_OK=0b0111

"""
    uv_access(path, mode)

Check access to a path.
Returns 2 results, first an error code (0 for all good), and second an error message.
https://stackoverflow.com/a/47126837/179081
"""
function uv_access(path, mode::AccessMode)
    local ret
    req = Libc.malloc(Base._sizeof_uv_fs)
    try
        ret = ccall(:uv_fs_access, Cint,
                (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}),
                Base.eventloop(), req, path, mode, C_NULL)
        ccall(:uv_fs_req_cleanup, Cvoid, (Ptr{Cvoid},), req)
    finally
        Libc.free(req)
    end
    return ret, ret==0 ? "OK" : Base.struverror(ret)
end

can_read_file(path) = uv_access(path, R_OK)[1] == 0

##########################################################################################################################
## Actually determining path being used (/going to be used) by a given datadep

"""
    determine_save_path(name)

Determines the location to save a datadep with the given name to.
"""
function determine_save_path(name, rel=nothing)::String
    cands = preferred_paths(rel; use_package_dir=false) #TODO Consider removing `rel` argument, it is not used
    path_ind = findfirst(cands) do path
        0 == first(uv_access(path, W_OK))
    end
    if path_ind === nothing
            @error """
            No writable path exists to save the data. Make sure there exists as writable path in your DataDeps Load Path.
            See https://www.oxinabox.net/DataDeps.jl/stable/z10-for-end-users/#The-Load-Path-1
            The current load path contains:
            """ cands
        throw(NoValidPathError("No writable path exists to save the data."))
    end
    return joinpath(cands[path_ind], name)
end

"""
    try_determine_load_path(name)

Trys to find a local path to the datadep with the given name.
If it fails then it returns nothing.
"""
function try_determine_load_path(name::String, rel)
    paths = list_local_paths(name, rel)
    paths = paths[first.(uv_access.(paths, Ref(R_OK))) .== 0] # 0 means passes
    length(paths)==0 ? nothing :  first(paths)
end

"""
    list_local_paths( name|datadep, [calling_filepath|module|nothing])

Lists all the local paths to a given datadep.
This may be an empty list
"""
function list_local_paths(name::String, rel)
    cands = preferred_paths(rel)
    joinpath.(cands, name)
end

list_local_paths(dd::AbstractDataDep, rel) = list_local_paths(dd.name, rel)


text:
function perform_blurring(M, sigma, bound="sym")
    """
        perform_blurring - gaussian blurs an image

        M = perform_blurring(M, sigma, options);

        M is the original data
        sigma is the width of blurs (in pixels)

        Copyright (c) 2007 Gabriel Peyre
    """

    if sigma == zeros(size(sigma))
        return M
    end

    if ndims(M) > 2
        for i in 1 : size(M)[3]
             M[:, :, i] = perform_blurring(M[:, :, i], sigma, bound)
        end
    end

    n = maximum(size(M))

    eta = 4
    p = round((sigma.*eta)./2.).*2 + 1
    p = min(p, (round(n/2.)*2 - 1).*ones(length(p)))

    A = [1., 1.]
    if ndims(M) == 1
        A = 1 #1D
    end

    h = compute_gaussian_filter(p.*A, sigma./(4.0 * n), n*A)
    M = perform_convolution(M, h, bound)

    return M
end

function compute_gaussian_filter(n, s, N)
    """
        compute_gaussian_filter - compute a 1D or 2D Gaussian filter.

          f = compute_gaussian_filter(n,s,N);

          'n' is the size of the filter, odd for no phase in the filter.
              (if too small it will alterate the filter).
              use n=[n1,n2] for a 2D filter or n = [n1] for a 1D filter
          's' is the standard deviation of the filter.
          'N' is the size of the big signal/image (supposed to lie in [0,1] or [0,1]x[0,1]).
              use N=[N1,N2] for a 2D filter or N = [N1] for a 1D filter

          The equation (in 1D) is
              f[k] = exp( -(x(k)^2/(2*s^2)) );
          where x spans [-1/2,1/2].

          The filter is normalised so that it sums to 1.

          Copyright (c) 2004 Gabriel Peyre
    """
    nd = 1
    if (length(n) > 1) & (n[2] > 1)
        nd = 2
    end

    if (nd == 2) & (length(s) == 1)
        s = vcat(s, s)
    end

    if (nd == 2) & (length(N) == 1)
        N = vcat(N, N)
    end

    if nd == 1
        f = build_gaussian_filter_1d(n, s, N)
    else
        f = build_gaussian_filter_2d(n, s, N)
    end
    return f
end

function build_gaussian_filter_2d(n, s, N = [])
    """
        build_gaussian_filter_2d - compute a 2D Gaussian filter.

        f = build_gaussian_filter_2d(n,s,N);

        'n' is the size of the filter, odd for no phase in the filter.
            (if too small it will alterate the filter).
        's' is the standard deviation of the filter.
        'N' is the size of the big image (supposed to lie in [0,1]x[0,1]).

        The filter is normalised so that it sums to 1.

        Copyright (c) 2004 Gabriel Peyre
    """

    # n = np.asarray(n)
    # s = np.asarray(s)
    # N = np.asarray(N)

    if length(N) == 0
        N = n
    end

    if (length(N) == 1) | (N[1] == 1)
        N = vcat(N, N)
    end

    if (length(s) == 1) | (s[1] == 1)
        s = vcat(s,s)
    end

    if length(collect(x for x in s if x <= 0)) > 0
        f = zeros(n)
        f[Int(round((n - 1)/2))] = 1
        return f
    end

    x = (collect(0 : n[1]) .- (n[1] - 1)/2.)/(N[1] - 1)
    y = (collect(0 : n[2]) .- (n[2] - 1)/2.)/(N[2] - 1)
    (Y, X) = meshgrid(y, x) #meshgrid is in the file ndgrid.jl
    f = exp(-(X.^2 ./ (2*s[1]^2)) - (Y.^2 ./ (2*s[2]^2)))
    f = f./sum(f)
    return f
end

function build_gaussian_filter_1d(n, s, N=[])
    """
        build_gaussian_filter_1d - compute a Gaussian filter.

        f = build_gaussian_filter_1d(n,s,N);

        Copyright (c) 2004 Gabriel Peyre
    """
    if length(N) == 0
        N = n
    end

    n = n[1]
    s = s[1]
    N = N[1]

    if s <= 0
        f = zeros(n)
        f[round((n - 1)/2)] = 1
        return f
    end

    x = (collect(0 : n) .- (n - 1)/2.)/(N - 1)
    f = exp(-x.^2 ./ (2*s^2))
    f = f./sum(f)
    return f
end

function perform_convolution(x, h, bound = "sym")
    """
        perform_convolution - compute convolution with centered filter.

        y = perform_convolution(x,h,bound);

        The filter 'h' is centred at 0 for odd
        length of the filter, and at 1/2 otherwise.

        This works either for 1D or 2D convolution.
        For 2D the matrix have to be square.

        'bound' is either 'per' (periodic extension)
        or 'sym' (symmetric extension).

        Copyright (c) 2004 Gabriel Peyre
    """

    if !(bound in ["sym", "per"])
        error("bound should be sym or per")
    end

    if ndims(x) == 3
      if size(x)[3] < 4
          #for color images
        y = x
        for i in 1 : size(x)[3]
          y[:, :, i] = perform_convolution(x[:, :, i], h, bound)
        end
      end
    return y
    end

    if ndims(x) == 3
      if size(x)[3] >= 4
        error("Not yet implemented for 3D array, use smooth3 instead.")
      end
    end

    n = size(x)
    p = size(h)

    nd = ndims(x)

    if nd == 1
        n = length(x)
        p = length(h)
    end

    if bound == "sym"

                #################################
        # symmetric boundary conditions #
        d1 = Array{Int64}(p)./2  # padding before
        d2 = p - d1 - 1    			    # padding after

        if nd == 1
        ################################# 1D #################################
            nx = length(x)
            xx = hcat(x[d1: -1: end], x, x[nx - 1 : -1 : nx - d2 - 1]) #A vérifier
            y = conv2(xx,h)
            y = y[p + 1 : nx - p - 1] #Il faut revérifier les indices

        elseif nd == 2
        ################################# 2D #################################
            #double symmetry
            nx = size(x)
            ny = size(x)
            xx = x
            xx = hcat(xx[d1[1] + 1 : -1 : end, :], xx, xx[nx : -1 : nx - d2[2] - 1, :])
            xx = vcat(xx[:, d1[2] : -1 : end], xx, xx[:, ny : -1 : ny - d2[2] - 1])
            y = conv2(xx, h) #mode="same"
            y = y[(2*d1[1] + 1) : (2*d1[1] + n[1] + 1), (2*d1[2] + 1) : (2*d1[2] + n[2] + 1)]
        end

    else

        ################################
        # periodic boundary conditions #

        if p > n
            error("h filter should be shorter than x.")
        end
        n = [n[1], n[2]]
        p = [p[1], p[2]]
        d = Array{Int64, 1}(floor((p - 1)/2.))
        if nd == 1
            h = hcat(h[d : end], hcat(zeros(n - p), h[1 : d - 1]))
            y = real(ifft((fft(x)*x).*(fft(h)*h)*((fft(x)*x).*(fft(h)*h))))
        else
            h = vcat(h[d[1] + 1 : end, :], vcat(zeros(n[1] - p[1], p[2]), h[1 : d[1], :]))
            h = hcat(h[:, d[2] + 1 : end], hcat(zeros(n[1], n[2] - p[2]), h[:, 1 : d[2] ]))
            y = real(plan_ifft((plan_fft(x)*x).*(plan_fft(h)*h))*((plan_fft(x)*x).*(plan_fft(h)*h)))
        end
    end
    return y
end


text:
# Interface

import Base: replace
import Manifolds: identity_element

export makeCoordsFromPoint, makePointFromCoords, getNumberCoords
export identity_element
export setPointPartial!, setPointsMani!
export replace

# Deprecate in favor of TranslationGroup instead, also type piracy
Manifolds.identity_element(::Euclidean{Tuple{N}}, val::AbstractVector{T}) where {N, T <: Number} = zeros(T, N)
Manifolds.identity_element(::Circle, val::AbstractVector{T}) where {T <: Real} = zeros(T, 1)

"""
    $SIGNATURES

Helper function to convert coordinates to a desired on-manifold point.

DevNotes
- FIXME need much better consolidation or even removal of this function entirely.
  - This function makes too strong an assumption of groups

Notes
- `u0` is used to identify the data type for a point
- Pass in a different `exp` if needed.
"""
makePointFromCoords(M::MB.AbstractManifold,  # Manifolds.AbstractGroupManifold
                    coords::AbstractVector{<:Real},
                    u0=zeros(manifold_dimension(M)),
                    ϵ=identity_element(M,u0),
                    retraction_method::AbstractRetractionMethod=ExponentialRetraction()  ) = retract(M, ϵ, hat(M, ϵ, coords), retraction_method)
#

# should perhaps just be dispatched for <:AbstractGroupManifold
# only works for AbstractGroupManifold (have an identity)
function makeCoordsFromPoint( M::MB.AbstractManifold,
                              pt::P,
                              ϵ = identity_element(M, pt) ) where P
  #

  vee(M, ϵ, log(M, ϵ, pt))
end

# Sphere(2) has 3 coords, even though the manifolds only has 2 dimensions (degrees of freedom)
getNumberCoords(M::MB.AbstractManifold, p) = length(makeCoordsFromPoint(M,p))

# TODO DEPRECATE
# related _pointsToMatrixCoords
function _matrixCoordsToPoints( M::MB.AbstractManifold,
                                pts::AbstractMatrix{<:Real},
                                u0  )
  #
  # ptsArr = Vector{Vector{Float64}}(undef, size(pts, 2))
  # @cast ptsArr[j][i] = pts[i,j]
  vecP = Vector{typeof(u0)}(undef, size(pts, 2))
  for j in 1:size(pts,2)
    pt = pts[:,j]
    vecP[j] = makePointFromCoords(M, pt, u0)
  end
  return vecP
end

function _pointsToMatrixCoords(M::MB.AbstractManifold, pts::AbstractVector{P}) where P
  mat = zeros(manifold_dimension(M), length(pts))
  ϵ = identity_element(M, pts[1])
  for (j,pt) in enumerate(pts)
    mat[:,j] = vee(M, ϵ, log(M, ϵ, pt))
  end

  return mat
end

# asPartial=true indicates that src coords are smaller than dest coords, and false implying src has dummy values in placeholder dimensions
function setPointPartial!(Mdest::AbstractManifold,
                          dest,
                          Msrc::AbstractManifold,
                          src,
                          partial::AbstractVector{<:Integer},
                          asPartial::Bool=true )
  #

  dest_ = AMP.makeCoordsFromPoint(Mdest,dest)
  # e0 = identity_element(Mdest, dest)
  # dest_ = vee(Mdest, e0, log(Mdest, e0, dest))

  src_ = AMP.makeCoordsFromPoint(Msrc,src)
  # e0s = identity_element(Msrc, src)
  # src_ = vee(Msrc, e0s, log(Msrc, e0s, src))

  # do the copy in coords
  dest_[partial] .= asPartial ? src_ : view(src_, partial)

  # update points base in original
  dest__ = makePointFromCoords(Mdest, dest_, dest)
  # dest__ = exp(Mdest, e0, hat(Mdest, e0, dest_))
  setPointsMani!(dest, dest__)

  #
  return dest
end

setPointsMani!(dest::AbstractVector, src::AbstractVector) = (dest .= src)
setPointsMani!(dest::AbstractMatrix, src::AbstractMatrix) = (dest .= src)
function setPointsMani!(dest::AbstractVector, src::AbstractMatrix)
  @assert size(src,2) == 1 "Workaround setPointsMani! currently only allows size(::Matrix, 2) == 1"
  setPointsMani!(dest, src[:])
end
function setPointsMani!(dest::AbstractMatrix, src::AbstractVector)
  @assert size(dest,2) == 1 "Workaround setPointsMani! currently only allows size(::Matrix, 2) == 1"
  setPointsMani!(view(dest,:,1), src)
end

function setPointsMani!(dest::AbstractVector, src::AbstractVector{<:AbstractVector})
  @assert length(src) == 1 "Workaround setPointsMani! currently only allows Vector{Vector{P}}(...) |> length == 1"
  setPointsMani!(dest, src[1])
end

function setPointsMani!(dest::ProductRepr, src::ProductRepr)
  for (k,prt) in enumerate(dest.parts)
    setPointsMani!(prt, src.parts[k])
  end
end

# default replace non-partial/non-marginal values
# Trivial case where no information from destination is kept, only from src.
function Base.replace(::ManifoldKernelDensity{M,<:BallTreeDensity,Nothing},
                      src::ManifoldKernelDensity{M,<:BallTreeDensity,Nothing}
                      ) where {M<:AbstractManifold}
  #
  src
end

# replace dest non-partial with incoming partial values
function Base.replace( dest::ManifoldKernelDensity{M,<:BallTreeDensity,Nothing},
                        src::ManifoldKernelDensity{M,<:BallTreeDensity,<:AbstractVector}
                      ) where {M<:AbstractManifold}
  #
  pl = src._partial
  oldPts = getPoints(dest.belief)
  # get source partial points only
  newPts = getPoints(src.belief)
  @assert size(newPts,2) == size(oldPts,2) "this replace currently requires the number of points to be the same, dest=$(size(oldPts,2)), src=$(size(newPts,2))"
  for i in 1:size(oldPts, 2)
    oldPts[pl,i] .= newPts[pl,i]
  end
  # and new bandwidth
  oldBw = getBW(dest.belief)[:,1]
  oldBw[pl] .= getBW(src.belief)[pl,1]

  # finaly update the belief with a new container
  newBel = kde!(oldPts, oldBw)

  # also set the metadata values
  ipc = deepcopy(dest.infoPerCoord)
  ipc[pl] .= src.infoPerCoord[pl]

  # and _u0 point is a bit more tricky
  c0 = vee(dest.manifold, dest._u0, log(dest.manifold, dest._u0, dest._u0))
  c_ = vee(dest.manifold, dest._u0, log(dest.manifold, dest._u0, src._u0))
  c0[pl] .= c_[pl]
  u0 = exp(dest.manifold, dest._u0, hat(dest.manifold, dest._u0, c0))

  # return the update destimation ManifoldKernelDensity object
  ManifoldKernelDensity(dest.manifold, newBel, nothing, u0, infoPerCoord=ipc)
end

# replace partial/marginal with different incoming partial values
function Base.replace(dest::ManifoldKernelDensity{M,<:BallTreeDensity,<:AbstractVector},
                      src::ManifoldKernelDensity{M,<:BallTreeDensity,<:AbstractVector}
                      ) where {M<:AbstractManifold}
  #
  pl = src._partial
  oldPts = getPoints(dest.belief)
  # get source partial points only
  newPts = getPoints(src.belief)
  @assert size(newPts,2) == size(oldPts,2) "this replace currently requires the number of points to be the same, dest=$(size(oldPts,2)), src=$(size(newPts,2))"
  for i in 1:size(oldPts, 2)
    oldPts[pl,i] .= newPts[pl,i]
  end
  # and new bandwidth
  oldBw = getBW(dest.belief)[:,1]
  oldBw[pl] .= getBW(src.belief)[pl,1]

  # finaly update the belief with a new container
  newBel = kde!(oldPts, oldBw)

  # also set the metadata values
  ipc = deepcopy(dest.infoPerCoord)
  ipc[pl] .= src.infoPerCoord[pl]

  # and _u0 point is a bit more tricky
  c0 = vee(dest.manifold, dest._u0, log(dest.manifold, dest._u0, dest._u0))
  c_ = vee(dest.manifold, dest._u0, log(dest.manifold, dest._u0, src._u0))
  c0[pl] .= c_[pl]
  u0 = exp(dest.manifold, dest._u0, hat(dest.manifold, dest._u0, c0))

  # and update the partial information
  pl_ = union(dest._partial, pl)

  # return the update destimation ManifoldKernelDensity object
  if length(pl_) == manifold_dimension(dest.manifold)
    # no longer a partial/marginal
    return ManifoldKernelDensity(dest.manifold, newBel, nothing, u0, infoPerCoord=ipc)
  else
    # still a partial
    return ManifoldKernelDensity(dest.manifold, newBel, pl_, u0, infoPerCoord=ipc)
  end
end

##============================================================================================================
## New Manifolds.jl aware API -- TODO find the right file placement
##============================================================================================================

# TODO, hack, use the proper Manifolds.jl intended vectoration methods instead
_makeVectorManifold(::MB.AbstractManifold, arr::AbstractArray{<:Real}) = arr
_makeVectorManifold(::MB.AbstractManifold, val::Real) = [val;]
_makeVectorManifold(::M, prr::ProductRepr) where {M <: typeof(SpecialEuclidean(2))} = coords(M, prr)
_makeVectorManifold(::M, prr::ProductRepr) where {M <: typeof(SpecialEuclidean(3))} = coords(M, prr)

#


text:
using Test
using CGP

const global D = 2
const global smax = 1000
const global snum = 4

files = ["cfg/default.yaml", "cfg/base.yaml", "cfg/test.yaml", "cfg/classic.yaml", "cfg/atari.yaml"]

function test_inps(f, inps)
    out = f(inps...)
    @test all(out == f(inps...))
    minout = minimum(out); maxout = maximum(out)
    # @printf("%0.2f, %0.2f, %s, ", minout, maxout, typeof(out))
    @test (minout >= -1.0) && (maxout <= 1.0)
end

@testset "Functions" begin
    for file in files
        CGP.Config.init(file)
        for f in CGP.Config.functions
            @testset "$file $f" begin
                println("$file $f")
                for d in 1:D
                    for s in Int64.(round.(range(1, stop=smax, length=snum)))
                        for cval in -1:1:1
                            cval = Float64(cval)
                            test_inps(f, [cval, cval, cval])
                            test_inps(f, [cval,
                                          cval * ones(repeat([s], inner=d)...),
                                          cval])
                            test_inps(f, [cval * ones(repeat([s], inner=d)...),
                                          cval, cval])
                            for dy in 1:D
                                test_inps(f, [cval * ones(repeat([s], inner=d)...),
                                              cval * ones(repeat([s], inner=dy)...),
                                              cval])
                            end
                        end
                        test_inps(f, [2*rand()-1,
                                      2*rand()-1,
                                      2*rand()-1])
                        test_inps(f, [2*rand()-1,
                                      2*rand(repeat([s], inner=d)...) .- 1,
                                      2*rand()-1])
                        test_inps(f, [2*rand(repeat([s], inner=d)...) .- 1,
                                      2*rand()-1,
                                      2*rand()-1])
                        for dy in 1:D
                            test_inps(f, [2*rand(repeat([s], inner=d)...) .- 1,
                                          2*rand(repeat([s], inner=dy)...) .- 1,
                                          2*rand()-1])
                        end
                    end
                end
            end
        end
        CGP.Config.reset()
    end
end


text:
@testset "mutable graph" begin

g = complete_digraph(5)
abm = ABM(Agent5, GraphSpace(g))
for n in 1:5
    for _ in 1:5
        add_agent!(n, abm, rand(abm.rng))
    end
end
@test nv(abm) == 5
@test abm.space.s[1] == 1:5
@test abm.space.s[2] == 6:10

rem_node!(abm, 2)
@test nv(abm) == 4
@test nagents(abm) == 20
# Last node became 2nd node (swapped places as per LightGraphs.jl)
@test abm.space.s[2] == 21:25
@test length(abm.space.s) == 4

n = add_node!(abm)
@test n == 5
@test length(abm.space.s) == 5
@test nv(abm) == 5
add_edge!(abm, n, 2)

a = add_agent!(n, abm, rand(abm.rng))
ids = nearby_ids(a, abm)
@test sort(ids) == 21:25

model = ABM(Agent5, GraphSpace(SimpleGraph()))
@test_throws ArgumentError add_agent!(model, rand(model.rng))
add_node!(model)
@test nv(model) == 1
@test add_edge!(model, 1, 2) == false
add_node!(model)
@test add_edge!(model, 1, 2) == true

end


text:
using Test, StatsModels
using DataFrames
using StatsBase

using Unfold
include("test_utilities.jl")

data,evts = loadtestdata("test_case_3a") #
data_r = reshape(data,(1,:))

data_e,times = Unfold.epoch(data=data_r,tbl=evts,τ=(0.,1),sfreq=20) # cut the data into epochs
basisfunction = firbasis(τ=(0.,1),sfreq=20,name="basisA")

# f  = @formula 0~ 1 * beta_0 + contA*beta_1
# fit => beta_0 = -1.5, beta_1 = 2
# predict => (contA=0.5) y_hat = 1 * -1.5 + 0.5 * 2 = -0.5

f  = @formula 0~1 # 1
m_mul,m_mul_results = fit(UnfoldLinearModel,f,evts,data_e,times)
m_tul,m_tul_results = fit(UnfoldLinearModel,f,evts,data_r,basisfunction)

conditionA  = [0,1.]
continuousA = [-1.,0,1.]

tmp = reshape([ [x,y]  for x=conditionA, y=continuousA ],length(conditionA)*length(continuousA))
evts_grid = DataFrame(hcat(tmp...)')
rename!(evts_grid,["conditionA","continuousA"])

yhat_tul = predict(m_tul,evts_grid)
yhat_mul = predict(m_mul,evts_grid,times)

@test unique(yhat_mul.times) == times

@test size(m_tul_results)[1] * size(evts_grid)[1] == size(yhat_tul)[1]
@test size(m_mul_results)[1] * size(evts_grid)[1] == size(yhat_mul)[1]

@test yhat_mul[end-3,:yhat] ≈ mean(data[data.!=0])
@test yhat_tul[end-3,:yhat] ≈ mean(data[data.!=0])

@test yhat_mul.times[1] == 0.

f  = @formula 0~1+conditionA+continuousA# 1
m_tul,m_tul_results = fit(UnfoldLinearModel,f,evts,data_r,basisfunction)
m_mul,m_mul_results = fit(UnfoldLinearModel,f,evts,data_e,times)
yhat_tul = predict(m_tul,evts_grid)
yhat_mul = predict(m_mul,evts_grid,times)
#@test predict(m_mul,evts).yhat ≈ yhat.yhat

@test isapprox(sum(m_tul.beta[[5,25,45]]),  yhat_tul.yhat[end-3],atol=0.001)

@test isapprox(yhat_tul[yhat_tul.times.==0.5,:yhat], yhat_mul[yhat_mul.times.==0.5,:yhat],atol=0.001)
@test isapprox(yhat_tul[yhat_tul.times.==0.5,:yhat], [-2,1,2,5,6,9.],atol=0.0001)

if 1 == 0
    #for visualization
    using AlgebraOfGraphics
    yhat_mul.conditionA = categorical(yhat_mul.conditionA)
    m = mapping(:times,:yhat,color=:continuousA,linestyle=:conditionA)
    df = yhat_mul
    AlgebraOfGraphics.data(df) * visual(Lines) * m  |> draw
end