text:
"""
    naive_mccs(treelist)

Find sets of nodes which are:
- clades in all trees of `treelist`,
- all subclades of nodes are clades in all trees of `treelist`
  (both of these properties define consistency),
- maximal: adding a node to a set results it in not being a clade in at least
  one of the trees.

All the trees of `treelist` should share the same leaf nodes.
"""
function naive_mccs(treelist)
    # Check that trees share leaves
    sh = mapreduce(t->share_labels(t[1],t[2]), *, zip(treelist[1:end-1], treelist[2:end]))
    !sh && error("Can only be used on trees that share leaf nodes.")

    # List of splits in trees
    S = Tuple(SplitList(t) for t in treelist)

    # List of already visited nodes
    tref = treelist[1]
    checklist = Dict(k=>false for k in keys(tref.lleaves))

    # Explore one leaf at a time
    mc_clades = []
    for (cl,v) in checklist
        if !v # If leave not already visited
            # Root of current maximal clade, in all trees
            croot = [t.lleaves[cl] for t in treelist]
            # Initial individual, always a common clade in all trees since it's a leaf.
            clabel = [cl]

            # We're going to go up in all trees at the same time
            flag = true
            while flag &&  mapreduce(x->!x.isroot, *, croot; init=true) #prod(!x.isroot for x in croot)
                # Ancestors of current maximal clade in all trees
                nroot = [x.anc for x in croot]
                # Each element of `nroot` defines a set of labels corresponding to a tree.
                # There are two possibilites
                # (i) Those sets of labels match.
                ## In this case, we have a potential consistent clade.
                ## To check further, call `is_coherent_clade`.
                ## This is checked using the previously computed `SplitList`
                # (ii) Else,
                ## The topology of trees is inconsistent above `croot`.
                ## `croot` is an MCC, break.
                if mapreduce(
                        i->S[1].splitmap[nroot[1].label] == S[i].splitmap[nroot[i].label],
                        *,
                        2:length(nroot)
                    )
                    # --> `r ∈ nroot` is the same split in all trees
                    # check if children of `r` are also same splits in all trees
                    if is_coherent_clade(nroot,S)
                        if nroot == croot
                            # Singleton in the tree, or clade with a single node
                            # --> the algorithm is stuck on this node
                            croot = [x.anc for x in croot]
                        else
                            croot = nroot
                        end
                        clabel = S[1].leaves[S[1].splitmap[nroot[1].label].dat]
                    else
                        flag = false
                    end
                else
                    flag = false
                end
            end
            #
            map(x->checklist[x]=true, [c for c in clabel])
            push!(mc_clades, sort([c for c in clabel]))
        end
    end
    return sort(mc_clades, lt = clt)
end
naive_mccs(t...) = naive_mccs(t)

#= Custom order for MCCs =#
function clt(x,y)
    if length(x) < length(y)
        return true
    elseif length(x) > length(y)
        return false
    else
        return x[1] < y[1]
    end
end
function sort_mccs(mccs)
    return sort([sort(x) for x in mccs], lt=clt)
end

"""
    is_coherent_clade(roots::Array{<:TreeNode,1}, S::Tuple{<:SplitList,n} where n)

Do all children of `r` correspond to the same splits?

In practice, check that the splits `S[i].splitmap[c]` for `i` in `1:length(S)` and `c` in children of `r` are all the same.
"""
function is_coherent_clade(roots::Array{<:TreeNode,1}, S::Tuple{<:SplitList,n} where n)
    # All `r` in `roots` should at least have the same number of children
    if !mapreduce(ic-> length(roots[ic].child) == length(roots[1].child), *, 2:length(roots))
        return false
    end

    # Checking clade consistency
    for cref in roots[1].child
        if cref.isleaf
            for i in 2:length(roots)
                found = false
                for c in roots[i].child
                    if c.label == cref.label
                        found = true
                        break
                    end
                end
                if !found
                    return false
                end
            end
        else
            children = [cref]
            sref = S[1].splitmap[cref.label]
            # Looking at child `cref` for `roots[1]`
            # For every `r` in roots, there has to be a child `c` with the same split
            for i in 2:length(roots)
                for c in roots[i].child
                    if !c.isleaf && S[i].splitmap[c.label] == sref
                        push!(children, c)
                        break
                    end
                end
                if length(children) != i # nothing found
                    return false
                end
            end

            # Clade below all children found
            if !is_coherent_clade(children, S)
                return false
            end
        end
    end
    return true
end

"""
    name_mcc_clades!(treelist, MCC)

For each clade `m` in `MCC`:
- Rename the root `r` of `m` to `MCC_\$(i)` or (`\$(r.label)` if `r` is a leaf) where `i`
  is an integer starting at `label_init`.
- Rename each non-leaf internal node of `m` to `shared_\$i_\$j` where `j` is an index
  specific to `m`.

"""
function name_mcc_clades!(treelist, MCC)
    # Finding initial label
    label_init = 1
    for t in treelist
        for n in values(t.lnodes)
            if match(r"MCC", n.label)!=nothing && parse(Int, n.label[5:end]) >= label_init
                label_init = parse(Int, n.label[5:end]) + 1
            end
        end
    end

    nd = Dict()
    for (i,m) in enumerate(MCC)
        cl = i + label_init - 1
        # Renaming root
        for t in treelist
            r = lca(t, m)
            old_label = r.label
            new_label = r.isleaf ? "$(old_label)" : "MCC_$(cl)"
            r.label = new_label
            delete!(t.lnodes, old_label)
            t.lnodes[new_label] = r
            nd[new_label] = m
        end
    end

    return nd
end

"""
    reduce_to_mcc(tree, MCC)

Reduce `tree` to its MCC by grouping leaves. Returns a tree with `length(MCC)` leaves.
"""
function reduce_to_mcc(tree::Tree, MCC)
    out = copy(tree, TreeTools.EmptyData)
    reduce_to_mcc!(out, MCC)
    return out
end
"""
    reduce_to_mcc!(tree, MCC)

Reduce `tree` to `MCCs` by grouping leaves.
"""
function reduce_to_mcc!(tree::Tree, MCC)
    for m in MCC
        r = lca(tree, m)
        if r.isroot
            node2tree!(tree, TreeNode(
                    r.data;
                    isleaf=true, isroot = true, label=r.label, r.tau
                ))
        elseif !r.isleaf
        	for c in reverse(r.child)
        		# prunenode!(c)
        		prunesubtree!(tree, c; remove_singletons=false)
        	end
        	r.isleaf = true
        	r.isroot = false
        	tree.lleaves[r.label] = r
        end
    end
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

struct TypeRegistry{T}
    by_name::IdDict{Symbol,Type{<:T}}
    by_type::IdDict{Type{<:T},Symbol}
end

function TypeRegistry{T}() where T
    TypeRegistry{T}(
        IdDict{Symbol,Type{<:T}}(),
        IdDict{Type{<:T},Symbol}()
    )
end

@noinline Base.getindex(registry::TypeRegistry{T}, name::Symbol) where T =
    registry.by_name[name]

@noinline Base.getindex(registry::TypeRegistry{T}, @nospecialize(U::Type{<:T})) where T =
    registry.by_type[U]

@noinline function Base.setindex!(registry::TypeRegistry{T}, @nospecialize(U::Type{<:T}), name::Symbol) where T
    registry.by_name[name] = U
    registry.by_type[U] = name
end

ext: jl


text:
export catch_error

import Base: show

"""
    catch_error(selectorFn::F) where F

Creates a `CatchErrorOperator`, which catches errors on the observable to be handled by returning a new observable or throwing an error.

# Arguments:
- `selectorFn::F`: a callable object that takes as arguments err, which is the error, and caught, which is the source observable, in case you'd like to "retry" that observable by returning it again. Whatever observable is returned by the selector will be used to continue the observable chain.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from(1:5) |> safe() |> map(Int, (d) -> d == 4 ? error(4) : d) |> catch_error((err, obs) -> of(1))

subscribe!(source, logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Data: 1
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`rerun`](@ref), [`logger`](@ref), [`safe`](@ref)
"""
catch_error(selectorFn::F) where F = CatchErrorOperator{F}(selectorFn)

struct CatchErrorOperator{F} <: InferableOperator
    selectorFn :: F
end

function on_call!(::Type{L}, ::Type{L}, operator::CatchErrorOperator{F}, source) where { L, F }
    return proxy(L, source, CatchErrorProxy{F}(operator.selectorFn))
end

operator_right(operator::CatchErrorOperator, ::Type{L}) where L = L

struct CatchErrorProxy{F} <: ActorSourceProxy
    selectorFn :: F
end

actor_proxy!(::Type{L}, proxy::CatchErrorProxy{F},  actor::A)  where { L, A, F } = CatchErrorActor{L, A, F}(proxy.selectorFn, actor, false, nothing, nothing)
source_proxy!(::Type{L}, proxy::CatchErrorProxy{F}, source::S) where { L, S, F } = CatchErrorSource{L, S}(source)

mutable struct CatchErrorActor{L, A, F} <: Actor{L}
    selectorFn           :: F
    actor                :: A
    is_completed         :: Bool
    current_source       :: Union{Nothing, Any}
    current_subscription :: Union{Nothing, Teardown}
end

function on_next!(actor::CatchErrorActor{L}, data::L) where L
    if !actor.is_completed
        next!(actor.actor, data)
    end
end

function on_error!(actor::CatchErrorActor, err)
    if !actor.is_completed
        if actor.current_subscription !== nothing
            unsubscribe!(actor.current_subscription)
        end

        fallback_source = actor.selectorFn(err, actor.current_source)

        actor.current_source       = fallback_source
        actor.current_subscription = subscribe!(fallback_source, actor)
    end
end

function on_complete!(actor::CatchErrorActor)
    if !actor.is_completed
        actor.is_completed = true
        complete!(actor.actor)
    end
end

@subscribable struct CatchErrorSource{L, S} <: Subscribable{L}
    source :: S
end

struct CatchErrorSubscription <: Teardown
    catch_error_actor
end

as_teardown(::Type{<:CatchErrorSubscription}) = UnsubscribableTeardownLogic()

function on_subscribe!(source::CatchErrorSource, actor::CatchErrorActor)
    actor.current_source       = source.source
    actor.current_subscription = subscribe!(source.source, actor)

    return CatchErrorSubscription(actor)
end

function on_unsubscribe!(subscription::CatchErrorSubscription)
    current_subscription = subscription.catch_error_actor.current_subscription

    subscription.catch_error_actor.current_source       = nothing
    subscription.catch_error_actor.current_subscription = nothing

    return unsubscribe!(current_subscription)
end

Base.show(io::IO, ::CatchErrorOperator)             = print(io, "CatchErrorOperator()")
Base.show(io::IO, ::CatchErrorProxy)                = print(io, "CatchErrorProxy()")
Base.show(io::IO, ::CatchErrorActor{L})     where L = print(io, "CatchErrorActor($L)")
Base.show(io::IO, ::CatchErrorSource{L})    where L = print(io, "CatchErrorSource($L)")
Base.show(io::IO, ::CatchErrorSubscription)         = print(io, "CatchErrorSubscription()")


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SelfConvolutive

using Nemo

export ModuleSelfConvolutive, SelfConvRec
export L000698, L001710, L003319, L005411, L005412, L006012, L006318
export L047891, L062980, L082298, L082301, L082302, L105523, L107716
export L111529, L111530, L111531, L111532, L111533, L146559, L167872

"""

Some exactly solvable self-convolutive recurrences.

* SelfConvRec, L000698, L001710, L003319, L005411, L005412, L006012, L006318, L047891, L062980, L082298, L082301, L082302, L105523, L107716, L111529, L111530, L111531, L111532, L111533, L146559, L167872
"""
const ModuleSelfConvolutive = ""

# R. J. Martin and M. J. Kearney,
# An exactly solvable self-convolutive recurrence
# http://arXiv.org/abs/1103.4936
# Integral representation of certain combinatorial recurrences
# see also http://oeis.org/A172455

Convolution(A, n) = sum(A[j+1] * A[n-j+1] for j ∈ 0:n)

"""

An exactly solvable self-convolutive recurrence studied by R. J. Martin and M. J. Kearney.
"""
function SelfConvRec(len, a, b, c)
    A = fill(fmpz(0), len)
    len == 0 && return A
    A[1] = 1

    for n ∈ 1:len-1
        A[n+1] = (a * (n + 1) + b) * A[n] + c * Convolution(A, n - 1)
    end
    A
end

# """
# Return the powers of 2.
# """
# L000079(len::Int) = SelfConvRec(len, 1, 1, -1)
# """
# Return the Catalan numbers.
# """
# L000108(len::Int) = SelfConvRec(len, 0, 0, 1)
# """
# Return the factorial numbers.
# """
# L000142(len::Int) = SelfConvRec(len, 1, -1, 0)
# """
# Return the powers of 3.
# """
# L000244(len::Int) = SelfConvRec(len, 2, 1, -2)
# """
# Return the powers of 5.
# """
# L000351(len::Int) = SelfConvRec(len, 4, 1, -4)
# """
# Return the powers of 6.
# """
# L000400(len::Int) = SelfConvRec(len, 5, 1, -5)
# """
# Return the powers of 7.
# """
# L000420(len::Int) = SelfConvRec(len, 6, 1, -6)

"""

Return the number of indecomposable perfect matchings on ``[2n]``.
"""
L000698(len::Int) = SelfConvRec(len, 2, -3, 1)
"""

Return the order of alternating group ``A_n``, or number of even permutations of ``n`` letters.
"""
L001710(len::Int) = SelfConvRec(len, 1, 1, 0)

# """
# Return n!/3!.
# """
# L001715(len::Int) = SelfConvRec(len, 1, 2, 0)
# """
# Return n!/4!.
# """
# L001720(len::Int) = SelfConvRec(len, 1, 3, 0)
# """
# Return n!/5!.
# """
# L001725(len::Int) = SelfConvRec(len, 1, 4, 0)
# """
# Return n!/6!.
# """
# L001730(len::Int) = SelfConvRec(len, 1, 5, 0)

"""

Return the number of connected permutations of ``[1..n]``. Also called indecomposable permutations.
"""
L003319(len::Int) = SelfConvRec(len, 1, -2, 1)
"""

Return the number of non-vanishing Feynman diagrams of order ``2n`` for the electron or the photon propagators in quantum electrodynamics.
"""
L005411(len::Int) = SelfConvRec(len, 2, -4, 1)
"""

Return the number of non-vanishing Feynman diagrams of order ``2n`` for the vacuum polarization in quantum electrodynamics.
"""
L005412(len::Int) = SelfConvRec(len, 2, -2, 1)
"""

Counting some sets of permutations.
"""
L006012(len::Int) = SelfConvRec(len, -1, 2, 2)
"""

Return the large Schröder numbers.
"""
L006318(len::Int) = SelfConvRec(len, 0, 1, 1)
"""

Return the number of planar rooted trees with ``n`` nodes and tricolored end nodes.
"""
L047891(len::Int) = SelfConvRec(len, 0, 2, 1)
# """
# Return (n+7)!/7!.
# """
# L049388(len::Int) = SelfConvRec(len, 1, 6, 0)
# """
# Return (3*n+4)!!!/4!!!.
# """
# L051604(len::Int) = SelfConvRec(len, 3, 1, 0)
# """
# Return (3*n+5)!!!/5!!!.
# """
# L051605(len::Int) = SelfConvRec(len, 3, 2, 0)
# """
# Return (3*n+6)!!!/6!!!.
# """
# L051606(len::Int) = SelfConvRec(len, 3, 3, 0)
# """
# Return (3*n+7)!!!/7!!!.
# """
# L051607(len::Int) = SelfConvRec(len, 3, 4, 0)
# """
# Return (3*n+8)!!!/8!!!.
# """
# L051608(len::Int) = SelfConvRec(len, 3, 5, 0)
# """
# Return (3*n+9)!!!/9!!!.
# """
# L051609(len::Int) = SelfConvRec(len, 3, 6, 0)
# """
# Return (4*n+5)(!^4)/5(!^4).
# """
# L051617(len::Int) = SelfConvRec(len, 4, 1, 0)
# """
# Return  (4*n+6)(!^4)/6(!^4).
# """
# L051618(len::Int) = SelfConvRec(len, 4, 2, 0)
# """
# Return (4*n+7)(!^4)/7(!^4).
# """
# L051619(len::Int) = SelfConvRec(len, 4, 3, 0)
# """
# Return (4*n+8)(!^4)/8(!^4).
# """
# L051620(len::Int) = SelfConvRec(len, 4, 4, 0)
# """
# Return (4*n+9)(!^4)/9(!^4).
# """
# L051621(len::Int) = SelfConvRec(len, 4, 5, 0)
# """
# Return (4*n+10)(!^4)/10(!^4).
# """
# L051622(len::Int) = SelfConvRec(len, 4, 6, 0)
# """
# Return (5*n+6)(!^5)/6.
# """
# L051687(len::Int) = SelfConvRec(len, 5, 1, 0)
# """
# Return (5*n+7)(!^5)/7(!^5).
# """
# L051688(len::Int) = SelfConvRec(len, 5, 2, 0)
# """
# Return (5*n+8)(!^5)/8(!^5).
# """
# L051689(len::Int) = SelfConvRec(len, 5, 3, 0)
# """
# Return (5*n+9)(!^5)/9(!^5).
# """
# L051690(len::Int) = SelfConvRec(len, 5, 4, 0)
# """
# Return (5*n+10)(!^5)/10(!^5).
# """
# L051691(len::Int) = SelfConvRec(len, 5, 5, 0)
# """
# Return ((6*n+7)(!^6))/7.
# """
# L053100(len::Int) = SelfConvRec(len, 6, 1, 0)
# """
# Return ((6*n+8)(!^6))/8(!^6).
# """
# L053101(len::Int) = SelfConvRec(len, 6, 2, 0)
# """
# Return ((6*n+9)(!^6))/9(!^6).
# """
# L053102(len::Int) = SelfConvRec(len, 6, 3, 0)
# """
# Return ((6*n+10)(!^6))/10(!^6).
# """
# L053103(len::Int) = SelfConvRec(len, 6, 4, 0)
# """
# Return ((7*n+8)(!^7))/8.
# """
# L053104(len::Int) = SelfConvRec(len, 7, 1, 0)
# """
# Return ((7*n+9)(!^7))/9(!^7).
# """
# L053105(len::Int) = SelfConvRec(len, 7, 2, 0)
# """
# Return ((7*n+10)(!^7))/10(1^7).
# """
# L053106(len::Int) = SelfConvRec(len, 7, 3, 0)

"""

Return the number of rooted unlabeled connected triangular maps on a compact closed oriented surface with ``2n`` faces.
"""
L062980(len::Int) = SelfConvRec(len, 6, -8, 1)
"""

Return the number of lattice paths from ``(0,0)`` to ``(n+1,n+1)`` that consist of steps ``(i,0)`` and ``(0,j)`` with ``i,j≥1`` and stay strictly below the diagonal line ``y=x`` except at the endpoints.
"""
L082298(len::Int) = SelfConvRec(len, 0, 3, 1)
"""

Return the number of Schröder paths of semilength n in which the (2,0)-steps come in 4 colors.
"""
L082301(len::Int) = SelfConvRec(len, 0, 4, 1)
"""

Return the coefficients in the expansion of ``(1-5x-√(25x^2-14x+1))/(2x)``.
"""
L082302(len::Int) = SelfConvRec(len, 0, 5, 1)

# """
# Return the coefficients in the expansion of ``(1-6x-√(36x^2-16x+1))/(2x)``.
# """
# L082305(len::Int) = SelfConvRec(len, 0, 6, 1)
# """
# Return the coefficients in the expansion of ``(1-7x-√(49x^2-18x+1))/(2x)``.
# """
# L082366(len::Int) = SelfConvRec(len, 0, 7, 1)
# """
# Return the coefficients in the expansion of ``(1-8x-√(64x^2-20x+1))/(2*x)``.
# """
# L082367(len::Int) = SelfConvRec(len, 0, 8, 1)

"""

Return the coefficients in the expansion of ``(1+2x-√(1+4x^2))/(2x)``.
"""
L105523(len::Int) = SelfConvRec(len, 0, -2, 1)
"""

A transform the of triple factorial numbers.
"""
L107716(len::Int) = SelfConvRec(len, 3, -4, 1)
"""

A convolutory inverse of the factorial sequence.
"""
L111529(len::Int) = SelfConvRec(len, 1, -3, 2)
"""

A convolutory inverse of the factorial sequence.
"""
L111530(len::Int) = SelfConvRec(len, 1, -4, 3)
"""

A convolutory inverse of the factorial sequence.
"""
L111531(len::Int) = SelfConvRec(len, 1, -5, 4)
"""

A convolutory inverse of the factorial sequence.
"""
L111532(len::Int) = SelfConvRec(len, 1, -6, 5)
"""

A convolutory inverse of the factorial sequence.
"""
L111533(len::Int) = SelfConvRec(len, 1, -7, 6)

#"""
#"""
#L111546(len::Int) = SelfConvRec(len, 1, 0, 1)
#"""
#"""
#L111556(len::Int) = SelfConvRec(len, 1, 1, 1)

"""

Return the expansion of ``(1-x)/(1 - 2x + 2x^2)``.
"""
L146559(len::Int) = SelfConvRec(len, 1, 1, -2)
"""

Return half the number of Feynman diagrams of order ``2(n+1)``, for the electron self-energy in quantum electrodynamics.
"""
L167872(len::Int) = SelfConvRec(len, 2, -3, 2)

# """
# """
# L172450(len::Int) = SelfConvRec(len, 2, 0, -1)
# """
# """
# L172455(len::Int) = SelfConvRec(len, 6, -4, -1)
# """
# """
# L172485(len::Int) = SelfConvRec(len, -1, -2, 3)
# """
# """
# L177354(len::Int) = SelfConvRec(len, 1, 2, 1)

#START-TEST-########################################################
using Test, SeqTests

function test()
    @testset "SelfConvRec" begin
        if data_installed()

            L = [
                L005411,
                L005412,
                L006012,
                L006318,
                L047891,
                L062980,
                L082298,
                L082301,
                L082302,
                L105523,
                L107716,
                L111529,
                L111530,
                L111531,
                L111532,
                L111533,
                L146559,
                L167872
            ]
            SeqTest(L, 'L')
        end
    end
end

function demo()

    for n ∈ 0:6
        println(n, " ↦ ", L005411(n))
    end

    f(len) = SelfConvRec(len, 2, -1, -1)
    println(f(11))
end

"""

L005411(500) :: 0.238653 seconds (252.99 k allocations: 3.887 MB)
"""
function perf()
    GC.gc()
    @time L005411(500)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
L --> A000079
fmpz[1, 2, 4, 8, 16, 32]
fmpz[1, 2, 4, 8, 16, 32]
L --> A000108
fmpz[1, 1, 2, 5, 14, 42]
fmpz[1, 1, 2, 5, 14, 42]
L --> A000142
fmpz[1, 1, 2, 6, 24, 120]
fmpz[1, 1, 2, 6, 24, 120]
L --> A000244
fmpz[1, 3, 9, 27, 81, 243]
fmpz[1, 3, 9, 27, 81, 243]
L --> A000351
fmpz[1, 5, 25, 125, 625, 3125]
fmpz[1, 5, 25, 125, 625, 3125]
L --> A000400
fmpz[1, 6, 36, 216, 1296, 7776]
fmpz[1, 6, 36, 216, 1296, 7776]
L --> A000420
fmpz[1, 7, 49, 343, 2401, 16807]
fmpz[1, 7, 49, 343, 2401, 16807]
#############################
L --> A000698
fmpz[1, 1, 2, 10, 74, 706]
fmpz[1, 2, 10, 74, 706, 8162]
#####################
seq((n+2)!/2, n=0..9);
L --> A001710
fmpz[1, 1, 1, 3, 12, 60]
fmpz[1, 3, 12, 60, 360, 2520]
######################
L --> A001715
fmpz[1, 4, 20, 120, 840, 6720]
fmpz[1, 4, 20, 120, 840, 6720]
L --> A001720
fmpz[1, 5, 30, 210, 1680, 15120]
fmpz[1, 5, 30, 210, 1680, 15120]
L --> A001725
fmpz[1, 6, 42, 336, 3024, 30240]
fmpz[1, 6, 42, 336, 3024, 30240]
L --> A001730
fmpz[1, 7, 56, 504, 5040, 55440]
fmpz[1, 7, 56, 504, 5040, 55440]
##################
L --> A003319
fmpz[0, 1, 1, 3, 13, 71]
fmpz[1, 1, 3, 13, 71, 461]
##################
L --> A005411
fmpz[1, 1, 4, 25, 208, 2146]
fmpz[1, 1, 4, 25, 208, 2146]
L --> A005412
fmpz[1, 3, 18, 153, 1638, 20898]
fmpz[1, 3, 18, 153, 1638, 20898]
L --> A006012
fmpz[1, 2, 6, 20, 68, 232]
fmpz[1, 2, 6, 20, 68, 232]
L --> A006318
fmpz[1, 2, 6, 22, 90, 394]
fmpz[1, 2, 6, 22, 90, 394]
L --> A047891
fmpz[1, 3, 12, 57, 300, 1686]
fmpz[1, 3, 12, 57, 300, 1686]
L --> A049388
fmpz[1, 8, 72, 720, 7920, 95040]
fmpz[1, 8, 72, 720, 7920, 95040]
L --> A051604
fmpz[1, 7, 70, 910, 14560, 276640]
fmpz[1, 7, 70, 910, 14560, 276640]
L --> A051605
fmpz[1, 8, 88, 1232, 20944, 418880]
fmpz[1, 8, 88, 1232, 20944, 418880]
L --> A051606
fmpz[1, 9, 108, 1620, 29160, 612360]
fmpz[1, 9, 108, 1620, 29160, 612360]
L --> A051607
fmpz[1, 10, 130, 2080, 39520, 869440]
fmpz[1, 10, 130, 2080, 39520, 869440]
L --> A051608
fmpz[1, 11, 154, 2618, 52360, 1204280]
fmpz[1, 11, 154, 2618, 52360, 1204280]
L --> A051609
fmpz[1, 12, 180, 3240, 68040, 1632960]
fmpz[1, 12, 180, 3240, 68040, 1632960]
L --> A051617
fmpz[1, 9, 117, 1989, 41769, 1044225]
fmpz[1, 9, 117, 1989, 41769, 1044225]
L --> A051618
fmpz[1, 10, 140, 2520, 55440, 1441440]
fmpz[1, 10, 140, 2520, 55440, 1441440]
L --> A051619
fmpz[1, 11, 165, 3135, 72105, 1946835]
fmpz[1, 11, 165, 3135, 72105, 1946835]
L --> A051620
fmpz[1, 12, 192, 3840, 92160, 2580480]
fmpz[1, 12, 192, 3840, 92160, 2580480]
L --> A051621
fmpz[1, 13, 221, 4641, 116025, 3364725]
fmpz[1, 13, 221, 4641, 116025, 3364725]
L --> A051622
fmpz[1, 14, 252, 5544, 144144, 4324320]
fmpz[1, 14, 252, 5544, 144144, 4324320]
L --> A051687
fmpz[1, 11, 176, 3696, 96096, 2978976]
fmpz[1, 11, 176, 3696, 96096, 2978976]
L --> A051688
fmpz[1, 12, 204, 4488, 121176, 3877632]
fmpz[1, 12, 204, 4488, 121176, 3877632]
L --> A051689
fmpz[1, 13, 234, 5382, 150696, 4972968]
fmpz[1, 13, 234, 5382, 150696, 4972968]
L --> A051690
fmpz[1, 14, 266, 6384, 185136, 6294624]
fmpz[1, 14, 266, 6384, 185136, 6294624]
L --> A051691
fmpz[1, 15, 300, 7500, 225000, 7875000]
fmpz[1, 15, 300, 7500, 225000, 7875000]
L --> A053100
fmpz[1, 13, 247, 6175, 191425, 7082725]
fmpz[1, 13, 247, 6175, 191425, 7082725]
L --> A053101
fmpz[1, 14, 280, 7280, 232960, 8852480]
fmpz[1, 14, 280, 7280, 232960, 8852480]
L --> A053102
fmpz[1, 15, 315, 8505, 280665, 10945935]
fmpz[1, 15, 315, 8505, 280665, 10945935]
L --> A053103
fmpz[1, 16, 352, 9856, 335104, 13404160]
fmpz[1, 16, 352, 9856, 335104, 13404160]
L --> A053104
fmpz[1, 15, 330, 9570, 344520, 14814360]
fmpz[1, 15, 330, 9570, 344520, 14814360]
L --> A053105
fmpz[1, 16, 368, 11040, 408480, 17973120]
fmpz[1, 16, 368, 11040, 408480, 17973120]
L --> A053106
fmpz[1, 17, 408, 12648, 480624, 21628080]
fmpz[1, 17, 408, 12648, 480624, 21628080]
L --> A062980
fmpz[1, 5, 60, 1105, 27120, 828250]
fmpz[1, 5, 60, 1105, 27120, 828250]
L --> A082298
fmpz[1, 4, 20, 116, 740, 5028]
fmpz[1, 4, 20, 116, 740, 5028]
L --> A082301
fmpz[1, 5, 30, 205, 1530, 12130]
fmpz[1, 5, 30, 205, 1530, 12130]
L --> A082302
fmpz[1, 6, 42, 330, 2814, 25422]
fmpz[1, 6, 42, 330, 2814, 25422]
L --> A082305
fmpz[1, 7, 56, 497, 4760, 48174]
fmpz[1, 7, 56, 497, 4760, 48174]
L --> A082366
fmpz[1, 8, 72, 712, 7560, 84616]
fmpz[1, 8, 72, 712, 7560, 84616]
L --> A082367
fmpz[1, 9, 90, 981, 11430, 140058]
fmpz[1, 9, 90, 981, 11430, 140058]
L --> A105523
fmpz[1, -1, 0, 1, 0, -2]
fmpz[1, -1, 0, 1, 0, -2]
L --> A107716
fmpz[1, 3, 21, 219, 2973, 49323]
fmpz[1, 3, 21, 219, 2973, 49323]
L --> A111529
fmpz[1, 1, 4, 22, 148, 1156]
fmpz[1, 1, 4, 22, 148, 1156]
L --> A111530
fmpz[1, 1, 5, 33, 261, 2361]
fmpz[1, 1, 5, 33, 261, 2361]
L --> A111531
fmpz[1, 1, 6, 46, 416, 4256]
fmpz[1, 1, 6, 46, 416, 4256]
L --> A111532
fmpz[1, 1, 7, 61, 619, 7045]
fmpz[1, 1, 7, 61, 619, 7045]
L --> A111533
fmpz[1, 1, 8, 78, 876, 10956]
fmpz[1, 1, 8, 78, 876, 10956]
L --> A111546
fmpz[1, 3, 15, 99, 783, 7083]
fmpz[1, 3, 15, 99, 783, 7083]
L --> A111556
fmpz[1, 4, 24, 184, 1664, 17024]
fmpz[1, 4, 24, 184, 1664, 17024]
L --> A146559
fmpz[1, 1, 0, -2, -4, -4]
fmpz[1, 1, 0, -2, -4, -4]
L --> A167872
fmpz[1, 3, 21, 207, 2529, 36243]
fmpz[1, 3, 21, 207, 2529, 36243]
L --> A172450
fmpz[1, 3, 12, 63, 432, 3798]
fmpz[1, 3, 12, 63, 432, 3798]
L --> A172485
fmpz[1, -1, -1, 3, 3, -21]
fmpz[1, -1, -1, 3, 3, -21]
#################################
L --> A177354
fmpz[5, 35, 305, 3095, 35225, 439775]
fmpz[1, 5, 35, 305, 3095, 35225]
###########################

L = [ L000079, L000108, L000142, L000244, L000351, L000400,
 L000420, L001715, L001720, L001725, L001730,
 L005411, L005412, L006012, L006318, L047891, L049388, L051604,
 L051605, L051606, L051607, L051608, L051609, L051617, L051618,
 L051619, L051620, L051621, L051622, L051687, L051688, L051689,
 L051690, L051691, L053100, L053101, L053102, L053103, L053104,
 L053105, L053106, L062980, L082298, L082301, L082302, L082305,
 L082366, L082367, L105523, L107716, L111529, L111530, L111531,
 L111532, L111533, L111546, L111556, L146559, L167872,
 L172450, L172485 ]
# L000698, L001710, L003319, L177354
SeqTest(L, 'L')

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
@testset "kernels.jl" begin
    kernel = TVExponentialKernel()

    # traits
    @test KernelFunctions.metric(kernel) === TotalVariation()

    # simple evaluation
    x, y = rand(10), rand(10)
    @test kernel(x, y) == exp(-totalvariation(x, y))

    # transformations
    @test (kernel ∘ ScaleTransform(0.1))(x, y) == exp(-0.1 * totalvariation(x, y))
    ard = rand(10)
    @test (kernel ∘ ARDTransform(ard))(x, y) == exp(-totalvariation(ard .* x, ard .* y))
end


text:
@testset "Kernels" begin
    import LikelihoodfreeInference: bandwidth, update!, Bandwidth, RosenblattSchedule, Scheduled
    k = Kernel(type = Gaussian,
               bandwidth = Bandwidth(heuristic = nothing, gamma = 2.))
    @test bandwidth(k) == 2.
    update!(k, 0)
    @test bandwidth(k) == 2.
    @test k(1, 4) == exp(-(4 - 1)^2/(2 * 2^2))
    @test k([2, 0, 1], zeros(3)) ≈ exp(-5/(2*2^2))
    k = Kernel(type = Gaussian,
               bandwidth = Bandwidth(heuristic = nothing, gamma = [2., 6.]))
    @test bandwidth(k) == [2., 6.]
    @test k([3, 0], [4, 2]) == exp(-sum(abs2, ([3, 0] - [4, 2]) ./ [2., 6])/2)
    k = Kernel(type = ModifiedGaussian,
               bandwidth = Bandwidth(heuristic = nothing, gamma = .5))
    @test k(1, 0) == exp(-1/(2*.5))
    @test k(.3, 0) == exp(-.3^2/(2*.5^2))
    k = Kernel(bandwidth = Bandwidth(heuristic = MedianHeuristic()))
    update!(k, [2, 2])
    @test bandwidth(k) == 2
    k = Kernel(bandwidth = Bandwidth(heuristic = MedianHeuristic(.1)))
    update!(k, [2, 2])
    @test bandwidth(k) == .1 * 2
    k = Kernel(bandwidth = Bandwidth(heuristic = ScottsHeuristic()))
    update!(k, [2, 2])
    @test bandwidth(k) == 0
    k = Kernel(bandwidth = Scheduled(Bandwidth(heuristic = nothing, gamma = 2),
                                     RosenblattSchedule(2)))
    update!(k, nothing)
    update!(k, nothing)
    @test bandwidth(k) == 2 * 2^(-1/(2 + 4))
end


text:
# Autogenerated wrapper script for Exodus_jll for x86_64-linux-gnu
export libexodus

using Zlib_jll
using NetCDF_jll
using HDF5_jll
JLLWrappers.@generate_wrapper_header("Exodus")
JLLWrappers.@declare_library_product(libexodus, "libexodus.so.2")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, NetCDF_jll, HDF5_jll)
    JLLWrappers.@init_library_product(
        libexodus,
        "lib/libexodus.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Oceananigans.TimeSteppers: calculate_pressure_correction!, pressure_correct_velocities!

calculate_pressure_correction!(::HydrostaticFreeSurfaceModel, Δt) = nothing

#####
##### Barotropic pressure correction for models with a free surface
#####

const HFSM = HydrostaticFreeSurfaceModel
const ExplicitFreeSurfaceHFSM = HFSM{<:Any, <:Any, <:Any, <:ExplicitFreeSurface}
const ImplicitFreeSurfaceHFSM = HFSM{<:Any, <:Any, <:Any, <:ImplicitFreeSurface}

pressure_correct_velocities!(model::ExplicitFreeSurfaceHFSM, Δt; kwargs...) = nothing

#####
##### Barotropic pressure correction for models with a free surface
#####

function pressure_correct_velocities!(model::ImplicitFreeSurfaceHFSM, Δt;
                                      dependencies = device_event(model.architecture))

    event = launch!(model.architecture, model.grid, :xyz,
                    _barotropic_pressure_correction,
                    model.velocities,
                    model.grid,
                    Δt,
                    model.free_surface.gravitational_acceleration,
                    model.free_surface.η,
                    dependencies = dependencies)

    wait(device(model.architecture), event)

    return nothing
end

@kernel function _barotropic_pressure_correction(U, grid, Δt, g, η)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.u[i, j, k] -= g * Δt * ∂xᶠᶜᵃ(i, j, k, grid, η)
        U.v[i, j, k] -= g * Δt * ∂yᶜᶠᵃ(i, j, k, grid, η)
    end
end


text:
"""
General electronic course-graining platform
"""
module FermiCG

#####################################
# External packages
#
using Compat
using HDF5
using KrylovKit
using LinearAlgebra
using NDTensors
using PackageCompiler
using Parameters
using Printf
using TimerOutputs
using BenchmarkTools
using OrderedCollections
using IterTools
# using Unicode
#
#####################################

#####################################
# Local Imports
#
include("Utils.jl")
include("hosvd.jl")
include("StringCI/StringCI.jl");
include("SymDenseMats.jl");
include("Solvers.jl");
include("Hamiltonians.jl");
include("Clusters.jl")
include("Indexing.jl")
include("States.jl")

include("FockSparse_ElementSparse.jl")
include("FockSparse_BlockSparse.jl")
include("FockSparse_BlockSparseTucker.jl")

include("ClusteredTerms.jl")
#include("ClusteredStates.jl")

include("tucker_inner.jl")
include("tucker_outer.jl")
include("bst.jl")

include("tpsci_inner.jl")
include("tpsci_matvec_thread.jl")
include("tpsci_outer.jl")
include("tpsci_helpers.jl")

include("dense_inner.jl")
include("dense_outer.jl")

include("CompressedTucker.jl")
include("CMFs.jl")
include("pyscf/PyscfFunctions.jl");
#
#####################################

export StringCI
export InCoreInts
export Molecule
export Atom
export Cluster
export ClusterBasis
export ClusteredState
export ClusterConfig
export FockConfig

end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# Autogenerated wrapper script for eccodes_jll for i686-w64-mingw32
export eccodes

using libpng_jll
using OpenJpeg_jll
JLLWrappers.@generate_wrapper_header("eccodes")
JLLWrappers.@declare_library_product(eccodes, "libeccodes.dll")
function __init__()
    JLLWrappers.@generate_init_header(libpng_jll, OpenJpeg_jll)
    JLLWrappers.@init_library_product(
        eccodes,
        "bin\\libeccodes.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function varsassign(mystr::String, varname::String, n::Int, ustr, flg::Int)

if flg==0
    for i = 1:n
        str1 = varname * string(i);
        if contains(ustr, str1)
            str2 = varname * "[" * string(i-1) * "]";
            mystr = mystr * "\t\tT " * str1 * " = " * str2 * ";\n";
        end
    end
elseif flg==1
    for i = 1:n
        str1 = varname * string(i);
        if contains(ustr, str1)
            str2 = varname * "[" * string(i-1) * "]";
            mystr = mystr + "\t\tint " * str1 * " = " + str2 * ";\n";
        end
    end
elseif flg==2
    for i = 1:n
        str1 = varname * string(i);
        if contains(ustr, str1)
            str2 = varname * "[" * string(i-1) * " + i*" * string(n) * "]";
            mystr = mystr * "\t\tT " * str1 * " = " * str2 * ";\n";
        end
    end
elseif flg==3
    for i = 1:n
        str1 = varname * string(i);
        if contains(ustr, str1)
            str2 = varname * "[" * string(i-1) * " + i*" * string(n) * "]";
            mystr = mystr * "\t\tint " * str1 * " = " * str2 * ";\n";
        end
    end
end

return mystr;

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Documenter, BitBasis
import BitBasis

const PAGES = [
    "Home" => "index.md",
    "Tutorial" => "tutorial.md",
    "Manual" => "man.md",
]

makedocs(
    modules = [BitBasis],
    format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = ("deploy" in ARGS) ? "https://quantumbfs.github.io/BitBasis.jl/latest/" : nothing,
    ),
    assets = ["assets/favicon.ico"],
    clean = false,
    sitename = "BitBasis.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = PAGES
)

deploydocs(
    repo = "github.com/QuantumBFS/BitBasis.jl.git",
    target = "build",
)


text:
# ---
# title: 112. Path Sum
# id: problem112
# author: zhwang
# date: 2022-02-21
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/path-sum/description/>
# hidden: true
# ---
#
# Given a binary tree and a sum, determine if the tree has a root-to-leaf path
# such that adding up all the values along the path equals the given sum.
#
# **Note:**  A leaf is a node with no children.
#
# **Example:**
#
# Given the below binary tree and `sum = 22`,
#
#
#
#           **5**
#          **/** \
#         **4**   8
#        **/**   / \
#       **11**  13  4
#      /  **\**      \
#     7    **2**      1
#
#
# return true, as there exist a root-to-leaf path `5->4->11->2` which sum is 22.
#
#
## @lc code=start
using LeetCode

has_path_sum(::Nothing, ::Int) = false
function has_path_sum(root::TreeNode, target_sum::Int)::Bool
    root.val == target_sum && isnothing(root.left) && isnothing(root.right) && return true
    return any(has_path_sum.([root.left, root.right], target_sum - root.val))
end

## @lc code=end


text:
# Collection of functions and signals for measuring selections

using DataFrames

function avg_gray(selection)
  # returns avg gray value across array
  sum = 0
  for i in 1:length(selection)
    sum += Float64(selection[i])
  end
  sum = sum/length(selection)
  return ColorTypes.Gray{FixedPointNumbers.Normed{UInt8,8}}(sum).val
end

# Initializes a measurements table
function init_measure_tbl()
  return DataFrame(Area = @data([]),Mean = @data([]),Minimum = @data([]),Maximum = @data([]))
end

# Returns a row that can be pushed to a measurements table
function measure(selection) # also saves selection for labeling purposes?
  return @data([length(selection),avg_gray(selection),minimum(selection).val,maximum(selection).val])
end

#=
using TestImages, Gtk.ShortNames

img = testimage("cameraman.tif")

rectview = Signal(view(img,1:20,1:40))

dataf = map(rectview) do r
    df = init_measure_tbl()
    push!(df, measure(img))
    push!(df, measure(r))
    df
end

txt = TextView()
push!(txt,string(value(dataf))) # should update when dataf updates...
# monospace font?

win2 = Window("Measure", 300, 200)
nb = Notebook()
push!(nb,txt,"Name")
push!(win2,nb)
showall(win2)
nothing
=#


text:
# ---
# title: 1414. Find the Minimum Number of Fibonacci Numbers Whose Sum Is K
# id: problem1414
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/find-the-minimum-number-of-fibonacci-numbers-whose-sum-is-k/description/>
# hidden: true
# ---
#
# Given an integer `k`, _return the minimum number of Fibonacci numbers whose
# sum is equal to_`k`. The same Fibonacci number can be used multiple times.
#
# The Fibonacci numbers are defined as:
#
#   * `F1 = 1`
#   * `F2 = 1`
#   * `Fn = Fn-1 + Fn-2` for `n > 2.`
#
# It is guaranteed that for the given constraints we can always find such
# Fibonacci numbers that sum up to `k`.
#
#
#
# **Example 1:**
#
#
#
#     Input: k = 7
#     Output: 2
#     Explanation: The Fibonacci numbers are: 1, 1, 2, 3, 5, 8, 13, ...
#     For k = 7 we can use 2 + 5 = 7.
#
# **Example 2:**
#
#
#
#     Input: k = 10
#     Output: 2
#     Explanation: For k = 10 we can use 2 + 8 = 10.
#
#
# **Example 3:**
#
#
#
#     Input: k = 19
#     Output: 3
#     Explanation: For k = 19 we can use 1 + 5 + 13 = 19.
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using MultiResolutionIterators
using Test
using Base.Iterators

@testset "basic 1" begin
    eg =   [["aaaa", "bbbb", "ccc"], ["AAA", "BB", "CCC", "DDDDD"], ["111","222"]]

    @testset "join" begin
        @test full_consolidate(join_levels(eg, Dict(2=>" "))) ==
            full_consolidate(["aaaa bbbb ccc","AAA BB CCC DDDDD","111 222"])

        @test full_consolidate(join_levels(eg, Dict(1=>"+", 2=>" "))) ==
            full_consolidate("aaaa bbbb ccc+AAA BB CCC DDDDD+111 222")
    end

    @testset "consolidate" begin
        @test full_consolidate(eg) == eg
    end

    @testset "flatten_levels (flatten)" begin
        @test full_consolidate(flatten_levels(eg, 1)) ==
            full_consolidate(["aaaa", "bbbb", "ccc","AAA", "BB", "CCC", "DDDDD", "111","222"]) ==
            ["aaaa", "bbbb", "ccc","AAA", "BB", "CCC", "DDDDD", "111","222"]

        @test full_consolidate(flatten_levels(eg, 2)) == consolidate.(["aaaabbbbccc","AAABBCCCDDDDD","111222"])

        @test full_consolidate(flatten_levels(eg, 1:2)) ==
            full_consolidate(flatten_levels(eg, 1:2)) ==
            full_consolidate(flatten_levels(eg, ALL_LEVELS))
            "aaaabbbbcccAAABBCCCDDDDD111222"

    end
end

@testset "turtles" begin
    # From
    #  https://simple.wikipedia.org/wiki/Turtle ,
    #  https://simple.wikipedia.org/wiki/Sea_turtle ,
    #  https://simple.wikipedia.org/wiki/Green_turtle

    multi_tokenized = [
        [[["Turtles", "are", "the", "reptile", "order", "Testudines", "."], ["They", "have", "a", "special", "bony", "or", "cartilaginous", "shell", "developed", "from", "their", "ribs", "that", "acts", "as", "a", "shield", "."]], [["The", "order", "Testudines","includes", "both", "living", "and", "extinct", "species", "."], ["The", "earliest", "fossil", "turtles", "date","from", "about", "220", "million", "years", "ago", "."], ["So", "turtles", "are", "one", "of", "the", "oldest", "surviving", "reptile", "groups", "and", "a", "more", "ancient", "group", "than", "lizards", ",", "snakes", "and", "crocodiles", "."]], [["Turtle", "have", "been", "very", "successful", ",", "and", "have", "almost", "world-wide", "distribution", "."], ["But", ",", "of", "the", "many", "species", "alive", "today", ",", "some", "are","highly", "endangered", "."]]],
        [[["Sea", "turtles", "(", "Chelonioidea", ")", "are", "turtles", "found", "in", "all", "the", "world", "'s", "oceans", "except", "the", "Arctic", "Ocean", ",", "and", "some", "species", "travel", "between", "oceans", "."], ["The", "term", "is", "US", "English", "."],["In", "British", "English", "they", "are", "simply", "called", "``", "turtles", "''", ";", "fresh-water", "chelonians", "are", "called", "``", "terrapins", "''", "and", "land", "chelonians", "are", "called", "tortoises", "."]], [["There", "are", "seven", "types", "of", "sea", "turtles", ":", "Kemp", "'s", "Ridley", ",", "Flatback", ",", "Green", ",", "Olive", "Ridley", ",", "Loggerhead", ",", "Hawksbill", "and", "the", "leatherback", "."], ["All", "but", "the", "leatherback", "are", "in", "the", "family", "Chelonioidea", "."], ["The", "leatherback", "belongs", "to", "the", "family", "Dermochelyidae", "and", "is", "its", "only", "member", "."], ["The", "leatherback", "sea", "turtle", "is", "the", "largest", ",", "measuring", "six", "or", "seven", "feet", "(", "2", "m", ")", "in", "length", "at", "maturity", ",", "and", "three", "to", "five", "feet", "(", "1", "to", "1.5", "m", ")", "in", "width", ",", "weighing", "up", "to", "2000", "pounds", "(", "about", "900", "kg", ")", "."], ["Most", "other", "species", "are", "smaller", ",", "being", "two", "to", "four", "feet", "in", "length", "(", "0.5", "to", "1", "m", ")", "and", "proportionally", "less", "wide", "."], ["The", "Flatback", "turtle", "is", "found", "solely", "on", "the", "northern", "coast", "of", "Australia", "."]]],
        [[["Chelonia", "mydas", ",", "commonly", "known", "as", "the", "green", "turtle", ",", "is", "a", "large", "sea", "turtle", "belonging", "to", "the", "family", "Cheloniidae", "."], ["It", "is", "the", "only", "species", "in", "its", "genus", "."], ["It", "is", "one", "of", "the", "seven", "marine", "turtles", ",", "which", "are", "all", "endangered", "."]], [["Although", "it", "might", "have", "some", "green", "on", "its", "carapace", "(", "shell", ")", ",", "the", "green", "turtle", "is", "not", "green", "."], ["It", "gets", "its", "name", "from", "the", "fact", "that", "its", "body", "fat", "is", "green","."], ["It", "can", "grow", "up", "to", "1", "m", "(", "3", "ft", ")", "long", "and", "weigh", "up", "to", "160","kg", "(", "353", "lb", ")", "."], ["They", "are", "an", "endangered", "species", ",", "especially", "in", "Florida", "and", "the", "Pacific", "coast", "of", "Mexico", "."], ["They", "can", "also", "be", "found", "in", "warm", "waters", "around", "the", "world", "and", "are", "found", "along", "the", "coast", "of", "140", "countries", "."]], [["The", "female", "turtle", "lays", "eggs", "in", "nests", "she", "builds", "in", "the", "sand", "on", "the", "beaches", "."], ["She", "uses", "the", "same", "beach", "that", "she", "was", "born", "on", "."], ["During", "the", "nesting", "season", "in", "summer", "she", "can", "make", "up", "to", "five", "nests", "."], ["She", "can", "lay", "as", "many", "as", "135", "eggs", "in", "a", "nest", "."], ["The", "eggs", "take","about", "two", "months", "to", "hatch", "."], ["The", "baby", "turtles", "are", "about", "50", "mm", "(", "2", "in", ")", "in", "length", "."]]]
        ]

    documents = full_consolidate(flatten_levels(multi_tokenized, 0))
    #corpus, doc, para, sent, word, char
    @test length(documents)==3
    @test typeof(documents[1][1][1][1])==String # words
    @test typeof(documents[1][1][1][1][1])==Char

    paras = full_consolidate(flatten_levels(multi_tokenized, 1))
    #corpus, para, sent, word, char
    @test length(paras) == sum(length.(documents))
    @test typeof(paras[1][1][1])==String # words
    @test typeof(paras[1][1][1][1])==Char

    words = full_consolidate(flatten_levels(multi_tokenized, 1:3))
    #corpus, word, char
    @test typeof(words[1])==String
    @test typeof(words[end])==String

    docs_of_words = full_consolidate(flatten_levels(multi_tokenized, 2:3))
    #corpus, doc, word, char
    @test length(docs_of_words)==3
    @test typeof(docs_of_words[1][1])==String

    @test full_consolidate(flatten_levels(docs_of_words, 1)) ==
          full_consolidate(flatten_levels(paras, 1:2)) ==
          words

    chars = full_consolidate(flatten_levels(multi_tokenized, ALL_LEVELS))
    @test typeof(chars[1]) == Char

    @test full_consolidate(flatten_levels(docs_of_words, ALL_LEVELS)) ==
          full_consolidate(flatten_levels(paras, ALL_LEVELS)) ==
          full_consolidate(flatten_levels(words, ALL_LEVELS)) ==
          chars
end


text:
using Test, OrdinaryDiffEq, HyperbolicDiffEq

u₀func(x) = sinpi(x)

balance_law = Burgers()
@inferred UniformPeriodicMesh1D(0., 2., 50)
meshx = UniformPeriodicMesh1D(0., 2., 50)
@inferred LobattoLegendre(4)
basis = LobattoLegendre(4)
tspan = (0.0, 0.3)

@inferred UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, EnergyConservativeFlux(), GodunovFlux())
semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, EnergyConservativeFlux(), GodunovFlux())
semidiscretise(semidisc, u₀func, tspan)
ode = semidiscretise(semidisc, u₀func, tspan)
sol = solve(ode, SSPRK104(), dt=max_dt(tspan[1], ode.u0, semidisc), save_everystep=false)
@inferred evaluate_coefficients(sol.u[end], semidisc)
x, u = evaluate_coefficients(sol.u[end], semidisc)

@inferred UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, L2L4ConservativeFlux(), GodunovFlux())
semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, L2L4ConservativeFlux(), GodunovFlux())
ode = semidiscretise(semidisc, u₀func, tspan)
sol = solve(ode, SSPRK104(), dt=max_dt(tspan[1], ode.u0, semidisc), save_everystep=false)
@inferred evaluate_coefficients(sol.u[end], semidisc)
x, u = evaluate_coefficients(sol.u[end], semidisc)

@inferred UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, L2L2sConservativeFlux(Val{3}()), GodunovFlux())
semidisc = UniformPeriodicFluxDiffDisc1D(balance_law, meshx, basis, L2L2sConservativeFlux(Val{3}()), GodunovFlux())
ode = semidiscretise(semidisc, u₀func, tspan)
sol = solve(ode, SSPRK104(), dt=max_dt(tspan[1], ode.u0, semidisc), save_everystep=false)
@inferred evaluate_coefficients(sol.u[end], semidisc)
x, u = evaluate_coefficients(sol.u[end], semidisc)

balance_law = HyperbolicDiffEq.Euler{Float64,2}(1.4)
ϱ₀(x,y)  = 1.
vx₀(x,y) = sin(x)*cos(y)
vy₀(x,y) = -cos(x)*sin(y)
p₀(x,y)  = 100/1.4 + (cos(2x)+cos(2y))*3/16
u₀(x,y) = conserved_variables(ϱ₀(x,y), vx₀(x,y), vy₀(x,y), p₀(x,y), balance_law)
xmin=0.; xmax=2π; ymin=0.; ymax=2π
fvol = KennedyGruberFlux()
fnum = fvol
p = 9
Nx = 5
Ny = 5
meshx = UniformPeriodicMesh1D(xmin, xmax, Nx)
meshy = UniformPeriodicMesh1D(ymin, ymax, Ny)
basis = LobattoLegendre(p)
tspan = (0., 14.)
semidisc = UniformPeriodicFluxDiffDisc2D(balance_law, meshx, meshy, basis, fvol, fnum, true)
ode = semidiscretise(semidisc, u₀, tspan)
tstops = range(tspan[1], tspan[end], length=100)
sol = solve(ode, SSPRK33(), dt=0.1/((2p+1)*max(Nx,Ny)), saveat=tstops, tstops=tstops)


text:
# Main entry point to the function

# MCEB +-

# mceb(DiscretizedNormalSamples, Target;
#        PriorFamily = ,
#        alpha=0.1,
#        nbhood=nothing)

# -

"""
    _split_data(n_total::Integer)
    _split_data(Zs::AbstractVector)
	_split_data(Zs::AbstractVector, ::Symbol)

Helper functions to partition all observations
into two disjoint folds.
"""
function _split_data(n_total::Integer)
    n_half = ceil(Int, n_total/2)
    idx_test = sample(1:n_total, n_half, replace=false)
    idx_train = setdiff(1:n_total, idx_test)
    idx_train, idx_test
end

function _split_data(Zs::AbstractVector)
    n = length(Zs)
    idx_train, idx_test = _split_data(n)
    Zs[idx_train], Zs[idx_test], idx_train, idx_test
end

function _split_data(Zs::AbstractVector, sym::Symbol)
    if sym==:random
    	_split_data(Zs)
	else
		throw(DomainError("Only :random implemented."))
	end
end

function _default_tuner(delta_tuner::DeltaTuner, Zs_test_discr, fkde_train)
	delta_tuner
end

function _default_tuner(delta_tuner::Type{<:DeltaTuner}, Zs_test_discr, fkde_train::KDEInfinityBand)
	n = nobs(Zs_test_discr)
	δ_min = default_δ_min(n, fkde_train.C∞)
	delta_tuner(n, δ_min)
end

function _default_tuner(delta_tuner::Type{<:DeltaTuner}, Zs_test_discr, ::Nothing)
	n = nobs(Zs_test_discr)
	δ_min = 0.0
	delta_tuner(n, δ_min)
end

function _default_tuner(Zs_test_discr, fkde_train)
	_default_tuner(RMSE, Zs_test_discr, fkde_train)
end

# train => first fold
# test -> second fold

"""
	MinimaxCalibratorSetup(;Zs_test,
                           Zs_train,
                           fkde_train,
                           prior_class,
                           Zs_test_discr,
                           delta_tuner,
                           pilot_method,
                           cache_target)

## Description:
This creates the object that holds all preliminary calculations necessary to conduct
inference for empirical Bayes estimands. It can be either constructed manually, or more typically
through `fit(opt::MinimaxCalibratorOptions, Zs)` where `Zs` comprises all the data and `opt`
is a [`MinimaxCalibratorOptions`](@ref) object. See `data_analysis.html` for examples of
both approaches.

## Arguments:
* `Zs_test`: Samples in the second data fold.
* `Zs_train`: Samples in the first data fold.
* `fkde_train`: A [`KDEInfinityBand`](@ref) instance computed based on `Zs_train` or external data (but not `Zs_test`).
* `Zs_test_discr`: A discretization of `Zs_test` into [`DiscretizedStandardNormalSamples`](@ref).
* `delta_tuner`: Specifies how the bias-variance tradeoff (i.e., ``\\delta``) will be navigated, for example [`RMSE`](@ref) or [`HalfCIWidth`](@ref).
* `pilot_method`: Method for computing the pilot estimator for the numerator and denominator of the empirical Bayes estimand. Should support the interface `estimate(target::LinearEBayesTarget, pilot_method, Zs_train)`.
* `cache_target`: Defaults to `false`. Setting it to `true` will speed up computations when multiple similar targets are being estimated (e.g., PosteriorMeans along a grid).
"""
mutable struct MinimaxCalibratorSetup{DS <: DiscretizedStandardNormalSamples,
                                           IDX,
										   ESTR,
                                           ESTE,
                                           Gcal <: ConvexPriorClass,
                                           DT <: DeltaTuner,
                                           P,
                                           NB}
    Zs_train::ESTR
    Zs_test::ESTE
    idx_train::IDX
    idx_test::IDX
    prior_class::Gcal
    fkde_train::NB
    Zs_test_discr::DS
    delta_tuner::DT
    pilot_method::P
	cache_target::Bool
	previous_target
	delta_cached::Float64
end

function MinimaxCalibratorSetup(;Zs_test,
	                             prior_class,
								 Zs_test_discr,
								 delta_tuner = _default_tuner(Zs_test_discr, fkde_train),
								 Zs_train = nothing,
								 idx_train = nothing,
								 idx_test = nothing,
								 fkde_train = nothing,
								 pilot_method = nothing,
								 cache_target = false,
								 previous_target = nothing,
								 delta_cached = NaN) #default_pilot?

	delta_tuner = _default_tuner(delta_tuner, Zs_test_discr, fkde_train)

	MinimaxCalibratorSetup(Zs_train, Zs_test, idx_train, idx_test,
	                       prior_class, fkde_train, Zs_test_discr,
						   delta_tuner, pilot_method,
						   cache_target, previous_target, delta_cached)

end

#update_delta_tuner()

broadcastable(setup::MinimaxCalibratorSetup) = Ref(setup)

function StatsBase.fit(mceb_setup::MinimaxCalibratorSetup, target::LinearEBayesTarget)
	@unpack Zs_test_discr, prior_class, delta_tuner = mceb_setup
	model = initialize_modulus_problem(Zs_test_discr, prior_class, target)
	opt_δ = optimal_δ(model, delta_tuner)
	modulus_at_δ!(model, opt_δ)
	sm  = SteinMinimaxEstimator(Zs_test_discr,
	                            prior_class,
	                            target,
	                            model)
end

"""
	CalibratedMinimaxEstimator(target::PosteriorTarget, sme::SteinMinimaxEstimator,
	                           pilot, pilot_num pilot_denom)

Object that stores a `pilot` estimator (including both the numerator `pilot_num`
and the denominator `pilot_denom`, where `pilot = pilot_num/pilot_denom`), as well
the quasi-minimax estimator for the calibrated target.
"""
Base.@kwdef struct CalibratedMinimaxEstimator{T<:PosteriorTarget,
	                                          SME<:SteinMinimaxEstimator,
											  S<:Real}
	target::T
	sme::SME
	pilot::S
	pilot_num::S
	pilot_denom::S
end

function _update_setup(mceb_setup::MinimaxCalibratorSetup, target::EBayesTarget)
	@unpack cache_target, previous_target, delta_cached = mceb_setup

	if (cache_target) && (typeof(previous_target) == typeof(target))
		updated_setup = @set mceb_setup.delta_tuner = FixedDelta(delta_cached)
	else
		updated_setup = mceb_setup
	end
	updated_setup
end

function StatsBase.fit(mceb_setup::MinimaxCalibratorSetup, target::PosteriorTarget)
	@unpack Zs_train, Zs_test_discr,
	        prior_class, delta_tuner, pilot_method,
			cache_target = mceb_setup

	numerator_target = target.num_target
	denominator_target = MarginalDensityTarget(location(target.num_target)) #	perhaps expose this as MarginalDensityTarget(target)

	num_hat = estimate(numerator_target, pilot_method, Zs_train)
	denom_hat = estimate(denominator_target, pilot_method, Zs_train)
	theta_hat = num_hat/denom_hat

	calib_target = CalibratedTarget(posterior_target = target,
	                                θ̄ = theta_hat,
									denom = denom_hat)

	tmp_setup = _update_setup(mceb_setup, target)
    calib_fit = fit(tmp_setup, calib_target)

	if cache_target
		mceb_setup.delta_cached = calib_fit.δ
		mceb_setup.previous_target = target
	end

	CalibratedMinimaxEstimator(target=target,
                               sme=calib_fit,
							   pilot=theta_hat,
							   pilot_num=num_hat,
							   pilot_denom=denom_hat)

end

function target_bias_std(target::PosteriorTarget,
	                     calib::CalibratedMinimaxEstimator,
						 Zs=calib.sme.Z; kwargs...)

	calib_target = calib.target
	cres = target_bias_std(calib_target, calib.sme, Zs; kwargs...)

	ctarget = cres.estimated_target
	cbias = cres.estimated_bias
	cstd = cres.estimated_std

	pilot_denom = calib.pilot_denom
	# TODO: sanity check this is the same as target of denominator.
	pilot = calib.pilot

	(estimated_target = pilot + ctarget/pilot_denom,
	estimated_bias = cbias/pilot_denom,
	estimated_std = cstd/pilot_denom)
end

"""
	MinimaxCalibratorOptions(split = :random
	                         prior_class::ConvexPriorClass,
	                         marginal_grid,
	                         infinity_band_options,
	                         pilot_options,
	                         tuner = HalfCIWidth
	                         cache_target = false)

Struct that wraps all options required to automatically
create a [`MinimaxCalibratorSetup`](@ref) object. The options are as follows

* `split`: Defaults to `:random`, i.e., a random 50-50 split of the data.
* `prior_class`: A [`ConvexPriorClass`](@ref) in which the true ``G`` is believed to lie.
* `marginal_grid`: Marginal domain discretization grid.
* `infinity_band_options`: Options for constructing the localization, cf. [`KDKDEInfinityBandOptions`](@ref).
* `pilot_options`: 	Options that describe how to derive the pilot estimators, for example [`ButuceaComteOptions`](@ref).
* `tuner`: Specifies how the bias-variance tradeoff (i.e., ``\\delta``) will be navigated, for example [`RMSE`](@ref) or [`HalfCIWidth`](@ref).
* `cache_target`: Defaults to `false`. Setting it to `true` will speed up computations when multiple similar targets are being estimated (e.g., PosteriorMeans along a grid).
"""
Base.@kwdef struct MinimaxCalibratorOptions{SPL,
                                     GCAL <: ConvexPriorClass,
									 GR <: AbstractVector,
									 PS,
									 OPT,
									 D}
   split::SPL = :random
   prior_class::GCAL
   marginal_grid::GR
   infinity_band_options::OPT = KDEInfinityBandOptions(a_min = minimum(marginal_grid),
                                              a_max = maximum(marginal_grid))
   pilot_options::PS
   tuner::D = HalfCIWidth
   cache_target::Bool = false
end

function fit(mceb_opt::MinimaxCalibratorOptions, Zs::AbstractVector{<:StandardNormalSample})
	@unpack split, prior_class, marginal_grid,
	        infinity_band_options, pilot_options, tuner, cache_target = mceb_opt
	Zs_train, Zs_test, idx_train, idx_test = _split_data(Zs, split)

	a_min, a_max = extrema(marginal_grid)
	infinity_band_options = KDEInfinityBandOptions(a_min=a_min, a_max=a_max)
	# TODO: let fkde_train be better at getting dispatched
	fkde_train = fit(infinity_band_options, response.(Zs_train))

	pilot_fit = fit(pilot_options, Zs_train)

	Zs_test_discr = DiscretizedStandardNormalSamples(Zs_test, marginal_grid)
	Zs_test_discr = set_neighborhood(Zs_test_discr, fkde_train)

	mceb_setup = MinimaxCalibratorSetup(Zs_train = Zs_train,
	                                    Zs_test = Zs_test,
										idx_train = idx_train,
										idx_test = idx_test,
							            prior_class = prior_class,
							            fkde_train = fkde_train,
							            Zs_test_discr = Zs_test_discr,
							            pilot_method = pilot_fit,
										delta_tuner = tuner,
										cache_target= cache_target)
end


text:
Tables.istable(::Type{<:Result}) = true

# Rows

Tables.rowaccess(::Type{<:Result}) = true
Tables.rows(jl_result::Result) = jl_result

Base.eltype(jl_result::Result) = Row
Base.length(jl_result::Result) = num_rows(jl_result)

function Base.iterate(jl_result::Result, (len, row)=(length(jl_result), 1))
    row > len && return nothing
    return Row(jl_result, row), (len, row + 1)
end

function Tables.schema(jl_result::Result)
    types = map(jl_result.not_null, column_types(jl_result)) do not_null, col_type
        not_null ? col_type : Union{col_type, Missing}
    end
    return Tables.Schema(map(Symbol, column_names(jl_result)), types)
end

struct Row
    result::Result
    row::Int
end

result(pqrow::Row) = getfield(pqrow, :result)
row_number(pqrow::Row) = getfield(pqrow, :row)

Base.propertynames(pqrow::Row) = map(Symbol, column_names(result(pqrow)))

function Base.getproperty(pqrow::Row, name::Symbol)
    jl_result = result(pqrow)
    row = row_number(pqrow)
    col = column_number(jl_result, name)
    return jl_result[row, col]
end

function Base.getindex(pqrow::Row, col::Integer)
    row = row_number(pqrow)
    return result(pqrow)[row, col]
end

Base.length(pqrow::Row) = num_columns(result(pqrow))

function Base.iterate(pqrow::Row, (len, col)=(length(pqrow), 1))
    col > len && return nothing
    return (result(pqrow)[row_number(pqrow), col], (len, col + 1))
end

# Columns

struct Column{T} <: AbstractVector{T}
    result::Result
    col::Int
    col_name::Symbol
    oid::Oid
    not_null::Bool
    typ::Type
    func::Base.Callable
end

struct Columns <: AbstractVector{Column}
    result::Result
end

result(cs::Columns) = getfield(cs, :result)

Base.propertynames(cs::Columns) = map(Symbol, column_names(result(cs)))

Base.getproperty(cs::Columns, name::Symbol) = Column(result(cs), name)
Base.getindex(cs::Columns, col::Integer) = Column(result(cs), col)
Base.IndexStyle(::Type{Columns}) = IndexLinear()
Base.length(cs::Columns) = num_columns(result(cs))
Base.size(cs::Columns) = (length(cs),)

Tables.columnaccess(::Type{<:Result}) = true
Tables.columns(jl_result::Result) = Columns(jl_result)

function Tables.schema(cs::Columns)
    jl_result = result(cs)
    types = map(jl_result.not_null, column_types(jl_result)) do not_null, col_type
        not_null ? col_type : Union{col_type, Missing}
    end
    return Tables.Schema(map(Symbol, column_names(jl_result)), types)
end

function Column(jl_result::Result, col::Integer, name=Symbol(column_name(jl_result, col)))
    @boundscheck if !checkindex(Bool, Base.OneTo(num_columns(jl_result)), col)
        throw(BoundsError(Columns(jl_result), col))
    end

    oid = column_oids(jl_result)[col]
    typ = column_types(jl_result)[col]
    func = jl_result.column_funcs[col]
    not_null = jl_result.not_null[col]
    element_type = not_null ? typ : Union{typ, Missing}
    return Column{element_type}(jl_result, col, name, oid, not_null, typ, func)
end

function Column(jl_result::Result, name::Symbol, col=column_number(jl_result, name))
    return Column(jl_result, col, name)
end

result(c::Column) = getfield(c, :result)
column_number(c::Column) = getfield(c, :col)
column_name(c::Column) = getfield(c, :col_name)

function Base.getindex(c::Column{T}, row::Integer)::T where T
    jl_result = result(c)
    col = column_number(c)
    if isnull(jl_result, row, col)
        return missing
    else
        return c.func(PQValue{c.oid}(jl_result, row, col))::c.typ
    end
end

Base.IndexStyle(::Type{<:Column}) = IndexLinear()
Base.length(c::Column) = num_rows(result(c))
Base.size(c::Column) = (length(c),)

"""
    LibPQ.load!(table, connection::LibPQ.Connection, query) -> LibPQ.Statement

Insert the data from `table` using `query`.
`query` will be prepared as a [`LibPQ.Statement`](@ref) and then [`execute`](@ref) is run
on every row of `table`.

For best performance, wrap the call to this function in a PostgreSQL transaction:

```jldoctest; setup = :(execute(conn, "CREATE TEMPORARY TABLE libpqjl_test (no_nulls varchar(10) PRIMARY KEY, yes_nulls varchar(10));"))
julia> execute(conn, "BEGIN;");

julia> LibPQ.load!(
           (no_nulls = ["foo", "baz"], yes_nulls = ["bar", missing]),
           conn,
           "INSERT INTO libpqjl_test (no_nulls, yes_nulls) VALUES (\\\$1, \\\$2);",
       );

julia> execute(conn, "COMMIT;");
```
"""
function load!(table::T, connection::Connection, query::AbstractString) where T
    rows = Tables.rows(table)
    stmt = prepare(connection, query)
    state = iterate(rows)
    state === nothing && return stmt
    row, st = state
    names = propertynames(row)
    sch = Tables.Schema(names, nothing)
    parameters = Vector{Parameter}(undef, length(names))
    while state !== nothing
        row, st = state
        Tables.eachcolumn(sch, row) do val, col, nm
            parameters[col] = if ismissing(val)
                missing
            elseif val isa AbstractString
                convert(String, val)
            else
                string_parameter(val)
            end
        end
        close(execute(stmt, parameters; throw_error=true))
        state = iterate(rows, st)
    end
    return stmt
end


text:
module SelfPropelledModel
using Random

struct Particle{S}
    x::S
    v::S
    ox::S
    ov::S
end

function Particle(x::A,v::A,ox::A,ov::A) where {N,T,A<:NTuple{N,T}}
    Particle{A}(x,v,ox,ov)
end

function Particle(a_x,a_v)
    x=Tuple(a_x)
    v=Tuple(a_v)
    Particle(x,v,x,v)
end

function sqnorm(a)
    foldl((x,y)->x+y*y,a,init=zero(first(a)))
end

function θ(a,b)
    n1=a*a
    n2=b*b
    ns=sqrt(n1+n2)
    return a/ns,b/ns
end

function dist(v)
    N=length(v)
    d=zero(v[1].x[1])
    m=zero(v[1].x[1])
    id=0
    for i in 2:(N-1)
        nd=sqnorm(v[i].x.-v[i+1].x)+sqnorm(v[i].x.-v[i-1].x)
        if nd>=d
            id=i
            d=nd
        end
        m+=nd
    end
    d,id,m/N
end

function optimize!(v)
    for i in 1:(length(v)*length(v))
        D,j,mD=dist(v)
        r=rand(eachindex(v))
        v[j],v[r]=v[r],v[j]
        nD,nj,nmD=dist(v)
        if nmD>mD && nD>D
            v[j],v[r]=v[r],v[j]
        end
        r1=rand(eachindex(v))
        r2=rand(eachindex(v))
        v[j+1],v[r1]=v[r1],v[j+1]
        v[j-1],v[r2]=v[r2],v[j-1]
        nD,nj,nmD=dist(v)
        if nmD>mD && nD>D
            v[j-1],v[r2]=v[r2],v[j-1]
            v[j+1],v[r1]=v[r1],v[j+1]
        end
    end
end

function update_velocity(i::Int,P,α,γ,ρ)
    p=P[i]
    N=length(P)
    v0=sqrt(sqnorm(p.ov))
    p_n=P[ifelse(i<N,i+1,1)]
    p_b=P[ifelse(i>1,i-1,N)]
    η=randn(length(p.x))
    nvx=v0 .* θ(@.(α * p.ov + γ * (p_n.ov+p_b.ov)/2 + ρ*η)...)
    P[i]=Particle(p.x,nvx,p.ox,p.ov)
end

function update_position(i::Int,P,Δt)
    p=P[i]
    nx=p.x .+ p.v .* Δt
    P[i]=Particle(nx,p.v,p.x,p.v)
end

function spm_step!(P,α,γ,ρ,Δt)
    Threads.@threads for i in eachindex(P)
        update_velocity(i,P,α,γ,ρ)
    end
    Threads.@threads for i in eachindex(P)
        update_position(i,P,Δt)
    end
end

function sample_particles(N,x_sampler,v_sampler)
    sv=[Particle(x_sampler(),v_sampler()) for i in 1:N]
    optimize!(sv)
    optimize!(sv)
    sv
end

export sample_particles, spm_step!
end # module


text:
# Zero Momentum
function lam4piflow(k, m2, T, Npi, lamda4pi)
    return lamda4pi^2 * (8 + Npi) * (dkF1TildeAll(k, m2, T)
    + dkF2TildeAll(k, m2, T)
    )
end
function dkVintqs_Zero(q0, k, kprim, m, T, lam4pik, Npi; kwargs...)
    return lam4pik^2 *(2 + Npi) *(
        k^3 *
        (2 + Npi) *
        (dkF1TildeAll(kprim, m, T)
        + dkF2TildeAll(kprim, m, T)
        ) / 3 +
        6 * dkF1TildeintqsAll(q0, k, kprim, m, T; kwargs...)
        +
        6 * dkF2TildeintqsAll(q0, k, kprim, m, T; kwargs...)
    )
end

function Vintqs_Zero(k, T, Npi, hfun, UVScale; kwargs...)
    -2 * hquadrature(
        kprim -> dkVintqs_Zero(
            Epi(k, -hfun(k)[2]),
            k,
            kprim,
            -hfun(kprim)[2],
            T,
            hfun(kprim)[1],
            Npi;
            kwargs...,
        ),
        k,
        UVScale;
        kwargs...,
        # initdiv=200,
    )[1] + (2 * k^3 * hfun(UVScale)[1] * (2 + Npi)) / 3
end

function Vintqs_LPA(k, Npi, lamdapik)
    (2 * k^3 * lamdapik * (2 + Npi)) / 3
end

function propReZeroflow(k, T, Npi, hfun, UVScale; kwargs...)
    2 *
    Coeffgamm2Simple(k, T, -hfun(k)[2]) *
    Vintqs_Zero(k, T, Npi, hfun, UVScale; kwargs...)
end

function propReLPAflow(k, T, Npi, hfun)
    2 * Coeffgamm2Simple(k, T, -hfun(k)[2]) * Vintqs_LPA(k, Npi, hfun(k)[1])
end

function dkV4piResimple(p0, ps, q0, qs, costhe, k, m, T, lam4pik, Npi)
    return lam4pik^2 *
           (2 + Npi) *
           (
               2 * (2 + Npi) * (dkF1TildeAll(k, m, T) + dkF2TildeAll(k, m, T)) +
               6 * dkF1TildeAll(
                   p0 - q0,
                   max(0.1, sqrt(ps^2 + 2 * costhe * ps * qs + qs^2)),
                   k,
                   m,
                   T,
               ) +
               6 * dkF1TildeAll(
                   p0 + q0,
                   max(0.1, sqrt(ps^2 + 2 * costhe * ps * qs + qs^2)),
                   k,
                   m,
                   T,
               ) +
               6 * (
                   dkF2TildeAll(
                       p0 - q0,
                       max(0.1, sqrt(ps^2 + 2 * costhe * ps * qs + qs^2)),
                       k,
                       m,
                       T,
                   ) + dkF2TildeAll(
                       p0 + q0,
                       max(0.1, sqrt(ps^2 + 2 * costhe * ps * qs + qs^2)),
                       k,
                       m,
                       T,
                   )
               )
           )
end

@doc raw"""
    dkVReintqs(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)

compute $\int_0^{qsmax}dq_s qs^2\int_{-1}^{1}d\cos\theta \tilde{\partial_k}\mathrm{Re}V(q_0)$.

`dkVImintqs` only contains $V(q_0)$, for $-q_0$, we have $\int d\cos\theta V(q_0)=\int d\cos\theta V(-q_0)$,
so we need an extra $2$ at somewhere.

`dkVReintqs` contains type-1 and type-2 delta function

# Arguments
- `qsmax`: we integrate $q_s$ from $0$ to $k$, `qsmax` will set to `k` when we do the integration $dk'$, it should be distinguished from $k'$
- `m`: mass square, it will be $m(k')$ when we do the integration $dk'$.
- `lam4pik`: $\lambda_{4\pi}$, it will be $\lambda_{4\pi}(k')$ when we do the integration $dk'$ .
"""
function dkVReintqs(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik; kwarg...)
    return lam4pik^2 * (2 + Npi) * (3 * (
            dkF1TildeintqsAll(p0 - q0, ps, qsmax, k, m, T; kwarg...) +
            dkF1TildeintqsAll(p0 + q0, ps, qsmax, k, m, T; kwarg...)
            +
            dkF2TildeintqsAll(p0 - q0, ps, qsmax, k, m, T; kwarg...) +
            dkF2TildeintqsAll(p0 + q0, ps, qsmax, k, m, T; kwarg...)
        ) +
        (Npi + 2) * 2 / 3 *
        qsmax^3 *
        (dkF1TildeAll(k, m, T)
        + dkF2TildeAll(k, m, T)
        )
    )
end

function dkVReintqs_Compensate(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    return lam4pik^2 *
           (2 + Npi) *
           (
               3 * (
                   dkF1TildeintqsAll_Compensate(p0 - q0, ps, qsmax, k, m, T) +
                   dkF1TildeintqsAll_Compensate(p0 + q0, ps, qsmax, k, m, T)
               )
           )
end

function dkVReintqs_delta1(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    return lam4pik^2 *
           (2 + Npi) *
           (
               3 * (
                   dkF1TildeintqsAll_delta1(p0 - q0, ps, qsmax, k, m, T) +
                   dkF1TildeintqsAll_delta1(p0 + q0, ps, qsmax, k, m, T)
               )
           )
end

function dkVReintqs_delta2(p0, ps, q0, qsmax, k, m, T, Npi, lam4pik)
    return lam4pik^2 *
           (2 + Npi) *
           (
               3 * (
                   dkF1TildeintqsAll_delta2(p0 - q0, ps, qsmax, k, m, T) +
                   dkF1TildeintqsAll_delta2(p0 + q0, ps, qsmax, k, m, T)
               )
           )
end

@doc raw"""
    VReintqs(p0, ps, k, T, Npi,IRScale,UVScale, mfun, lamfun)

compute $\int_0^{k}dq_s qs^2\int_{-1}^{1}d\cos\theta \mathrm{Re}V(q_0,k)$.
In our code, we perform integration over `kprim`, `q0` & `qs` does not involved,
so `qs=k`, `q0=Epi(k, mfun(k))`.

`VReintqs` contains type-1 and type-2 delta function.

# Arguments
- `mfun::Function`: $m^2(k)$, input from zero momentum result
- `lampifun::Function`: $\lambda_{4\pi}(k)$, input from zero momentum result.
"""
function VReintqs(p0, ps, k, T, Npi, IRScale, UVScale, mfun, lamfun; kwarg...)
    -hquadrature(
        kprim -> dkVReintqs(
            p0,
            ps,
            Epi(k, mfun(k)),
            k,
            kprim,
            mfun(kprim),
            T,
            Npi,
            lamfun(kprim);
            kwarg...,
        ),
        k,
        UVScale;
        kwarg...,
    )[1] + 2 * (k^3 * lamfun(UVScale) * (2 + Npi)) / 3
end

function propReintqs(p0, ps, T, IRScale, UVScale, Npi, mfun, lamfun; kwarg...)
    -hquadrature(
        k ->
            2 *
            VReintqs(
                p0,
                ps,
                k,
                T,
                Npi,
                IRScale,
                UVScale,
                mfun,
                lamfun;
                kwarg...,
            ) *
            Coeffgamm2(k, T, mfun),
        IRScale,
        UVScale,
        rtol=1e-8,
        atol=1e-8,
    )[1] + p0^2 - ps^2 - mfun(UVScale)
end

function fastpropReintqs(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVReintqs(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale);
                kwarg...,
            ) *
            Coeffgamm2(x[1], T, mfun) +
            2 *
            Coeffgamm2(x[1], T, mfun) *
            (2 * x[1]^3 * lamfun(UVScale) * (2 + Npi)) / 3,
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1] + p0^2 - ps^2 - mfun(UVScale)
end

function fastpropReintqs_Compensate(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVReintqs_Compensate(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropReintqs_delta1(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVReintqs_delta1(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropReintqs_delta2(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVReintqs_delta2(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale),
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1]
end

function fastpropReintqs_All(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return fastpropReintqs(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           ) +
           fastpropReintqs_Compensate(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           ) +
           fastpropReintqs_delta1(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           ) +
           fastpropReintqs_delta2(
               p0,
               ps,
               T,
               IRScale,
               UVScale,
               Npi,
               mfun,
               lamfun;
               kwarg...,
           )
end

function fastpropReintqs2(
    p0,
    ps,
    T,
    IRScale,
    UVScale,
    Npi,
    mfun,
    lamfun;
    kwarg...,
)
    return -hcubature(
        x ->
            (x[1] - UVScale) *
            2 *
            dkVReintqs(
                p0,
                ps,
                Epi(x[1], mfun(x[1])),
                x[1],
                x[2] * (x[1] - UVScale) + UVScale,
                mfun(x[2] * (x[1] - UVScale) + UVScale),
                T,
                Npi,
                lamfun(x[2] * (x[1] - UVScale) + UVScale);
                kwarg...,
            ) *
            Coeffgamm2(x[1], T, mfun),
        [IRScale, 0.0],
        [UVScale, 1.0];
        kwarg...,
    )[1] - hquadrature(
        k ->
            2 *
            Coeffgamm2(k, T, mfun) *
            (2 * (k^3 * lamfun(UVScale) * (2 + Npi)) / 3),
        IRScale,
        UVScale;
        kwarg...,
    )[1] + p0^2 - ps^2 - mfun(UVScale)
end


text:
# This file contains the functions and types related to the

abstract type  AdaptationAlgorithm end

type noAdaptation <: AdaptationAlgorithm
  Cov::Array{Float64}
end

type AMUpdate <: AdaptationAlgorithm
  C_0::Array{Float64}
  r_cov_m::Float64
  gamma_0::Float64
  k::Float64
  t_0::Int64
end

type AMUpdate_gen <: AdaptationAlgorithm
  C_0::Array{Float64}
  r_cov_m::Float64
  P_star::Float64
  gamma_0::Float64
  k::Float64
  t_0::Int64
end

# Function for different updating algorithms

# set up functions

"Set parameters for noAdaptation"
function set_adaptive_alg_params(algorithm::noAdaptation, nbr_of_unknown_parameters::Int64, theta::Vector,R::Int64)

  return (algorithm.Cov, NaN)

end

"Set parameters for AMUpdate"
function set_adaptive_alg_params(algorithm::AMUpdate, nbr_of_unknown_parameters::Int64, theta::Vector,R::Int64)

  # define m_g m_g_1
  m_g = zeros(nbr_of_unknown_parameters,1)
  m_g_1 = zeros(nbr_of_unknown_parameters,1)

  return (algorithm.C_0, algorithm.gamma_0, algorithm.k, algorithm.t_0, [algorithm.r_cov_m], m_g, m_g_1)

end

"Set parameters for AMUpdate_gen"
function set_adaptive_alg_params(algorithm::AMUpdate_gen, nbr_of_unknown_parameters::Int64, theta::Vector,R::Int64)

  Cov_m = algorithm.C_0
  log_r_cov_m = log(algorithm.r_cov_m)
  log_P_star = log(algorithm.P_star)
  gamma_0 = algorithm.gamma_0
  k = algorithm.k
  t_0 = algorithm.t_0
  vec_log_r_cov_m = zeros(1,R)
  vec_log_r_cov_m[1:t_0] = log_r_cov_m
  # define m_g m_g_1
  m_g = zeros(nbr_of_unknown_parameters,1)
  m_g_1 = zeros(nbr_of_unknown_parameters,1)
  #diff_a_log_P_start = 1.
  #factor = 1.
  #gain = 1.

  return (Cov_m, log_P_star, gamma_0, k, t_0, vec_log_r_cov_m, [log_r_cov_m], m_g, m_g_1)

end

# return covariance functions

doc"""
    return_covariance_matrix(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)
"""
function return_covariance_matrix(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)

  return adaptive_update_params[1]

end

doc"""
    return_covariance_matrix(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)

Print function for AMUpdate.
"""
function return_covariance_matrix(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return r_cov_m^2*Cov_m
end

doc"""
    return_covariance_matrix(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)

Print function for AMUpdate_gen.
"""
function return_covariance_matrix(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return exp(log_r_cov_m)^2*Cov_m
end

# print_covariance functions

doc"""
    print_covariance(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)
"""
function print_covariance(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)

  println(adaptive_update_params[1])

end

doc"""
    print_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)
"""
function print_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    println(r_cov_m^2*Cov_m)
end

doc"""
    print_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
"""
function print_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    println(exp(log_r_cov_m)^2*Cov_m)
end

# get_covariance functions

doc"""
    get_covariance(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)
"""
function get_covariance(algorithm::noAdaptation, adaptive_update_params::Tuple,r::Int64)

  return adaptive_update_params[1]

end

doc"""
    get_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)
"""
function get_covariance(algorithm::AMUpdate, adaptive_update_params::Tuple,r::Int64)
    r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return r_cov_m^2*Cov_m
end

doc"""
    get_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
"""
function get_covariance(algorithm::AMUpdate_gen, adaptive_update_params::Tuple,r::Int64)
    log_r_cov_m = adaptive_update_params[end-2][1]
    Cov_m = adaptive_update_params[1]
    return exp(log_r_cov_m)^2*Cov_m
end

# gaussian random walk functions

doc"""
    gaussian_random_walk(algorithm::noAdaptation, adaptive_update_params::Tuple, Theta::Vector, r::Int64)
"""
function gaussian_random_walk(algorithm::noAdaptation, adaptive_update_params::Tuple, Theta::Vector, r::Int64)

  return rand(MvNormal(Theta, adaptive_update_params[1])), zeros(size(Theta))

end

doc"""
    gaussian_random_walk(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Vector, r::Int64)
"""
function gaussian_random_walk(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Vector, r::Int64)
  r_cov_m = adaptive_update_params[end-2][1]
  Cov_m = adaptive_update_params[1]
  return rand(MvNormal(Theta, r_cov_m^2*Cov_m)), zeros(size(Theta))
end

doc"""
    gaussian_random_walk(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Vector, r::Int64)

RW for AMUpdate_gen.
"""
function gaussian_random_walk(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Vector, r::Int64)
  log_r_cov_m = adaptive_update_params[end-2][1]
  Cov_m = adaptive_update_params[1]
  return rand(MvNormal(Theta, exp(log_r_cov_m)^2*Cov_m)), zeros(size(Theta))
end

# functions for adaptation of parameters
doc"""
    adaptation(algorithm::noAdaptation, adaptive_update_params::Tuple, Theta::Array, r::Int64,a_log::Float64)
"""
function adaptation(algorithm::noAdaptation, adaptive_update_params::Tuple, Theta::Array, r::Int64,a_log::Float64)

  # do nothing

end

doc"""
    adaptation(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Array, r::Int64,a_log::Float64)
"""
function adaptation(algorithm::AMUpdate, adaptive_update_params::Tuple, Theta::Array, r::Int64,a_log::Float64)

  Cov_m = adaptive_update_params[1]
  m_g = adaptive_update_params[end-1]
  m_g_1 = adaptive_update_params[end]
  k = adaptive_update_params[3]
  gamma_0 = adaptive_update_params[2]
  t_0 = adaptive_update_params[4]

  g = r-1
  if r-1 == t_0
      m_g = mean(Theta[:,1:r-1],2)
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:, g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
  elseif r-1 > t_0
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:, g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
  end

end

doc"""
    adaptation(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Array, r::Int64,a_log::Float64)
"""
function adaptation(algorithm::AMUpdate_gen, adaptive_update_params::Tuple, Theta::Array, r::Int64, a_log::Float64)

  Cov_m = adaptive_update_params[1]
  m_g = adaptive_update_params[end-1]
  m_g_1 = adaptive_update_params[end]
  log_P_star = adaptive_update_params[2]
  log_r_cov_m = adaptive_update_params[end-2][1]
  k = adaptive_update_params[4]
  gamma_0 = adaptive_update_params[3]
  t_0 = adaptive_update_params[5]

  g = r-1;
  if r-1 >= t_0
    diff_a_log_P_start = abs(min(1, exp(a_log)) - exp(log_P_star)) #abs( min(log(1), a_log) - log_P_star)
    factor = 1 #diff_a_log_P_start/abs(log_r_cov_m)
    gain = diff_a_log_P_start/factor
    if min(1, exp(a_log)) < exp(log_P_star) #min(log(1), a_log) < log_P_star
      sign_val = -1
    else
      sign_val = 1
    end
    adaptive_update_params[end-2][1] = log_r_cov_m + sign_val * ( gamma_0/( (g+1)^k ) ) * gain
    adaptive_update_params[6][g+1] = log_r_cov_m
    if r-1 == t_0
      m_g = mean(Theta[:,1:r-1],2)
      m_g_1 = m_g + ( gamma_0/( (g+1)^k ) ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:,g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
    elseif r-1 > t_0
      m_g_1 = m_g + gamma_0/( (g+1)^k ) * (Theta[:,g+1] - m_g)
      adaptive_update_params[1][:] = Cov_m + ( gamma_0/( (g+1)^k ) ) * ( (Theta[:,g+1] - m_g)*(Theta[:,g+1] - m_g)'     -  Cov_m)
      adaptive_update_params[end-1][:] = m_g_1
    end
  end

end

max_stars_repo_name: SamuelWiqvist/abcmcmcrickermodel

max_issues_repo_name: SamuelWiqvist/abcmcmcrickermodel

max_forks_repo_name: SamuelWiqvist/abcmcmcrickermodel


text:
"""
    ARCHLMTest <: HypothesisTest
Engle's (1982) LM test for autoregressive conditional heteroskedasticity.
"""
struct ARCHLMTest{T<:Real} <: HypothesisTest
    n::Int         # number of observations
    p::Int         # number of lags
    LM::T    # test statistic
end

"""
    ARCHLMTest(am::UnivariateARCHModel, p=max(o, p, q, ...))
Conduct Engle's (1982) LM test for autoregressive conditional heteroskedasticity with
p lags in the test regression.
"""
ARCHLMTest(am::UnivariateARCHModel, p=presample(typeof(am.spec))) = ARCHLMTest(residuals(am), p)

"""
    ARCHLMTest(u::Vector, p::Integer)
Conduct Engle's (1982) LM test for autoregressive conditional heteroskedasticity with
p lags in the test regression.
"""
function ARCHLMTest(u::Vector{T}, p::Integer) where T<:Real
    @assert p>0
    n = length(u)
    u2 = u.^2
    X = zeros(T, (n-p, p+1))
    X[:, 1] .= one(eltype(u))
    for i in 1:p
        X[:, i+1] = u2[p-i+1:n-i]
    end
    y = u2[p+1:n]
    B = X \ y
    e = y - X*B
    ybar = y .- mean(y)
    LM = n * (1 - (e'e)/(ybar'ybar)) #T*R^2
    ARCHLMTest(n, p, LM)
end

testname(::ARCHLMTest) = "ARCH LM test for conditional heteroskedasticity"
population_param_of_interest(x::ARCHLMTest) = ("T⋅R² in auxiliary regression", 0, x.LM)
function show_params(io::IO, x::ARCHLMTest, ident)
    println(io, ident, "sample size:                    ", x.n)
    println(io, ident, "number of lags:                 ", x.p)
    println(io, ident, "LM statistic:                   ", x.LM)
end

pvalue(x::ARCHLMTest) = pvalue(Chisq(x.p), x.LM; tail=:right)

"""
    DQTest <: HypothesisTest
Engle and Manganelli's (2004) out-of-sample dynamic quantile test.
"""
struct DQTest{T<:Real} <: HypothesisTest
    n::Int         # number of observations
    p::Int         # number of lags
    level::T       # VaR level
    DQ::T          # test statistic
end

"""
    DQTest(data, vars, level, p=1)
Conduct Engle and Manganelli's (2004) out-of-sample dynamic quantile test with
p lags in the test regression. `vars` shoud be a vector of out-of-sample Value at Risk
predictions at level `level`.
"""
function DQTest(data::Vector{T}, vars::Vector{T}, level::AbstractFloat, p::Integer=1) where T<:Real
    @assert p>0
    @assert length(data) == length(vars)
    n = length(data)
    hit = (data .< -vars).*1 .- level
    y = hit[p+1:n]
    X = zeros(T, (n-p, p+2))
    X[:, 1] .= one(T)
    for i in 1:p
        X[:, i+1] = hit[p-i+1:n-i]
    end
    X[:, p+2] = vars[p+1:n]
    B = X \ y
    DQ = B' * (X'*X) *B/(level*(1-level)) # y'X * inv(X'X) * X'y / (level*(1-level)); note 2 typos in the paper
    DQTest(n, p, level, DQ)
end

testname(::DQTest) = "Engle and Manganelli's (2004) DQ test (out of sample)"
population_param_of_interest(x::DQTest) = ("Wald statistic in auxiliary regression", 0, x.DQ)
function show_params(io::IO, x::DQTest, ident)
    println(io, ident, "sample size:                    ", x.n)
    println(io, ident, "number of lags:                 ", x.p)
    println(io, ident, "VaR level:                      ", x.level)
    println(io, ident, "DQ statistic:                   ", x.DQ)
end

pvalue(x::DQTest) = pvalue(Chisq(x.p+2), x.DQ; tail=:right)


text:
# # Variable selection with the VIF

using SimpleSDMLayers
using GLM
using Statistics

# **Justification for this use case:** a lot of predictive variables are
# auto-correlated, and therefore one might argue that maybe, we may eventually
# build better models by removing some of them. This is generally refered to as
# variable selection, and sharing one's opinion on this is the fastest way to
# start a brawl at any gathering of ecologists.

# We will illustrate variable selection with the Variance Inflation Factor using
# the bioclim data from the region in which *Hypomyces lactifluorum* is found.

layers = SimpleSDMPredictor(
    WorldClim, BioClim, 1:19; left=-169.0, right=-50.0, bottom=24.0, top=71.0
);

# We will first gather everything in a matrix:

x = hcat([layer[keys(layer)] for layer in layers]...);
size(x)

# Because of the spread of some values, we will center and reduce this matrix to
# give every variable a mean of 0 and unit variance:

X = (x .- mean(x; dims=1)) ./ std(x; dims=1)
round.(Int, mean(X; dims=1))

# The VIF is simply measured as 1/(1-R²) by regressing every variable against
# all others. Let's have an illustration with the first predictor:

function vif(model)
    R² = r2(model)
    return 1 / (1 - R²)
end

vif(lm(X[:, 2:end], X[:, 1]))

# The generally agreed threshold for a good VIF is 2, or 5, or 10 (so both
# "generally" and "agreed" are overstatements here), and as this one is higher,
# it suggests that we may not need all of these data. It is also entirely
# possible to have an *infinite* VIF, in case two variables are perfectly
# correlated (this happens a fair amount with bioclim, in fact).

# For this reason, we will go through an iterative process to get rid of
# variables one by one until the largest VIF is no larger than some threshold.
# Specifically, we get rid of the variable with the largest VIF first. Let's try
# this with the full sample:

vifs = zeros(Float64, length(layers))
for i in eachindex(layers)
    vifs[i] = vif(lm(X[:, setdiff(eachindex(selection), i)], X[:, i]))
end
findmax(vifs)

# This result suggests that this variable, because it has the highest VIF (and
# one that is above our threshold), should be dropped. Repeating the process is
# a good use case for a recursive function:

function stepwisevif(
    layers::Vector{T}, selection=collect(1:length(layers)), threshold::Float64=5.0
) where {T<:SimpleSDMLayer}
    x = hcat([layer[keys(layer)] for layer in layers[selection]]...)
    X = (x .- mean(x; dims=1)) ./ std(x; dims=1)
    vifs = zeros(Float64, length(selection))
    for i in eachindex(selection)
        vifs[i] = vif(lm(X[:, setdiff(eachindex(selection), i)], X[:, i]))
    end
    all(vifs .<= threshold) && return selection
    drop = last(findmax(vifs))
    popat!(selection, drop)
    @info "Variables remaining: $(selection)"
    return stepwisevif(layers, selection, threshold)
end

# This function will operate on a collection of layers, and starting from a
# selection (of indices), iterate until a subset satisfying max(VIF) < threshold
# is found. We can run this function on our entire dataset:

selected_variables_id = stepwisevif(layers)

# Which variables are these?

layernames(WorldClim, BioClim)[selected_variables_id]

# Finally, we can select the variables this process recommends:

layers[selected_variables_id]

# Before we move on -- variable selection, especially stepwise, is not
# necessarilly a *good* practice. Alternatives are model selection,
# dimensionality reduction using *e.g.* PCA, or other methods to remove the
# problematic covariance structure in the data. In Julia, a lot of this can be
# done using the `MultivariateStats` package. This can be an interesting
# exercise: rather than relying on VIF, what would the dimensionsality of the
# results look like with a PCA cutoff at 99% of explained variance?


text:
@testset "960.delete-columns-to-make-sorted-iii.jl" begin
    @test min_deletion_size(["babca","bbazb"]) == 3
    @test min_deletion_size(["edcba"]) == 4
    @test min_deletion_size(["ghi","def","abc"]) == 0
end


text:
using Plots
Plots.gr()

include("simpleModel.jl")
include("demo.jl")

setRNGs(12345)
model, lM = setupSimpleModel(1000, 1.0) # box is [-1.0,1.0], n = 1000
@time b11 = CCSMCBoundaries(model, lM, 2, :BS, 10000); plot(b11)
@time b12 = CCSMCBoundaries(model, lM, 4, :BS, 10000); plot!(b12)
@time b13 = CCSMCBoundaries(model, lM, 8, :BS, 10000); plot!(b13)
@time b14 = CCSMCBoundaries(model, lM, 16, :BS, 10000); plot!(b14)
@time b15 = CCSMCBoundaries(model, lM, 32, :BS, 10000); plot!(b15)
plot(b11, label="N=2")
plot!(b12, label="N=4")
plot!(b13, label="N=8")
plot!(b14, label="N=16")
plot!(b15, label="N=32")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for boxSize = 1.0")
savefig("output/simpleDemoBoundaries1.pdf")

setRNGs(12345)
model, lM = setupSimpleModel(1000, 5.0) # box is [-5.0,5.0], n = 1000
@time b21 = CCSMCBoundaries(model, lM, 2, :BS, 10000); plot(b21)
@time b22 = CCSMCBoundaries(model, lM, 4, :BS, 10000); plot!(b22)
@time b23 = CCSMCBoundaries(model, lM, 8, :BS, 10000); plot!(b23)
@time b24 = CCSMCBoundaries(model, lM, 16, :BS, 10000); plot!(b24)
@time b25 = CCSMCBoundaries(model, lM, 32, :BS, 10000); plot!(b25)
plot(b21, label="N=2")
plot!(b22, label="N=4")
plot!(b23, label="N=8")
plot!(b24, label="N=16")
plot!(b25, label="N=32")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for boxSize = 5.0")
savefig("output/simpleDemoBoundaries2.pdf")

setRNGs(12345)
model, lM = setupSimpleModel(1000, 10.0) # box is [-10.0,10.0], n = 1000
@time b35 = CCSMCBoundaries(model, lM, 32, :BS, 1000); plot(b35)
@time b36 = CCSMCBoundaries(model, lM, 64, :BS, 1000); plot!(b36)
@time b37 = CCSMCBoundaries(model, lM, 128, :BS, 1000); plot!(b37)
@time b38 = CCSMCBoundaries(model, lM, 256, :BS, 1000); plot!(b38)
plot(b35, label="N=32")
plot!(b36, label="N=64")
plot!(b37, label="N=128")
plot!(b38, label="N=256")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for boxSize = 10.0")
savefig("output/simpleDemoBoundaries3.pdf")


text:
# `SymBand` for `FiniteGP`
function Makie.convert_arguments(
    ::Type{P}, x::AbstractVector, gp::FiniteGP
) where {P<:SymBand}
    y, var = mean_and_var(gp)
    Δy = sqrt.(var)
    return convert_arguments(P, x, y, Δy)
end

# `GPSample` for `FiniteGP`
Makie.convert_arguments(::Type{<:GPSample}, x::AbstractVector, gp::FiniteGP) = (x, gp)

# Default fallback: plot mean of GP
function Makie.convert_arguments(P::Makie.PlotFunc, x::AbstractVector, gp::FiniteGP)
    return convert_arguments(P, x, mean(gp))
end

# `Band` with scaling for `FiniteGP`
Makie.used_attributes(::Type{<:Band}, ::AbstractVector, ::FiniteGP) = (:bandscale,)
function Makie.convert_arguments(
    ::Type{<:Band}, x::AbstractVector, gp::FiniteGP; bandscale=1
)
    y, var = mean_and_var(gp)
    Δy = bandscale .* sqrt.(var)
    return (Point2f.(x, y .- Δy), Point2f.(x, y .+ Δy))
end

# default conversions for `FiniteGP` and `AbstractGP`
function Makie.convert_arguments(
    P::Makie.PlotFunc, x::AbstractVector, gp::AbstractGPs.AbstractGP
)
    return convert_arguments(P, gp(x, 1e-9))
end
function Makie.convert_arguments(P::Makie.PlotFunc, gp::FiniteGP)
    return convert_arguments(P, gp.x, gp)
end


text:
using Multilane
using StatsBase
using MCTS
using JLD

# from "Agents" paper
# bsafe = b
# a_thr = 0.1 # from the "General Lane-Changing" paper
# p in [0.2, 1] from Agents Paper, set to 0.5 for all
normal_idm = IDMParam(1.4, 2.0, 1.5, 120/3.6, 2.0, 4.0)
timid_idm = IDMParam(1.0, 1.0, 1.8, 100/3.6, 4.0, 4.0)
aggressive_idm = IDMParam(2.0, 3.0, 1.0, 140/3.6, 1.0, 4.0)
agents_behavior(idm, idx) = IDMMOBILBehavior(idm, MOBILParam(0.5, idm.b, 0.1), idx)

normal = agents_behavior(normal_idm, 1)
timid = agents_behavior(timid_idm, 2)
aggressive = agents_behavior(aggressive_idm, 3)

behaviors = Dict{String, Any}(
    "all_normal" => DiscreteBehaviorSet([normal], WeightVec([1])),
)

for p in linspace(0., 1., 5)
    behaviors[@sprintf("agents_%03d", 100*p)] =
            DiscreteBehaviorSet([normal, timid, aggressive],
                                WeightVec([p, (1-p)/2, (1-p)/2]))
end

# @show behaviors

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())
solvers = Dict{String, Any}(
    "dpw"=>dpws,
    "assume_normal"=>SingleBehaviorSolver(dpws, normal)
)

curve = TestSet(lambda=[1.0, 2.15, 4.64, 10., 21.5, 46.4, 100.,], N=1)

tests = []
for p in linspace(0., 1., 5)
    push!(tests, TestSet(curve, solver_key="dpw",
                         behaviors=@sprintf("agents_%03d", p*100),
                         key=@sprintf("upper_bound_%03d", p*100)))
    push!(tests, TestSet(curve, solver_key="assume_normal",
                         behaviors=@sprintf("agents_%03d", p*100),
                         key=@sprintf("assume_normal_%03d", p*100)))
end

objects = gen_initials(tests, behaviors=behaviors, generate_physical=true)

@show objects["param_table"]
objects["solvers"] = solvers
for is in values(objects["initial_states"])
    @assert nothing !== (is.cars[1].behavior)
end

results = evaluate(tests, objects, parallel=true)

filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results["histories"] = nothing
save(filename, results)
println("results saved to $filename")


text:
# Includes all Julia code except for the file  src/run.jl.
using Statistics
using Random
using Distributed
using Dates
#module ContVarEvolution
include("types.jl")
include("propsel.jl")   # Proportional selection
include("cont_var.jl")  # Primary code for simulation
include("run_cont_var.jl")
include("bin_data.jl")
include("Propsel.jl")
#end
#using Main.ContVarEvolution