text:
"""
ULMFiT - Custom layers

This file contains the custom layers defined for this model:
    AWD_LSTM
    WeightDroppedLSTM
    VarDrop
    PooledDense
"""

import Flux: gate, testmode!, _dropout_kernel

reset_masks!(entity) = nothing
reset_probability!(entity) = nothing

"""
    drop_mask(x, p)

Drop mask generator

This function generates dropout mask for given 'x' with `p` probability
    or
It can be used to generate the mask by giving the shape of the desired mask and probaility
"""
function drop_mask(x, p)
    y = similar(x, size(x))
    Flux.rand!(y)
    y .= Flux._dropout_kernel.(y, p, 1 - p)
    return y
end

drop_mask(shape::Tuple, p; type = Float32) = (mask = rand(type, shape...);mask .= _dropout_kernel.(mask, p, 1 - p))

#################### Weight-Dropped LSTM Cell ######################
"""
Weight-Dropped LSTM Cell

This is an LSTM layer with dropped weights functionality, that is, DropConnect technique

cite this paper to know about DropConnec:
http://yann.lecun.com/exdb/publis/pdf/wan-icml-13.pdf

Moreover this also follows the Vartional DropOut citeria, that is,
the drop mask is remains same for a whole training pass.
This is done by saving the masks in 'maskWi' and 'maskWh' fields
"""
mutable struct WeightDroppedLSTMCell{A, V, S, M}
    Wi::A
    Wh::A
    b::V
    h::S
    c::S
    p::Float64
    maskWi::M
    maskWh::M
    active::Bool
end

function WeightDroppedLSTMCell(in::Integer, out::Integer, p::Float64=0.0;
    init = Flux.glorot_uniform)
    @assert 0 ≤ p ≤ 1
    cell = WeightDroppedLSTMCell(
        init(out*4, in),
        init(out*4, out),
        init(out*4),
        reshape(zeros(Float32, out), out, 1),
        reshape(zeros(Float32, out), out, 1),
        p,
        drop_mask((out*4, in), p),
        drop_mask((out*4, out), p),
        true
    )
    cell.b[gate(out, 2)] .= 1
    return cell
end

function (m::WeightDroppedLSTMCell)((h, c), x)
    b, o = m.b, size(h, 1)
    Wi = m.active ? m.Wi .* m.maskWi : m.Wi
    Wh = m.active ? m.Wh .* m.maskWh : m.Wh
    g = Wi*x .+ Wh*h .+ b
    input = σ.(gate(g, o, 1))
    forget = σ.(gate(g, o, 2))
    cell = tanh.(gate(g, o, 3))
    output = σ.(gate(g, o, 4))
    c = forget .* c .+ input .* cell
    h′ = output .* tanh.(c)
    return (h′, c), h′
end

Flux.@functor WeightDroppedLSTMCell

Flux.trainable(m::WeightDroppedLSTMCell) = (m.Wi, m.Wh, m.b, m.h, m.c)

testmode!(m::WeightDroppedLSTMCell, mode=true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

"""
    WeightDroppedLSTM(in::Integer, out::Integer, p::Float64=0.0)

WeightDroppedLSTM layer

This makes the WeightDroppedLSTMCell stateful, by wrapping the layer in a Recur field.

Defining an instance:

julia> wd = WeightDroppedLSTM(4, 5, 0.3);
"""
function WeightDroppedLSTM(a...; kw...)
    cell = WeightDroppedLSTMCell(a...;kw...)
    hidden = (cell.h, cell.c)
    return Flux.Recur(cell, hidden)
end

"""
    reset!(m)

Resets the h, c parameters of the LSTM Cell.

For more refer [`Flux.reset`](@ref https://fluxml.ai/Flux.jl/stable/models/layers/#Flux.reset!)
"""
function reset!(m)
    try		# to accomodate the definition in previously trained Language Model
        (m.state = (m.cell.h, m.cell.c))
    catch
    	Flux.reset!(m)
    end
end

"""
    reset_masks!(layer)

This is an important funciton since it used to reset the masks
which are saved in WeightDroppedLSTMCell after every pass.

julia> wd = WeightDroppedLSTM()

julia> reset_masks!(wd)
"""
function reset_masks!(wd::T) where T <: Flux.Recur{<:WeightDroppedLSTMCell}
    wd.cell.maskWi = drop_mask(wd.cell.Wi, wd.cell.p)
    wd.cell.maskWh = drop_mask(wd.cell.Wh, wd.cell.p)
    return
end
####################################################################

################## ASGD Weight-Dropped LSTM Layer ##################
"""
    AWD_LSTM(in::Integer, out::Integer, p::Float64=0.0)

Average SGD Weight-Dropped LSTM

This custom layer is used for training the Language model,
instead of standard LSTM layer.
This layer carries two addtional functionality:
    Weight-dropping (DropConnect)
    Averaging of weights

AWD_LSTM is basically a wrapper aroung WeightDroppedLSTM layer,
it has three fields:
    layer : WeightDroppedLSTM layer
    T     : Trigger iteration, to trigger averaging
    accum : After triggring the accumlation of weights is saved here

cite this paper to know more:
https://arxiv.org/pdf/1708.02182.pdf
"""
mutable struct AWD_LSTM
    layer::Flux.Recur
    T::Integer
    accum
end

AWD_LSTM(in::Integer, out::Integer, p::Float64=0.0; kw...) = AWD_LSTM(WeightDroppedLSTM(in, out, p; kw...), -1, [])

Flux.@functor AWD_LSTM

Flux.trainable(m::AWD_LSTM) = (m.layer,)

(m::AWD_LSTM)(in) = m.layer(in)

# To set the trigger point for the AWD_LSTM layer
set_trigger!(t, m) = nothing
set_trigger!(trigger_point::Integer, m::AWD_LSTM) = m.T = trigger_point;

# resets mask of the WeightDroppedLSTM layer contained in AWD_LSTM
reset_masks!(awd::AWD_LSTM) = reset_masks!(awd.layer)

"""
    asgd_step!(i::Integer, layer::AWD_LSTM)

Averaged Stochastic Gradient Descent Step

This funciton performs the Averaging step to the given AWD_LSTM layer,
if the trigger point or trigger iteration is reached.
Arguments:
    i       : current iteration of the training loop
    layer   : AWD_LSTM layer
"""
asgd_step!(i, l) = nothing

function asgd_step!(iter::Integer, layer::AWD_LSTM)
    if (iter >= layer.T) & (layer.T > 0)
        p = get_trainable_params([layer])
        avg_fact = 1/max(iter - layer.T + 1, 1)
        if avg_fact != 1
            layer.accum = layer.accum .+ p
            for (ps, accum) in zip(p, layer.accum)
                ps .= avg_fact*accum
            end
        else
            layer.accum = deepcopy(p)   # Accumulator for ASGD
        end
    end
    return
end
####################################################################

########################## Varitional DropOut ######################
"""
    VarDrop(p::Float64=0.0)

Variational Dropout layer

Unlike standard dropout layer, which applies new dropout mask to every Array ot Matrix passed to it,
this layer saves the mask applied till it is explicitly set to reset mode using 'reset_masks!'.

Usage:
julia> vd = VarDrop()

To reset mask:
julia> reset_masks!(vd)
"""

mutable struct VarDrop{F}
    p::F
    mask
    reset::Bool
    active::Bool
end

VarDrop(p::Float64=0.0) = VarDrop(p, Array{Float32, 2}(UndefInitializer(), 0, 0), true, true)

function (vd::VarDrop)(x)
    vd.active || return x
    if vd.reset
        vd.mask = drop_mask(x, vd.p)
        vd.reset = false
    end
    return (x .* vd.mask)
end

testmode!(m::VarDrop, mode=true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

# method for reseting mask of VarDrop
reset_masks!(vd::VarDrop) = (vd.reset = true)

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

################# Varitional Dropped Embeddings ######################
"""
    DroppedEmbeddings(in::Integer, embed_size::Integer, p::Float64=0.0)

Embeddings with varitional dropout

This struct defines an embedding layer with Varitional Embedding dropout functionality.
Instead of randomly dropping values of embedding matrix,
this layer drops all values of a specific token, in other words,
that token is dropped from the embedding matrix for that particular pass.

Since this follows Variational DropOut criteria, it also saves the drop mask,
which should be reset to new mask explicilty using `reset_masks!` function

# Usage:
It takes input size, embedding size and dropout probability

julia> de = DroppedEmbeddings(1000, 20, 0.4)

To reset mask:

julia> reset_masks!(de)
"""
mutable struct DroppedEmbeddings{A, F}
    emb::A
    p::F
    mask
    active::Bool
end

function DroppedEmbeddings(in::Integer, embed_size::Integer, p::Float64=0.0;
    init = Flux.glorot_uniform)
        de = DroppedEmbeddings{AbstractArray, typeof(p)}(
            init(in, embed_size),
            p,
            drop_mask((in,), p),
            true
        )
    return de
end

function (de::DroppedEmbeddings)(x::AbstractArray, tying::Bool=false)
    dropped = de.active ? de.emb .* de.mask : de.emb
    return tying ? dropped * x : transpose(dropped[x, :])
end

Flux.@functor DroppedEmbeddings

Flux.trainable(m::DroppedEmbeddings) = (m.emb,)

testmode!(m::DroppedEmbeddings, mode=true) =
  (m.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)

function reset_masks!(de::DroppedEmbeddings)
    de.mask = drop_mask(de.mask, de.p)
    return
end
####################################################################

################# Concat Pooling Dense layer #######################
"""
    PooledDense(hidden_sz::Integer, out::Integer, σ = identity)

Concat-Pooled Dense layer

This is basically a modified version of the `Dense` layer.
It takes the `Vector` of outputs of RNN at all time-steps,
then it calculates the mean and max pools for those outputs and
concatenates output RNN at the last time-step with these max and mean pools.
Then this conatenated `Vector` is multiplied with weights and added with bias
and passes through specified activation function.

Usage:
The first argument `hidden_sz` takes length of the ouput of the preceding RNN layer.
Other two arguments are output size and activation function

# Example

julia> pd = PooledDense(40, 20)    # if the output size of the RNN layer is 40 in this case
"""
mutable struct PooledDense{F, S, T}
    W::S
    b::T
    σ::F
end

PooledDense(W, b) = PooledDense(W, b, identity)

function PooledDense(hidden_sz::Integer, out::Integer, σ = identity;
             initW = Flux.glorot_uniform, initb = (dims...) -> zeros(Float32, dims...))
return PooledDense(initW(out, hidden_sz*3), initb(out), σ)
end

Flux.@functor PooledDense

function (a::PooledDense)(x)
    W, b, σ = a.W, a.b, a.σ
    x = cat(x..., dims=3)
    maxpool = maximum(x, dims=3)[:, :, 1]
    meanpool = (sum(x, dims=3)/size(x, 3))[:, :, 1]
    hc = cat(x[:, :, 1], maxpool, meanpool, dims=1)
    σ.(W*hc .+ b)
end

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

"""
get_trainable_params(layers)

This funciton works same as `params` function except for `AWD_LSTM` layer.
While getting `Params` of the `AWD_LSTM` it does not include the `h` and `c` `params` of `AWD_LSTM`.
This is useful while calculating gradients because calculating gradients for `h` and `c` fields
in `AWD_LSTM` is unnecessary here.

# Example:

julia> layers = Chain(DroppedEmbeddings(4,5,0.2),
                    AWD_LSTM(5, 3),
                    Dense(3, 2),
                    softmax
                );
julia> p1 = params(layers);
julia> p2 = get_trainable_params(layers);

julia> length(p1)
8

julia> length(p2)
6

`Params` from all the other layers are included in p2 except for `h` and `c`
"""
function get_trainable_params(layers)
    p = []
    function get_awd_params(awd::AWD_LSTM)
        return [awd.layer.cell.Wi,
        awd.layer.cell.Wh,
        awd.layer.cell.b]
    end
    for layer in layers
        layer isa Array || (layer = [layer])
        for l in layer
            l isa AWD_LSTM && (append!(p, get_awd_params(l)); continue)
            push!(p, l)
        end
    end
    return params(p...)
end


text:
# using VimBindings
import VimBindings: word_next, word_end, word_big_next, Motion, is_punctuation, find_c, line_begin
# import VimBindings: line, TextObject
using Test
# const VB = VimBindings

# include("action.jl")
include("parse.jl")
include("textutils.jl")
include("command.jl")
include("motion.jl")
include("textobject.jl")

@testset "VimBindings.jl" begin
end
@testset "motion.jl" begin
    #               0123456789
    #               |-----|
    buf = IOBuffer("Hello worl")
    motion = word_next(buf)
    @test motion == Motion(0, 6)

    @assert is_punctuation('-')
    buf = IOBuffer("Hello-worl")
    motion = word_next(buf)
    @test motion == Motion(0, 5)

end

@testset "word end" begin

    buf = IOBuffer("using VimBindings")
    @test word_end(buf) == Motion(0, 4)

    motion = word_end(IOBuffer("VimBindings.jl"))
    @test motion == Motion(0, 10)

    motion = word_end(IOBuffer(" VimBindings.jl"))
    @test motion == Motion(0, 11)
end

@testset "big word" begin
    buf = IOBuffer("push!(LOAD_PATH, dirname(file))")
    motion = word_big_next(buf)
    @test motion == Motion(0, 17)
end

# @testset "begin big word" begin
#     buf = IOBuffer("push!(LOAD_PATH, dirname(file))")
#     seek(buf, 20)
#     @test motion == Motion(0, 17)
# end

@testset "find char" begin

    buf = IOBuffer("using VimBindings")
    @test find_c(buf, 'g') == Motion(0, 4)
end

@testset "line diffs" begin

end

@testset "line begin" begin
    s = """
First line
    second line
    third line
"""
    buf = IOBuffer(s)
    seek(buf, 17)
    # @show peek(buf)
    @test line_begin(buf) == Motion(17, 15)

    # in the second line's space
    seek(buf, 11)
    @test line_begin(buf) == Motion(11, 15)

    seek(buf, 4)
    @test line_begin(buf) == Motion(4, 0)
end

ext: jl


text:
for ForceSpaceMatrix in (:MomentumMatrix, :WrenchMatrix)
    @eval struct $ForceSpaceMatrix{A<:AbstractMatrix}
        frame::CartesianFrame3D
        angular::A
        linear::A

        @inline function $ForceSpaceMatrix{A}(frame::CartesianFrame3D, angular::AbstractMatrix, linear::AbstractMatrix) where {A<:AbstractMatrix}
            @boundscheck size(angular, 1) == 3 || throw(DimensionMismatch())
            @boundscheck size(linear, 1) == 3 || throw(DimensionMismatch())
            @boundscheck size(angular, 2) == size(linear, 2) || throw(DimensionMismatch())
            new{A}(frame, angular, linear)
        end
    end
end

"""
$(TYPEDEF)

A momentum matrix maps a joint velocity vector to momentum.

This is a slight generalization of the centroidal momentum matrix
(Orin, Goswami, "Centroidal momentum matrix of a humanoid robot: Structure and properties.")
in that the matrix (and hence the corresponding total momentum) need not be
expressed in a centroidal frame.
"""
MomentumMatrix

for ForceSpaceMatrix in (:MomentumMatrix, :WrenchMatrix)
    @eval begin
        @inline function $ForceSpaceMatrix(frame::CartesianFrame3D, angular::A, linear::A) where {A<:AbstractMatrix}
            $ForceSpaceMatrix{A}(frame, angular, linear)
        end

        @inline function $ForceSpaceMatrix(frame::CartesianFrame3D, angular::A1, linear::A2) where {A1<:AbstractMatrix, A2<:AbstractMatrix}
            $ForceSpaceMatrix(frame, promote(angular, linear)...)
        end

        @inline function $ForceSpaceMatrix{A}(mat::$ForceSpaceMatrix) where A
            $ForceSpaceMatrix(mat.frame, A(angular(mat)), A(linear(mat)))
        end

        function Base.show(io::IO, m::$ForceSpaceMatrix)
            print(io, "$($(string(ForceSpaceMatrix))) expressed in \"$(string(m.frame))\":\n$(Array(m))")
        end

        @inline function transform(mat::$ForceSpaceMatrix, tf::Transform3D)
            @framecheck(mat.frame, tf.from)
            R = rotation(tf)
            Av = R * linear(mat)
            Aω = R * angular(mat) + colwise(×, translation(tf), Av)
            $ForceSpaceMatrix(tf.to, Aω, Av)
        end
    end
end

"""
$(TYPEDEF)

A `Momentum` is the product of a `SpatialInertia` and a `Twist`, i.e.
```math
h^i =
\\left(\\begin{array}{c}
k^{i}\\\\
l^{i}
\\end{array}\\right) =
I^i T^{i, j}_k
```
where ``I^i`` is the spatial inertia of a given body expressed in frame ``i``,
and ``T^{i, j}_k`` is the twist of frame ``k`` (attached to the body) with
respect to inertial frame ``j``, expressed in frame ``i``. ``k^i`` is the
angular momentum and ``l^i`` is the linear momentum.
"""
struct Momentum{T}
    frame::CartesianFrame3D
    angular::SVector{3, T}
    linear::SVector{3, T}

    @inline function Momentum{T}(frame::CartesianFrame3D, angular::AbstractVector, linear::AbstractVector) where T
        new{T}(frame, angular, linear)
    end
end

"""
$(TYPEDEF)

A wrench represents a system of forces.

The wrench ``w^i`` expressed in frame ``i`` is defined as
```math
w^{i} =
\\left(\\begin{array}{c}
\\tau^{i}\\\\
f^{i}
\\end{array}\\right) =
\\sum_{j}\\left(\\begin{array}{c}
r_{j}^{i}\\times f_{j}^{i}\\\\
f_{j}^{i}
\\end{array}\\right)
```
where the ``f_{j}^{i}`` are forces expressed in frame ``i``, exerted at
positions ``r_{j}^{i}``. ``\\tau^i`` is the total torque and ``f^i`` is the
total force.
"""
struct Wrench{T}
    frame::CartesianFrame3D
    angular::SVector{3, T}
    linear::SVector{3, T}

    @inline function Wrench{T}(frame::CartesianFrame3D, angular::AbstractVector, linear::AbstractVector) where T
        new{T}(frame, angular, linear)
    end
end

for ForceSpaceElement in (:Momentum, :Wrench)
    @eval begin
        # Construct with possibly eltype-heterogeneous inputs
        @inline function $ForceSpaceElement(frame::CartesianFrame3D, angular::AbstractVector{T1}, linear::AbstractVector{T2}) where {T1, T2}
            $ForceSpaceElement{promote_type(T1, T2)}(frame, angular, linear)
        end

        """
        $(SIGNATURES)

        Create a $($(string(ForceSpaceElement))) given the angular and linear
        components, which should be expressed in the same frame.
        """
        function $ForceSpaceElement(angular::FreeVector3D, linear::FreeVector3D)
            @framecheck angular.frame linear.frame
            $ForceSpaceElement(angular.frame, angular.v, linear.v)
        end

        function Base.show(io::IO, f::$ForceSpaceElement)
            print(io, "$($(string(ForceSpaceElement))) expressed in \"$(string(f.frame))\":\nangular: $(angular(f)), linear: $(linear(f))")
        end

        function Base.zero(::Type{$ForceSpaceElement{T}}, frame::CartesianFrame3D) where {T}
            $ForceSpaceElement(frame, zero(SVector{3, T}), zero(SVector{3, T}))
        end

        Base.zero(f::$ForceSpaceElement) = zero(typeof(f), f.frame)

        function Random.rand(::Type{$ForceSpaceElement{T}}, frame::CartesianFrame3D) where {T}
            $ForceSpaceElement(frame, rand(SVector{3, T}), rand(SVector{3, T}))
        end

        """
        $(SIGNATURES)

        Transform the $($(string(ForceSpaceElement))) to a different frame.
        """
        @inline function transform(f::$ForceSpaceElement, tf::Transform3D)
            @framecheck(f.frame, tf.from)
            rot = rotation(tf)
            lin = rot * linear(f)
            ang = rot * angular(f) + translation(tf) × lin
            $ForceSpaceElement(tf.to, ang, lin)
        end

        @inline function Base.:+(f1::$ForceSpaceElement, f2::$ForceSpaceElement)
            @framecheck(f1.frame, f2.frame)
            $ForceSpaceElement(f1.frame, angular(f1) + angular(f2), linear(f1) + linear(f2))
        end

        @inline function Base.:-(f1::$ForceSpaceElement, f2::$ForceSpaceElement)
            @framecheck(f1.frame, f2.frame)
            $ForceSpaceElement(f1.frame, angular(f1) - angular(f2), linear(f1) - linear(f2))
        end

        @inline Base.:-(f::$ForceSpaceElement) = $ForceSpaceElement(f.frame, -angular(f), -linear(f))

        function Base.isapprox(x::$ForceSpaceElement, y::$ForceSpaceElement; atol = 1e-12)
            x.frame == y.frame && isapprox(angular(x), angular(y), atol = atol) && isapprox(linear(x), linear(y), atol = atol)
        end
    end
end

# Wrench-specific functions
"""
$(SIGNATURES)

Create a Wrench from a force, ``f``, and the application point of the force, ``r``.
The torque part of the wrench will be computed as ``r \\times f``. The force
and the application point should be expressed in the same frame.
"""
Wrench(application_point::Point3D, force::FreeVector3D) = Wrench(application_point × force, force)


text:
module CmdStanExtract

# Write your package code here.

export extract

"""
# extract(chns::Array{Float64,3}, cnames::Vector{String})

RStan/PyStan style extract

chns: Array: [draws, vars, chains], cnames: ["lp__", "accept_stat__", "f.1", ...]
Output: name -> [size..., draws, chains]
"""
function extract(chns::Array{Float64,3}, cnames::Vector{String})
    draws, vars, chains = size(chns)

    ex_dict = Dict{String, Array}()

    group_map = Dict{String, Array}()
    for (i, cname) in enumerate(cnames)
        sp_arr = split(cname, ".")
        name = sp_arr[1]
        if length(sp_arr) == 1
            ex_dict[name] = chns[:,i,:]
        else
            if !(name in keys(group_map))
                group_map[name] = Any[]
            end
            push!(group_map[name], (i, [Meta.parse(i) for i in sp_arr[2:end]]))
        end
    end

    for (name, group) in group_map
        max_idx = maximum(hcat([idx for (i, idx) in group]...), dims=2)[:,1]
        ex_dict[name] = similar(chns, max_idx..., draws, chains)
    end

    for (name, group) in group_map
        for (i, idx) in group
            ex_dict[name][idx..., :, :] = chns[:,i,:]
        end
    end

    return ex_dict
end

end


text:
using ZigZagBoomerang
using SparseArrays
using StructArrays
using ForwardDiff
using ForwardDiff: Dual, value, partials
const D1 = Dual{Nothing, Float64, 1}
const D𝕏 = typeof(D1.(zero(𝕏)))

function sfwguidtree!(X, guidedsegs, messages, tree::Tree, f, g, θ, Z; apply_time_change=false)
    ll = zeros(D1, tree.n)
    l = rand(tree.lids)
    chance = zeros(tree.n)
    chance[tree.lids] .= 1/length(tree.lids)
    active = zeros(Bool, tree.n)
    active[l] = true
    for i in reverse(eachindex(tree.T))
        i == 1 && continue
        active[tree.Par[i]] = active[i]
        chance[tree.Par[i]] += chance[i]
    end
    @assert chance[1] ≈ 1.0
    for i in eachindex(tree.T)
        i == 1 && continue  # skip root-node (has no parent)
        ipar = tree.Par[i]
        if active[i]
            κ = MSDE.SDEKernel(f, g, messages[i].ts, θ)
            solfw, llnew = MSDE.forwardguiding(κ, messages[i], (X[ipar], 0.0), Z[i-1]; inplace=false, save_noise=true, apply_timechange=apply_time_change)
            X[i] = D𝕏(solfw[end][1:end-1])
            ll[i] = llnew/chance[i] + ll[ipar] * tree.lastone[i]
            guidedsegs[i] = solfw
        else
            ll[i] = ll[ipar] * tree.lastone[i]
        end

    end
    𝐋 = sum(ll[tree.lids])
    X, guidedsegs, ll, 𝐋
end

θinit = [4.5, 0.1], σ0
iters = 5_000
#θs, guidedsegs, frac_accepted = mcmc(tree, Xd, f, g, θlin, θinit; iters=iters)
σprop = 0.05
precisionatleaves=1e-4
dt = 0.01

Q = [i in tree.lids ? WGaussian{(:μ,:Σ,:c)}(Vector(Xd[i]), precisionatleaves*Matrix(I(d)), 0.0) : missing for i in tree.ids]
Q, messages = bwfiltertree!(Q, tree, θlin, dt)

f(u,θ,t) = SVector((tanh.(Diagonal(θ[1]) * M * u))...)  # f(u,θ,t) = Diagonal(θ[1]) * M * u

guidedsegs = Vector{Any}(undef, tree.n) # save all guided segments
X = zeros(D𝕏, tree.n)  # values at nodes
for id in tree.lids
    X[id] = Xd[id]
end

Z = [innov(messages[i].ts) for i ∈ 2:tree.n]

F(x, Z) = sfwguidtree!(X, guidedsegs, messages, tree, f, g, (x, σ0), Z)[4]

t0 = 0.0
x0 = [1.5, 0.5]
N = length(x0)

function mk_∇ϕi(ρ)
    ith = zeros(N)
    Z = [innov(messages[i].ts) for i ∈ 2:tree.n]
    function (x,i)
        Z[:] = [pcn_innov(Z[i], ρ) for i ∈ eachindex(Z)]
        ith[i] = 1
        sa = StructArray{ForwardDiff.Dual{}}((x, ith))
        δ = F(sa, Z).partials[]
        ith[i] = 0
        return 0.1*x[i] - δ
    end
end

v0 = rand((-1.0,1.0), N)
c = 1.0
T = 300.
Γ = sparse(1.0*I(N))
tr, _, (acc, num), cs = @time ZigZagBoomerang.spdmp(mk_∇ϕi(0.99), t0, x0, v0, T, c*ones(N), ZigZag(Γ, 0*x0); adapt=true, structured=false, progress=true)

ts, θs = ZigZagBoomerang.sep(tr)
## trace plots of pars
fig1 = plot(ts, getindex.(θs, 1), color=:blue)
fig2 = plot(ts, getindex.(θs, 2), color=:blue)
plot!(fig1, ts, fill(θ0[1][1], length(ts)), color=:green)
plot!(fig2, ts, fill(θ0[1][2], length(ts)), color=:green)
pl = plot(fig1, fig2, layout = (2, 1), legend = false)
display(pl)
png("zz_trace_theta0.png")

pl


text:
module Expronicon

using MLStyle
using Markdown
using OrderedCollections
using MLStyle.MatchImpl
using MLStyle.AbstractPatterns

export
    # types
    NoDefault, JLExpr, JLFor, JLIfElse, JLMatch,
    JLFunction, JLField, JLKwField, JLStruct, JLKwStruct,
    # analysis
    @expr, @test_expr, compare_expr, compare_vars,
    AnalysisError, is_fn, is_kw_fn, split_function, split_function_head, split_struct,
    split_struct_name, split_ifelse, annotations, uninferrable_typevars, has_symbol,
    is_literal, has_kwfn_constructor, has_plain_constructor,
    # transformations
    no_default, prettify, rm_lineinfo, flatten_blocks, name_only,
    rm_annotations, replace_symbol, subtitute, eval_interp, eval_literal,
    # codegen
    codegen_ast,
    codegen_ast_kwfn,
    codegen_ast_kwfn_plain,
    codegen_ast_kwfn_infer,
    codegen_ast_struct,
    codegen_ast_struct_head,
    codegen_ast_struct_body,
    codegen_match,
    construct_method_plain,
    construct_method_inferable,
    struct_name_plain,
    struct_name_without_inferable,
    # printings
    with_marks, with_parathesis, with_curly, with_brackets, within_line, within_indent,
    with_begin_end, indent, no_indent, no_indent_first_line, indent_print, indent_println

include("patches.jl")
include("types.jl")
include("transform.jl")
include("analysis.jl")
include("codegen.jl")
include("printing.jl")

end


text:
export getStiffAndForce, getInternalForce,
getBodyForce, getEdgeForce, getStrain, getMassMatrix,
getNodes, getGaussPoints, getEdgeGaussPoints, commitHistory

abstract type Continuum end

@doc raw"""
    getBodyForce(elem::Continuum, fvalue::Array{Float64,2})

Returns the body force.

- `fvalue` is a $n_{gauss}\times 2$ matrix, which is ordered the same as
    Gaussian points in the undeformed parent element.

Returns the nodal force due to the body force
$$\int_{e} \mathbf{f}(\mathbf{x})\cdot \delta \mathbf{u}(\mathbf{x}) d \mathbf{x}
  = \int_{e} \mathbf{f}(\mathbf{\xi})\cdot \delta \mathbf{u}(\mathbf{\xi})
  |\frac{\partial \mathbf{x}}{\partial \mathbf{\xi}}| d \mathbf{\xi}$$

!!! todo
    Add force in the deformed domain.
"""
function getBodyForce(elem::Continuum, fvalue::Array{Float64,2})
    n = dofCount(elem)
    fbody = zeros(Float64,n)

    nnodes = length(elem.elnodes)
    for k = 1:length(elem.weights)
        fbody[1:nnodes] += elem.hs[k] * fvalue[k,1] * elem.weights[k]
        fbody[nnodes+1:2*nnodes] += elem.hs[k] * fvalue[k,2] * elem.weights[k]
    end
    return fbody
end

@doc raw"""
    getBodyForce(elem::Continuum, fvalue::Array{Float64, 1})

Returns
```math
\int_A f \delta v dx
```
on a specific element $A$
`fvalue` has the same length as number of Gauss points.
"""
function getBodyForce(elem::Continuum, fvalue::Array{Float64, 1})
    nnodes = length(elem.elnodes)
    fbody = zeros(nnodes)
    for k = 1:length(elem.weights)
        fbody += elem.hs[k] * fvalue[k] * elem.weights[k]
    end
    return fbody
end

@doc raw"""
    getEdgeForce(elem::Continuum, iedge::Float64, fvalue::Array{Float64,2})

Returns the force imposed by boundary tractions.

`fvalue` is a $n_{edge_gauss}\times 2$ matrix, which is ordered the same as the
    Gaussian points in undeformed parent edge element.

```
    The element nodes are ordered as
    #   4 ---- 3             #   4 --7-- 3
    #                        #   8   9   6
    #   1 ---- 2             #   1 --5-- 2
    for porder=1     or          porder=2
    iedge 1, 2, 3, 4 are (1,2), (2,3), (3,4), (4,1)
                    are (1,2,5), (2,3,6), (3,4,7), (4,1,8)
```

Returns the nodal force due to the traction on the iedge-th edge of the element
$$\int_{s} \mathbf{f}(\mathbf{x})\cdot \delta \mathbf{u}(\mathbf{x}) d s
  = \int_{e} \mathbf{f}(\xi)\cdot \delta \mathbf{u}(\xi)
  |\frac{\partial \mathbf{x}}{\partial \xi}| d \xi$$

!!! todo
    This function imposes force in the undeformed domain. Add force in the deformed domain in the future.
"""
function getEdgeForce(elem::Continuum, iedge::Int64, fvalue::Array{Float64,2})
    n = length(elem.elnodes)
    ngp = Int64(sqrt(length(elem.weights)))
    @assert(n == 4 || n == 9)

    n1, n2 = iedge, ((iedge+1)==5 ? 1 : iedge+1)
    loc_id = (n == 4 ? [n1, n2] : [n1, n2, iedge+4])

    x = elem.coords[loc_id,:]

    weights, hs = get1DElemShapeData(x, ngp)

    fedge = zeros(Float64,2n)
    for igp = 1:ngp
        fedge[loc_id] += hs[igp] * fvalue[igp,1] * weights[igp]
        fedge[n .+ loc_id] += hs[igp] * fvalue[igp,2] * weights[igp]
    end
    return fedge

end

"""
    getMassMatrix(elem::Continuum)

Returns the mass matrix and lumped mass matrix of the element `elem`.
"""
function getMassMatrix(elem::Continuum)
    ndofs = dofCount(elem)
    nnodes = length(elem.elnodes)
    mass = zeros(ndofs,ndofs)
    for k = 1:length(elem.weights)
        rho = elem.mat[k].ρ
        mass += [elem.hs[k]*elem.hs[k]' zeros(nnodes, nnodes)
                 zeros(nnodes, nnodes)  elem.hs[k]*elem.hs[k]']  * rho * elem.weights[k]
    end
    lumped = sum(mass, dims=2)
    mass, lumped
end

"""
    getNodes(elem::Continuum)

Alias for `elem.elnodes`
"""
function getNodes(elem::Continuum)
    return elem.elnodes
end

"""
    getGaussPoints(elem::Continuum)

Returns the Gauss quadrature nodes of the element in the undeformed domain
"""
function getGaussPoints(elem::Continuum)
    x = elem.coords'
    gnodes = zeros(length(elem.weights),2)
    for k = 1:length(elem.weights)
        gnodes[k,:] = x * elem.hs[k]
    end
    return gnodes
end

"""
    getEdgeGaussPoints(elem::Continuum, iedge::Int64)

```
    The element nodes are ordered as
    #   4 ---- 3             #   4 --7-- 3
    #                        #   8   9   6
    #   1 ---- 2             #   1 --5-- 2
    for porder=1     or          porder=2
    edge 1, 2, 3, 4 are (1,2), (2,3), (3,4), (4,1)
                    are (1,2,5), (2,3,6), (3,4,7), (4,1,8)
```

Returns the Gauss quadrature nodes of the element on its iedge-th edge in the undeformed domain
"""
function getEdgeGaussPoints(elem::Continuum, iedge::Int64)
    n = length(elem.elnodes)
    ngp = Int64(sqrt(length(elem.weights)))

    @assert(n == 4 || n == 9)

    n1, n2 = iedge, ((iedge+1)==5 ? 1 : iedge+1)
    loc_id = (n == 4 ? [n1, n2] : [n1, n2, iedge+4])

    x = elem.coords[loc_id, :]

    gnodes = zeros(ngp,2)

    _, hs = get1DElemShapeData(x, ngp)

    gnodes = zeros(ngp,2)
    for igp = 1:ngp
        gnodes[igp,:] = x' * hs[igp]
    end

    return gnodes
end

function dofCount(elem::Continuum)
    return 2length(elem.elnodes)
end

"""
    commitHistory(elem::Continuum)

Updates the historic parameters in the material properties.
"""
function commitHistory(elem::Continuum)
    for m in elem.mat
        commitHistory(m)
    end
end


text:
struct MoveGenCommon
    friends::Bitboard
    enemies::Bitboard
    empty::Bitboard
    occ::Bitboard
    ourking::Bitboard
    unpinned::Bitboard
end

function MoveGenCommon(board::Board)
    MoveGenCommon(friendly(board), enemy(board), empty(board), occupied(board), friendly(board) & kings(board), ~pinned(board))
end

function build_pawn_advances!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    if board.turn == WHITE
        build_wpawn_advances!(movestack, board, targets, common)
    else
        build_bpawn_advances!(movestack, board, targets, common)
    end
    return
end

function build_wpawn_advances!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = (pawns(board) & common.friends) & (common.unpinned | file(common.ourking))

    # Single pawn advances
    dests = (pwns << 8) & common.empty & ~RANK_8
    for dest in dests & targets
        # add a normal move
        push!(movestack, Move(dest - 8, dest))
    end

    # Double pawn advances
    doubledests = (dests << 8) & targets & RANK_4
    for dest in doubledests
        # add a double pawn move
        push!(movestack, Move(dest - 16, dest, __DOUBLE_PAWN))
    end
    return
end

function build_bpawn_advances!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = (pawns(board) & common.friends) & (common.unpinned | file(common.ourking))

    # Single pawn advances
    dests = (pwns >> 8) & common.empty & ~RANK_1
    for dest in dests & targets
        # add a normal move
        push!(movestack, Move(dest + 8, dest))
    end

    # Double pawn advances
    doubledests = (dests >> 8) & targets & RANK_5
    for dest in doubledests
        # add a double pawn move
        push!(movestack, Move(dest + 16, dest, __DOUBLE_PAWN))
    end
    return
end

function build_free_pawn_captures!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    if board.turn == WHITE
        build_free_wpawn_captures!(movestack, board, targets, common)
    else
        build_free_bpawn_captures!(movestack, board, targets, common)
    end
    return
end

function build_free_wpawn_captures!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = pawns(board) & common.friends & common.unpinned

    # "Left" pawn captures

    dests = (pwns << 9) & targets & ~FILE_H
    for dest in (dests & ~RANK_8)
        push!(movestack, Move(dest - 9, dest))
    end
    build_promo_internal!(movestack, dests & RANK_8, 9)

    # "Right" pawn captures

    dests = (pwns << 7) & targets & ~FILE_A
    for dest in (dests & ~RANK_8)
        push!(movestack, Move(dest - 7, dest))
    end
    build_promo_internal!(movestack, dests & RANK_8, 7)
    return
end

function build_free_bpawn_captures!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = pawns(board) & common.friends & common.unpinned

    # "Left" pawn captures

    dests = (pwns >> 7) & targets & ~FILE_H
    for dest in (dests & ~RANK_1)
        push!(movestack, Move(dest + 7, dest))
    end
    build_promo_internal!(movestack, dests & RANK_1, -7)

    # "Right" pawn captures

    dests = (pwns >> 9) & targets & ~FILE_A
    for dest in (dests & ~RANK_1)
        push!(movestack, Move(dest + 9, dest))
    end
    build_promo_internal!(movestack, dests & RANK_1, -9)
    return
end

function build_enpass_moves!(movestack::MoveStack, board::Board, common::MoveGenCommon)
    ep_sqr = board.enpass
    if ep_sqr !== zero(UInt8)
        king = square(common.ourking)
        for pwn in (pawnAttacks(!board.turn, ep_sqr) & pawns(board) & common.friends)
            occ = common.occ
            occ &= ~Bitboard(pwn) & ~Bitboard(ep_sqr - 24 + (board.turn.val << 4))
            occ |= Bitboard(ep_sqr)
            if isempty(bishopMoves(king, occ) & bishoplike(board) & common.enemies) && isempty(rookMoves(king, occ) & rooklike(board) & common.enemies)
                push!(movestack, Move(pwn, ep_sqr, __ENPASS))
            end
        end
    end
    return
end

function build_pawn_advance_promos!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    if board.turn == WHITE
        build_wpawn_advance_promos!(movestack, board, targets, common)
    else
        build_bpawn_advance_promos!(movestack, board, targets, common)
    end
    return
end

function build_wpawn_advance_promos!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = pawns(board) & common.friends & (common.unpinned | file(common.ourking))
    dests = (pwns << 8) & targets & RANK_8
    build_promo_internal!(movestack, dests, 8)
    return
end

function build_bpawn_advance_promos!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = pawns(board) & common.friends & (common.unpinned | file(common.ourking))
    dests = (pwns >> 8) & targets & RANK_1
    build_promo_internal!(movestack, dests, -8)
    return
end

function build_promo_internal!(movestack::MoveStack, dests::Bitboard, shift::Integer)
    for dest in dests
        _move_uint16 = Move(dest - shift, dest).val
        push!(movestack, Move(_move_uint16 | __KNIGHT_PROMO))
        push!(movestack, Move(_move_uint16 | __BISHOP_PROMO))
        push!(movestack, Move(_move_uint16 | __ROOK_PROMO))
        push!(movestack, Move(_move_uint16 | __QUEEN_PROMO))
    end
    return
end

function build_pinned_pawn_captures_and_promos!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    pwns = pawns(board) & common.friends & pinned(board)
    king = square(common.ourking)
    for pwn in pwns
        for move_to in (pawnAttacks(board.turn, pwn) & common.enemies & targets)
            if !isempty(Bitboard(pwn) & blockers(king, move_to))
                _move_uint16 = Move(pwn, move_to).val
                if !isempty(Bitboard(move_to) & RANK_18)
                    push!(movestack, Move(_move_uint16 | __KNIGHT_PROMO))
                    push!(movestack, Move(_move_uint16 | __BISHOP_PROMO))
                    push!(movestack, Move(_move_uint16 | __ROOK_PROMO))
                    push!(movestack, Move(_move_uint16 | __QUEEN_PROMO))
                else
                    push!(movestack, Move(_move_uint16))
                end
            end
        end
    end
    return
end

function push_normal!(movestack::MoveStack, sqr_from::Integer, dests::Bitboard)
    for dest in dests
        push!(movestack, Move(sqr_from, dest))
    end
    return
end

# internal function to build knight moves
function build_knight_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    for knight in (knights(board) & common.friends & common.unpinned)
        push_normal!(movestack, knight, knightMoves(knight) & targets)
    end
    return
end

# internal functions to build bishop moves
function build_bishop_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    build_free_bishop_moves!(movestack, board, targets, common)
    if !isempty(pinned(board))
        build_pinned_bishop_moves!(movestack, board, targets, common)
    end
    return
end

function build_free_bishop_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    occ = common.occ
    for bishop in (bishoplike(board) & common.friends & common.unpinned)
        push_normal!(movestack, bishop, bishopMoves(bishop, occ) & targets)
    end
    return
end

function build_pinned_bishop_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    king = square(common.ourking)
    occ = common.occ
    for bishop in (bishoplike(board) & common.friends & pinned(board))
        block_1 = blockers(king, bishop)
        for move_to in (bishopMoves(bishop, occ) & targets)
            if !isempty(Bitboard(move_to) & block_1) || !isempty(Bitboard(bishop) & blockers(king, move_to))
                push!(movestack, Move(bishop, move_to))
            end
        end
    end
    return
end

# internal functions to build rook moves
function build_rook_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    build_free_rook_moves!(movestack, board, targets, common)
    if !isempty(pinned(board))
        build_pinned_rook_moves!(movestack, board, targets, common)
    end
    return
end

function build_free_rook_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    occ = common.occ
    for rook in (rooklike(board) & common.friends & common.unpinned)
        push_normal!(movestack, rook, rookMoves(rook, occ) & targets)
    end
    return
end

function build_pinned_rook_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    king = square(common.ourking)
    occ = common.occ
    for rook in (rooklike(board) & common.friends & pinned(board))
        block_1 = blockers(king, rook)
        for move_to in (rookMoves(rook, occ) & targets)
            if !isempty(Bitboard(move_to) & block_1) || !isempty(Bitboard(rook) & blockers(king, move_to))
                push!(movestack, Move(rook, move_to))
            end
        end
    end
    return
end

# internal function to build king moves
function build_king_moves!(movestack::MoveStack, board::Board, targets::Bitboard, common::MoveGenCommon)
    king = square(common.ourking)
    for move_to in (kingMoves(king) & targets)
        if !isattacked_through_king(board, move_to)
            push!(movestack, Move(king, move_to))
        end
    end
    return
end

function build_castling!(movestack::MoveStack, board::Board, common::MoveGenCommon)
    if cancastlekingside(board)
        build_kingside_castling!(movestack, board, common)
    end
    if cancastlequeenside(board)
        build_queenside_castling!(movestack, board, common)
    end
    return
end

function build_kingside_castling!(movestack::MoveStack, board::Board, common::MoveGenCommon)
    king = square(common.ourking)
    rook = king - 3
    if isempty(common.occ & blockers(king, rook)) && !isattacked(board, king - 1) && !isattacked(board, king - 2)
        push!(movestack, Move(king, king - 2, __KING_CASTLE))
    end
    return
end

function build_queenside_castling!(movestack::MoveStack, board::Board, common::MoveGenCommon)
    king = square(common.ourking)
    rook = king + 4
    if isempty(common.occ & blockers(king, rook)) && !isattacked(board, king + 1) && !isattacked(board, king + 2)
        push!(movestack, Move(king, king + 2, __QUEEN_CASTLE))
    end
    return
end

"""
    gen_quiet_moves!(movestack::MoveStack, board::Board)

Generate all the quiet moves on a `board`, pushing each `Move` to the `MoveStack`.
"""
function gen_quiet_moves!(movestack::MoveStack, board::Board)
    common = MoveGenCommon(board)
    if isdoublecheck(board)
        build_king_moves!(movestack, board, common.empty, common)
        return
    end
    if ischeck(board)
        if isempty(checkers(board) & knights(board))
            targets = blockers(square(common.ourking), square(checkers(board)))
        else
            # if a knight gives check, only *quiet* move is for king to step away
            build_king_moves!(movestack, board, common.empty, common)
            return
        end
    else
        build_castling!(movestack, board, common)
        targets = common.empty
    end
    build_pawn_advances!(movestack, board, targets, common)
    build_king_moves!(movestack, board, common.empty, common)
    build_knight_moves!(movestack, board, targets, common)
    build_bishop_moves!(movestack, board, targets, common)
    build_rook_moves!(movestack, board, targets, common)
    return
end

"""
    gen_noisy_moves!(movestack::MoveStack, board::Board)

Generate all the noisy moves on a `board`, pushing each `Move` to the `MoveStack`.
"""
function gen_noisy_moves!(movestack::MoveStack, board::Board)
    common = MoveGenCommon(board)
    enemies = common.enemies
    if isdoublecheck(board)
        build_king_moves!(movestack, board, enemies, common)
        return
    end
    if ischeck(board)
        targets = checkers(board)
        build_pawn_advance_promos!(movestack, board, blockers(square(common.ourking), square(targets)), common)
    else
        targets = enemies
        build_pawn_advance_promos!(movestack, board, common.empty, common)
    end
    build_free_pawn_captures!(movestack, board, targets, common)
    build_pinned_pawn_captures_and_promos!(movestack, board, targets, common)
    build_enpass_moves!(movestack, board, common)
    build_king_moves!(movestack, board, enemies, common)
    build_knight_moves!(movestack, board, targets, common)
    build_bishop_moves!(movestack, board, targets, common)
    build_rook_moves!(movestack, board, targets, common)
    return
end

"""
    gen_moves!(movestack::Movestack, board::Board)

Generate all the possible moves in the position.
"""
function gen_moves!(movestack::MoveStack, board::Board)
    gen_quiet_moves!(movestack, board)
    gen_noisy_moves!(movestack, board)
    return
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
module Queryable # Octo

using ..Octo: Structured, FromItem, SubQuery, Field, SQLAlias, SQLExtract, SQLFunction, Predicate, SQLKeyword
using Dates: CompoundPeriod
using Dates # DatePeriod, TimePeriod, TimeType

"""
    from(M::Type, alias=nothing)::FromItem
"""
function from(M::Type, alias=nothing)::FromItem
    FromItem(M, alias)
end

function Base.getproperty(clause::FromItem, field::Symbol)
     if field in (:__octo_model, :__octo_alias)
         getfield(clause, field)
     else
         Field(clause, field)
     end
end

"""
    from(query::Structured, alias=nothing)::SubQuery
"""
function from(query::Structured, alias=nothing)::SubQuery
    SubQuery(query, alias)
end

function Base.getproperty(clause::SubQuery, field::Symbol)
     if field in (:__octo_query, :__octo_alias)
         getfield(clause, field)
     else
         Field(clause, field)
     end
end

"""
    as(field::Union{Field, SQLFunction, Predicate}, alias::Symbol)::SQLAlias
"""
function as(field::Union{Field, SQLFunction, Predicate}, alias::Symbol)::SQLAlias
    SQLAlias(field, alias)
end

"""
    extract(field::Union{SQLKeyword, Type{DP}, Type{TP}}, from::Union{Dates.DateTime, DP, TP, Dates.CompoundPeriod})::SQLExtract where DP <: Dates.DatePeriod where TP <: Dates.TimePeriod
"""
function extract(field::Union{SQLKeyword, Type{DP}, Type{TP}}, from::Union{DateTime, <:TimeType, DP, TP, CompoundPeriod})::SQLExtract where DP <: DatePeriod where TP <: TimePeriod
    SQLExtract(field, from)
end

end # module Octo.Queryable


text:
export readtexdir

function ssss(path::String)
    for dir in readdir(path)
        texpath, pdfpath = readtexdir(joinpath(path,dir))

    end
end

function readtexdir(dir::String)
    files = readdir(dir)
    texids = findall(f -> endswith(f,".tex"), files)
    length(texids) == 1 || return
    pdfids = findall(f -> endswith(f,".pdf"), files)
    length(pdfids) == 1 || return

    texpath = joinpath(dir, files[texids[1]])
    pdfpath = joinpath(dir, files[pdfids[1]])
    texstr = readtex(texpath)
    pdtexts = readpdf(pdfpath)
    idx2texts = PDText[]
    for t in pdtexts
        push!(t.tags, "O")
        for _ = 1:sizeof(t.str)
            push!(idx2texts, t)
        end
    end
    pdstr = join(map(t -> t.str, pdtexts))

    for m in eachmatch(r"\\section\**\{(.+)\}", texstr)
        query = replace(string(m[1]), " "=>"")
        r = findfirst(query, pdstr)
        idx2texts[first(r)].tags[1] = "B-section"
        idx2texts[last(r)].tags[1] = "E-section"
    end
    write("a.tex.out", pdtexts)
end

function readtex(path::String)
    str = String(open(read,path))
    str
end


text:
@inline convergenceλ(x, abstol, reltol) = max(abstol, max(one(real(x)), norm(x)) * reltol)

function no_convergence(xn0, xn1, fxn1, abstol, reltol)
    λ = convergenceλ(xn1, abstol, reltol)
    (norm(fxn1) <= λ || (check_approx(xn1, xn0, abstol, reltol) && norm(fxn1) <= cbrt(λ))) ? false : true
end #on ifelse: isn't there going to be a branch anyway when the calling program checks whether this returned true/faslse?
function convergence(xn0, xn1, fxn1, abstol, reltol)
    λ = convergenceλ(xn1, abstol, reltol)
    (norm(fxn1) <= λ || (check_approx(xn1, xn0, abstol, reltol) && norm(fxn1) <= cbrt(λ))) ? true : false
end #on ifelse: isn't there going to be a branch anyway when the calling program checks whether this returned true/faslse?

@inline function initial_conditions(ap::AbstractAsymptoticPosterior, i = profile_ind(ap))
    @inbounds begin
        fx0 = rstar(ap)
        x0 = mode(ap, i)
        x1 = x0 + fx0 * std_estimates(ap,i) # inverse slope
    end
    debug_rootsearch() && @show fx0, x0, x1
    fx1, s = fdf_adjrstar_p(ap, x1, i, Val{true}())
    x0, fx0, x1, fx1, s

end

function linear_search(ap::AbstractAsymptoticPosterior{P,T}, i = profile_ind(ap)) where {P,T}
    x0, fx0, x1, fx1, s = initial_conditions(ap, i)
    ϵ = eps(max(abs(x0),abs(x1)))
    abstol, reltol = √ϵ, ∛ϵ
    # if norm(fx1) <= convergenceλ(x1, abstol, reltol)
    #     return x1#, fx1
    # elseif signbit(fx0) != signbit(fx1)
    #     return custom_bisection(ap, x1, x0, fx1, fx0, true, cbrt(eps()))
    # end
    norm(fx1) <= convergenceλ(x1, abstol, reltol) && return x1

    debug_rootsearch() && @show ((x0 - x1) / (fx1 - fx0), 1 / s)
    debug_rootsearch() && @show x1 + fx1 / s
    debug_rootsearch() && @show x1 + (fx1 * (x0 - x1)) / (fx1 - fx0)
    # x1, x0 = x1 + fx1 * (x0 - x1) / (fx1 - fx0), x1
    x1, x0 = x1 - fx1 / s, x1
    debug_rootsearch() && @show x1

    debug_rootsearch() && @show rstar(ap)
    fx0 = fx1
    fx1, s = fdf_adjrstar_p(ap, x1, i, Val{false}())
    debug_rootsearch() && @show fx1, ap(x1, i)
    # fx1, fx0 = ap(x1, i), fx1

    not_converged = no_convergence(x0, x1, fx1, abstol, reltol)
    # @show not_converged

    # Keep using quadratic approximations until either convergence or signs flip.
    while not_converged && signbit(fx0) == signbit(fx1)
        # x1, x0 = x1 + fx1 * (x0 - x1) / (fx1 - fx0), x1
        x1, x0 = x1 - fx1 / s, x1
        # fx1, fx0 = ap(x1, i), fx1
        fx0 = fx1
        fx1, s = fdf_adjrstar_p(ap, x1, i, Val{false}())
        debug_rootsearch() && @show (fx1, fx0)
        not_converged = no_convergence(x0, x1, fx1, abstol, reltol)
    end
    if not_converged # If signs flipped, switch to Brent
        debug_rootsearch() && @show (x0, fx0, ap(x0, i))
        debug_rootsearch() && @show (x1, fx1, ap(x1, i))
        # reset_state!(ap.state, x1, x0, fx1, fx0)
        # x1 = ap.state.xn1
        x1 = custom_bisection(ap, x1, x0, Φ⁻¹(1 - α(ap) - fx1) + rstar(ap), Φ⁻¹(1 - α(ap) - fx0) + rstar(ap), abstol)
        # fx1 = ap.state.fx1
    end

    x1#, fx1
end

function update(x0, x1, x2, fx0, fx1, fx2)

    # δx01 = (x0 - x1)
    # δx02 = (x0 - x2)
    # δx12 = (x1 - x2)
    # denom = δx01*δx02*δx12

    δx01 = fx2 * (x0 - x1)
    δx02 = fx1 * (x0 - x2)
    δx12 = fx0 * (x1 - x2)
    a = (δx12         - δx02         + δx01)#/denom
    b = (δx12*(x1+x2) - δx02*(x0+x2) + δx01*(x0+x1))/2a#/denom
    c = (δx12*x1*x2   - δx02*x0*x2   + δx01*x0*x1)#/denom

    discriminant = abs2(b) - c/a
    if discriminant > 0 # Would projecting whether the upper or lower root is correct be faster?
        root_discrim = sqrt(discriminant)
        r_u = b + root_discrim
        r_l = b - root_discrim
        x0 = ifelse(norm(r_u - x2) < norm(r_l-x2), r_u, r_l)
    else
        x0 = x2 + fx2 * δx12 / (fx2-fx1)
    end

    # a = (δx12         - δx02         + δx01)#/denom
    # b = (δx12*(x1+x2) - δx02*(x0+x2) + δx01*(x0+x1))/2#/denom
    # c = (δx12*x1*x2   - δx02*x0*x2   + δx01*x0*x1)#/denom

    # discriminant = abs2(b) - c*a
    # if discriminant > 0 # Would projecting whether the upper or lower root is correct be faster?
    #     root_discrim = sqrt(discriminant)
    #     r_u = (b + root_discrim) / a
    #     r_l = (b - root_discrim) / a
    #     x0 = ifelse(norm(r_u - x2) < norm(r_l-x2), r_u, r_l)
    # else
    #     x0 = x2 + fx2 * δx12 / (fx2-fx1)
    # end

    x1, x2, x0
end

function quadratic_search(ap::AbstractAsymptoticPosterior{P,T}, i = profile_ind(ap)) where {P,T}
    x0, fx0, x1, fx1, s = initial_conditions(ap, i)

    abstol, reltol = 4eps(T), cbrt(eps(T))
    if norm(fx1) <= convergenceλ(x1, abstol, reltol)
        return x1#, fx1
    end

    # s = -s
    δf = fx1 - fx0
    δx = x1 - x0
    # sδx = δx*s
    # a = δf + sδx
    # b = (abs2(x1) - abs2(x0))*s + 2x0*δf
    # c = fx1*abs2(x0) + fx0*abs2(x1) - x0*x1*(2fx0 - sδx)
    # # @show (a,b,c)
    # # println("$x0 $c $b $a")
    # # println("$x1 $c $b $a")
    # # @show (x0, x1, fx0, fx1, s)

    # discriminant = abs2(b) -4a*c
    # if discriminant > 0
    #     x2 = (b-sqrt(discriminant)) / 2a
    # else #If there are no solutions to the quadratic problem, we will take a linear step.
        x2 = x1 - fx1 * δx / δf
    # end
    fx2 = ap(x2, i)
    not_converged = no_convergence(x1, x2, fx2, abstol, reltol)

    # Keep using quadratic approximations until either convergence or signs flip.
    while not_converged && signbit(fx2) == signbit(fx1)
        x0,  x1,  x2  = update(x0, x1, x2, fx0, fx1, fx2)
        fx0, fx1, fx2 = fx1, fx2, ap(x2, i)
        not_converged = no_convergence(x1, x2, fx2, abstol, reltol)
    end
    if not_converged # If signs flipped, switch to FalsePosition
        # reset_state!(ap.state, x2, x1, fx2, fx1)
        # x2 = ap.state.xn1
        x2 = custom_bisection(ap, x2, x1, fx2, fx1, true, cbrt(eps()))
        # fx2 = ap.state.fx1
    end

    x2#, fx2
end


text:
using DiffEqSensitivity, Test, QuasiMonteCarlo

function ishi_batch(X)
    A= 7
    B= 0.1
    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
end
function ishi(X)
    A= 7
    B= 0.1
    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
end

function linear_batch(X)
    A= 7
    B= 0.1
    @. A*X[1,:]+B*X[2,:]
end
function linear(X)
    A= 7
    B= 0.1
    A*X[1]+B*X[2]
end

n = 600000
lb = -ones(4)*π
ub = ones(4)*π
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,sampler)

res1 = gsa(ishi,Sobol(order=[0,1,2]),A,B)
res2 = gsa(ishi_batch,Sobol(),A,B,batch=true)

@test res1.S1 ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4
@test res2.S1 ≈ [0.3139335358797363, 0.44235918402206326, 0.0, 0.0] atol=1e-4

@test res1.ST ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4
@test res2.ST ≈ [0.5576009081644232, 0.44237102330046346, 0.24366241588532553, 0.0] atol=1e-4

res1 = gsa(ishi,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000)
res2 = gsa(ishi_batch,eFAST(),[[lb[i],ub[i]] for i in 1:4],n=15000,batch=true)

@test res1.S1 ≈ [0.307599  0.442412  3.0941e-25  3.42372e-28] atol=1e-4
@test res2.S1 ≈ [0.307599  0.442412  3.0941e-25  3.42372e-28] atol=1e-4

@test res1.ST ≈ [0.556244  0.446861  0.239259  0.027099] atol=1e-4
@test res2.ST ≈ [0.556244  0.446861  0.239259  0.027099] atol=1e-4

m = gsa(ishi, Morris(num_trajectory=500000), [[lb[i],ub[i]] for i in 1:4])
@test m.means_star[1,:] ≈ [2.25341,4.40246,2.5049,0.0] atol = 5e-2
@test m.means[1, :] ≈ [-0.416876, -0.0077712, -0.015714,  0.0] atol = 5e-2
@test m.means_star[1,:] ≈ [2.25341,4.40246,2.5049,0.0] atol = 5e-2


text:
module TestCovidModels

using CovidModels

println("running validation tests")
include("validation.jl")
println("running odemodel tests")
include("odemodel.jl")
println("running CT tests")
include("covid.jl")
end


text:
function shufflevector_instrs(W, T, I, two_operands)
    typ = LLVM_TYPES[T]
    vtyp1 = "<$W x $typ>"
    M = length(I)
    vtyp3 = "<$M x i32>"
    vtypr = "<$M x $typ>"
    mask = '<' * join(map(x->string("i32 ", x), I), ", ") * '>'
    v2 = two_operands ? "%1" : "undef"
    M, """
        %res = shufflevector $vtyp1 %0, $vtyp1 $v2, $vtyp3 $mask
        ret $vtypr %res
    """
end
@generated function shufflevector(v1::Vec{W,T}, v2::Vec{W,T}, ::Val{I}) where {W,T,I}
    M, instrs = shufflevector_instrs(W, T, I, true)
    quote
        $(Expr(:meta, :inline))
        Vec(llvmcall($instrs, _Vec{$M,$T}, Tuple{_Vec{$W,$T}, _Vec{$W,$T}}, data(v1), data(v2)))
    end
end
@generated function shufflevector(v1::Vec{W,T}, ::Val{I}) where {W,T,I}
    M, instrs = shufflevector_instrs(W, T, I, false)
    quote
        $(Expr(:meta, :inline))
        Vec(llvmcall($instrs, _Vec{$M,$T}, Tuple{_Vec{$W,$T}}, data(v1)))
    end
end


text:
using Pkg
Pkg.activate(pwd())

using titanic

traindata, trainlabels = loadtrain("../data/train.csv")
testdata = loadtest("../data/test.csv")

w = regression(traindata, trainlabels)

prediction = classify(testdata, w)

savepredictions("../results/predictions-logisticregress.csv", prediction)


text:
"""
    NearestNeighborElement <: NeutralLandscapeMaker

    NearestNeighborElement(; n=3, k=1)
    NearestNeighborElement(n, [k=1])

Assigns a value to each patch using a k-NN algorithmm with `n` initial clusters
and `k` neighbors. The default is to use three cluster and a single neighbor.
"""
@kwdef struct NearestNeighborElement <: NeutralLandscapeMaker
    n::Int64 = 3
    k::Int64 = 1
    function NearestNeighborElement(n::Int64, k::Int64 = 1)
        @assert n > 0
        @assert k > 0
        @assert k <= n
        new(n,k)
    end
end

function _coordinatematrix(mat)
    coordinates = Matrix{Float64}(undef, (2, prod(size(mat))))
    for (i, p) in enumerate(Iterators.product(axes(mat)...))
        coordinates[1:2, i] .= p
    end
    coordinates
end

function _landscape!(mat, alg::NearestNeighborElement)
    fill!(mat, 0.0)
    clusters = sample(eachindex(mat), alg.n; replace=false)
    for (i,n) in enumerate(clusters)
        mat[n] = i
    end
    coordinates = _coordinatematrix(mat)
    tree = KDTree(coordinates[:,clusters])
    if alg.k  == 1
        mat[:] .= nn(tree, coordinates)[1]
    else
        mat[:] .= map(mean, knn(tree, coordinates, alg.k)[1])
    end
    return mat
end


text:
using Euphony
using Base.Test

include("event.test.jl")
include("sequence.test.jl")
include("chord.test.jl")


text:
module bavarian_wood

using JuMP
using Cbc
using CSV

function solve_bavarian_wood()

    n = 10;

    Q = 25;
    h = 5;

    K = 8;
    o = [ 1  1  2  2  3  3  4  4  ];
    d = [ 9  10 9  10 9  10 9  10 ];
    L = [ 28 14 42 35 39 18 19 31 ];

    R = [
        (1,5) (1,7) (1,8) (1,9) (1,10) #=	why this #= =# stuff? because otherwise
    =#  (2,5) (2,7) (2,8) (2,9) (2,10) #=	Julia interprets it as 8 different rows
    =#  (3,5) (3,7) (3,8) (3,9) (3,10) #=	but it is just one row of pairs. I want
    =#  (4,6) (4,7) (4,8) (4,9) (4,10) #=	to display it in 8 rows in the code to
    =#  (5,7) (5,8) (5,9) (5,10)       #=	increase readability. See further dis-
    =#  (6,7) (6,8) (6,9) (6,10)       #=    cussion here: https://github.com/JuliaLang/julia/issues/27533
    =#  (7,8) (7,9) (7,10)             #=
    =#  (8,10)
    ];

    N = [5 6 7 8];

    N_out = [
    #=  1 =# (5, 7, 8, 9, 10)
    #=  2 =# (5, 7, 8, 9, 10)
    #=  3 =# (5, 7, 8, 9, 10)
    #=  4 =# (6, 7, 8, 9, 10)
    #=  5 =# (7, 8, 9, 10)
    #=  6 =# (7, 8, 9, 10)
    #=  7 =# (8, 9, 10)
    #=  8 =# (10)
    #=  9 =# ()
    #= 10 =# ()
    ];

    N_in = [
    #=  1 =# ()
    #=  2 =# ()
    #=  3 =# ()
    #=  4 =# ()
    #=  5 =# (1, 2, 3)
    #=  6 =# (4)
    #=  7 =# (1, 2, 3, 4, 5, 6)
    #=  8 =# (1, 2, 3, 4, 5, 6, 7)
    #=  9 =# (1, 2, 3, 4, 5, 6, 7)
    #= 10 =# (1, 2, 3, 4, 5, 6, 7, 8)
    ];

    c = [
        Inf	Inf	Inf	Inf	280	Inf	530	830	780	1000;
        Inf	Inf	Inf	Inf	245	Inf	495	795	745	965;
        Inf	Inf	Inf	Inf	175	Inf	425	725	675	895;
        Inf	Inf	Inf	Inf	Inf	220	670	680	920	850;
        Inf	Inf	Inf	Inf	Inf	Inf	250	550	500	720;
        Inf	Inf	Inf	Inf	Inf	Inf	450	460	700	630;
        Inf	Inf	Inf	Inf	Inf	Inf	Inf	300	250	470;
        Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	170;
        Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf;
        Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf	Inf;
    ]

    model = Model(with_optimizer(Cbc.Optimizer));

    @variable(model, x[1:n,1:n,1:K] >= 0)
    @variable(model, y[1:n,1:n] >= 0, Int)

    @objective(model, Min, sum( c[i,j]*y[i,j] for (i,j) in R ) +
        sum(
            sum(
                sum(
                    h*x[j,i,k]
                for j in N_in[i] )
            for i in N )
        for k in 1:K ));

    @constraint(model,[k in 1:K,i = o[k]], sum(x[i,j,k] for j in N_out[i]) == L[k]);
    @constraint(model,[k in 1:K,j = d[k]], sum(x[i,j,k] for i in N_in[j]) == L[k]);
    @constraint(model,[(i,j) in R], sum(x[i,j,k] for k in 1:K) <= Q*y[i,j]);

    f = open("bavarian_wood.lp", "w")
    print(f, model)
    close(f)

    #-------
    # SOLVE
    #-------

    optimize!(model)

    #------------------------
    # WRITE SOLUTION TO FILE
    #------------------------

    f = open("bavarian_wood.csv","w");

    for k in 1:K
        print(f,"commodity ",o[k],"-",d[k]," - qty = ",L[k]);
        for (i,j) in R
            if (value(x[i,j,k]) > 0)
                println(f,",route ",i,"-",j,",",value(x[i,j,k]));
            end
        end
    end

    for (i,j) in R
        if (value(y[i,j]) != 0)
            print(f,"route ",i,"-",j," - ",value(y[i,j])," trucks");
            for k in 1:K
                if (value(x[i,j,k]) != 0)
                    println(f,",commodity ",o[k],"-",d[k],",",value(x[i,j,k]));
                end
            end
        end
    end

    close(f);

end

solve_bavarian_wood();

end


text:
include(normpath(joinpath(@__FILE__, "..", "..", "exps", "common.jl")))
using IterTools
using CSV

function runconfig(setup::ScreenSetup,
                   lib::Library,
                   num_runs::Int64,
                   output_dir::Compat.String,
                   suppress_graph::Bool)

    info("Running config")
    methods = [venn, auprc]
    measures = [:inc, :dec, :incdec]
    genetypes = [:sigmoidal, :linear, :all]

    crisprtype = lib.cas9_behavior

    test_method_wrapper = (bc_counts, genes) -> test_methods_snr(bc_counts, genes, methods, measures, genetypes)

    # a little abuse of lexical scoping
    gen_plots = (bc_counts, genes) -> begin
        info("Generating plots")
        new_filename = joinpath(output_dir, "counts.svg")
        nopseudo = bc_counts[(bc_counts[:counts_bin1] .> 0.5) .& (bc_counts[:counts_bin2] .> 0.5), :]
        draw(SVG(new_filename, 10cm, 10cm), plot(nopseudo,
        x=:counts_bin1, y=:counts_bin2, color=:class, Scale.x_log10, Scale.y_log10,
        Theme(highlight_width=0pt), Coord.cartesian(fixed=true),
        Guide.title("$(typeof(crisprtype)) $(typeof(setup))"),
        Guide.xlabel("log counts t0"), Guide.ylabel("log counts t1")))
        new_filename = joinpath(output_dir, "volcano.svg")
        draw(SVG(new_filename, 10cm, 10cm), plot(genes,
        x=:mean_bin2_div_bin1, y=:pvalue_bin2_div_bin1, color=:class, Theme(highlight_width=0pt),
        Guide.title("$(typeof(crisprtype)) $(typeof(setup))"),
        Guide.xlabel("mean log2 fold change"), Guide.ylabel("-log10 pvalue")))
        test_methods_snr(bc_counts, genes, methods, measures, genetypes)
    end

    if suppress_graph
        funcs = [test_method_wrapper for i in 1:num_runs]
    else
        eval(:(using Gadfly))
        funcs = [gen_plots, [test_method_wrapper for i in 1:(num_runs - 1)]...]
    end

    # setup up runs
    runs = [(setup, lib, funcs[i], i) for i in 1:num_runs]

    results = pmap(args -> run_exp(args[1], args[2], args[3]; run_idx=args[4]), runs)

    info("Analyzing results")

    results = DataFrame(permutedims(hcat(results...), [2, 1]))
    mean_snr = mean(results[:x19]./results[:x20])
    std_snr = std(results[:x19]./results[:x20])
    snr_result = "SNR score = $(round(mean_snr, 3)) +/- $(round(std_snr, 3))"
    delete!(results, [:x19, :x20])
    hierarchy = vcat([hcat(item...) for item in IterTools.product(map(Symbol, methods), measures, genetypes)]...)
    new_names = [[:method, :measure, :genetype, :score]...; fieldnames(setup)...; :run_idx]
    results = construct_hierarchical_label(hierarchy, results, new_names)

    venn_result, auprc_result = "", ""

    grouped_info = by(results, [:method, :measure, :genetype]) do grouped_df
        n = size(grouped_df, 1)
        mean_score = mean(grouped_df[:score])
        std_score = std(float.(grouped_df[:score]))
        conf_int = 2.58 * std_score./sqrt(n)

        if (grouped_df[1, :measure] == :incdec) &&
           (grouped_df[1, :genetype] == :all)

            disp_score = round(mean_score, 3)
            conf_max = round(mean_score + conf_int, 3)
            conf_min = round(mean_score - conf_int, 3)

            if (grouped_df[1, :method] == :venn)
                venn_result = "Venn score = $disp_score, 95% conf int ($conf_min, $conf_max)"
            else
                auprc_result = "AUPRC score = $disp_score, 95% conf int ($conf_min, $conf_max)"
            end
        end

        DataFrame(
            std_score = std_score,
            mean_score = mean_score,
            conf_max = mean_score + conf_int,
            conf_min = mean_score - conf_int,
            n = n
        )
    end
    info("Saving results in $output_dir")
    CSV.write(joinpath(output_dir, "results_table.csv"), grouped_info)
    println("\nQuick results:\n##############\n$venn_result\n$auprc_result\n$snr_result")
end


text:
#  File for random projection tree implementation
# we instantiate with Tree and node with T=Vector{Float64}

# interesting distance are
#                        Euclidean     L2
#                        Cityblock     L1
#                        Jaccard       fuzzy nucleotide : 1 - sum min (a,b)/ sum max(a,b)
#                        CosineDist
# they are instantiated by calling Euclidean(),  ...
#
#  TODOs:
#     diameter estimation with kcover algorithm?
#     optimisation of threshold used in splitting
#

const rptreeDebugLevel = 1
debuglock=Threads.ReentrantLock()

using Distances
using Match
using Statistics
using LinearAlgebra
using Random
using Printf
using Distributed

include("Tree.jl")

const splitDiam = Int64(1)
const splitProj = Int64(2)
const splitNull = Int64(3)   # when a node is not splitted (leaf)

export
    RPTree,
    RPTreeArg,
    RPTProjParams,
    RPTreeEvent,
    randomProjection,
    analyzeSplittingInfo

"""
KeyVector

The base data of node is a vector indexed by its rank in original data
so it is a : Dict{Int64,Vector{Float64}}.

We store in a node the Dictionary associating the rank
in original data and the corresponding data vector.

"""
const KeyVector = Dict{Int64,Vector{Float64}}

"""

# RPTDiamSplit

This struct stores data necessary to describe a split by diameter

FIELDS
-----

* pivot : the vector aroud which we construct a ball to split
* radius : distance to pivot above which we go to right, under which we go to left
 """
mutable struct RPTDiamSplit
    pivot::Vector{Float64}
    radius::Float64
    #
    function RPTDiamSplit(center::Vector{Float64})
        new(center, -1.)
    end
end

"""
# RPTProjParams

This struct stores the data along which data are projected.

FIELDS
------

- projray : ray on which do project data
- cosineVector : a vector containing for each data the cosine with projray.

Data are split according to position of cosine above or below median of cosine)
"""
mutable struct RPTProjParams
    projray::Vector{Float64}
    cosineVector::Vector{Float64}
    #
    function RPTProjParams(ray::Vector{Float64}, cosVec::Vector{Float64})
        new(ray,cosVec)
    end
end

"""

# RPTreeEvent

Event on Node of Tree
We store the event as it enables propagating a new data across the tree if we
want to do any learning

 FIELDS

. split : the split event coded as constants `splitDiam` or `splitProj` or `splitNull`
. diameters : a vector of size 2 containing mean diameter and then max diameter
. projEvent  field of struct Union{RPTProjParams, Nothing}
. diamEvent  field of struct Union{RPTDiamSplit, Nothing}

  METHODS

    1. `function RPTreeEvent(s::Int64, diameters::Vector{Float64}, pivot::Union{RPTDiamSplit,Nothing})`
    constructor for split on diameter criteria

    2. `function RPTreeEvent(s::Int64, diameters::Vector{Float64}, projData::Union{RPTProjParams, Nothing})`
        constructor for split on projection against a random Vector
        as diameters are always estimated (the struct of split depends upon computed diameters),
        they always are in the constructor

    3. `function RPTreeEvent(diameters::Vector{Float64})`
         a constructor used to just store diameter estimation at leaves. Used with event splitNull
"""
mutable struct RPTreeEvent
    split::Int64
    # first mean diam then max diam
    diameters::Vector{Float64}
    projEvent::Union{RPTProjParams, Nothing}
    diamEvent::Union{RPTDiamSplit, Nothing}
    # constructor for split by diameter
    function RPTreeEvent(s::Int64, diameters::Vector{Float64}, diamEvent::Union{RPTDiamSplit, Nothing})
        # could check that s is splitDiam or splitProj
        new(s, diameters, nothing, diamEvent)
    end
    #
    function RPTreeEvent(s::Int64, diameters::Vector{Float64}, projData::Union{RPTProjParams, Nothing})
        new(s, diameters, projData)
    end
    #
    function RPTreeEvent(diameters::Vector{Float64})
        new(splitNull, diameters, nothing, nothing)
    end
end

"""

# RPTreeArg
 collects the argument for tree growing

 FIELDS

 *  D : the metric to use to compute all distances
 *  depth : The maximum depth of the tree
 *  threshold : the ratio between the mean (meanD) and the max (maxD) of distances between 2 objetcs in node.

    and the  maxDiameter and mean distance between 2 objects above which
    we split by diameter.
    If maxD / meanD > threshold we split by diameter to reach some sphericity in the node
    else we split by projection.
    So with threshold = 1. we always use split by diameter rule and the higher
    threshold the more split nodes by the Projection rule.

    The function analyzeSplittingInfo can be used to get information on the splitting after a run
    and adjust threshold depending on the information obtained.

 CONSTRUCTORS

 RPTreeArg(D::Distances.SemiMetric, depth::Int64, thresh::Float64)

"""
mutable struct RPTreeArg
    D::Distances.SemiMetric
    # the depth of tree
    depth::Int64
    # the ratio between maxDiameter and mean distance between 2 objects above which we split by diameter
    threshold::Float64
    #
    RPTreeArg(D::Distances.SemiMetric, depthArg::Int64, thresholdArg::Float64) = new(D, depthArg, thresholdArg)
end

"""
# RPTree

The struct RPTree stores the list of event that occurred during tree construction.
The parameters used to grow the tree are stored in field argument.

 FIELDS

* treedata : the Tree
* argument : the parameters describing tree growing
* eventDict : the event that occured in each node.
* treelock : a lock to manipulate eventDict

"""
mutable struct RPTree
    treedata::Tree{KeyVector, RPTreeEvent}
    argument::RPTreeArg
    eventDict::Dict{TreeNode{KeyVector,RPTreeEvent}, RPTreeEvent}
    #
    treelock::Threads.ReentrantLock
    #
    function RPTree(argumentArg::RPTreeArg, vdata::Array{Vector{Float64},1})
        eventDict = Dict{TreeNode{KeyVector,RPTreeEvent}, RPTreeEvent}()
        # transform data into indexed data
        keyData = KeyVector()
        for i in 1:length(vdata)
            keyData[i] = vdata[i]
        end
        new(Tree{KeyVector, RPTreeEvent}(keyData), argumentArg, eventDict, Threads.ReentrantLock())
    end
end

"""
diameterEstimation

compute some  estimates of distances between pairs of element in node

returns
* mean distance between pairs,
* max distance of pair
* key of data element with min sum of distances to others.
"""
function diameterEstimation(D::Distances.SemiMetric, node::TreeNode{KeyVector} , fraction::Float64)
    #
    # diameter max and mean estimation
    # We sample elements and compute distances between all couples
    #
    nb=length(node.data)
    if nb <= 1
        return 0.,0.,0
    end
    allKeys = collect(keys(node.data))
    keptI = filter(x -> (rand() < fraction) , allKeys)
    # could not use map!(node -> evaluate(D, node.data, center) , distances, collection)
    # as we do not use all indices. ??
    nbkept=length(keptI)
    costs=zeros(nbkept,nbkept)   # to store sum of distance of each kept index to all others.
    maxByI=zeros(nbkept)
    maxDiameter=0.
    meanDiameter=0.
    nbcouple = 0
    for i in 2:nbkept
        distancesToI = Vector{Float64}(undef,i-1)    # in fact this is cost free versus one maximal allocation
                                               # and maximum taken on a subrange
        map!(j->evaluate(D, node.data[keptI[i]], node.data[keptI[j]]), distancesToI, 1:i-1)
        costs[1:i-1,i]=distancesToI
    end
    Fcosts=Symmetric(costs,:U)
    maxDiameter = maximum(Fcosts)
    meanDiameter = sum(Fcosts)
    # as we have 0 on diagonal, we divide by full matrix minus diagonal
    nbcouple = (nbkept-1) * nbkept
    meanDiameter = meanDiameter/nbcouple
    costmin,imin = findmin(costs)
    #
    return meanDiameter, maxDiameter , keptI[imin]
end

# the following is faster than the preceding as soon as we have 1000 nodes after applying fraction
# for vectors of lenght 10000. So we need 10**5 distances to use parallelism
#

# in the following we should use remotecall_fetch with a

function diameterEstimationByRange(D::Distances.SemiMetric, node::TreeNode{KeyVector,RPTreeEvent},
                                   keptI::Array{Int64,1}, range::UnitRange{Int64})
    maxByTask = 0.
    sumByTask = 0.
    nbkept = length(keptI)
    costs=zeros(nbkept,nbkept)   # to store sum of distance of each kept index to all others.
    ikeys = keys(node.data)
    for i in range
        distancesToI = Vector{Float64}(undef, i-1)
        map!(j->evaluate(D, node.data[ikeys[keptI[i]]], node.data[ikeys[keptI[j]]]), distancesToI, 1:i-1)
        costs[1:i-1,i]=distancesToI
    end
    Fcosts=Symmetric(costs,:U)
    maxByTask = maximum(Fcosts)
    # sum and convert the 1*nbkept Array{Float64,2} to a Vector{Float64} !!
    colCost=convert(Vector{Float64},sum(Fcosts,1)[1,:])
    sumByTask = sum(colCost)
    #
    return sumByTask, maxByTask, colCost
end

#
# The following is not efficient
#

function diameterEstimation2Threads(D::Distances.SemiMetric, node::TreeNode{KeyVector,RPTreeEvent} , fraction::Float64)
    #
    # diameter max and mean estimation
    # We sample elements and compute distances between all couples
    #
    nb=length(node.data)
    keptI = filter(x -> (rand() < fraction) , 1:nb)
    # could not use map!(node -> evaluate(D, node.data, center) , distances, collection)
    # as we do not use all indices. ??
    nbkept=length(keptI)
    maxDiameter=0.
    meanDiameter=0.
    nbcouple = 0
    res = Vector{ Tuple{Float64,Float64, Vector{Float64}} }(undef, 2)
    # compute m so that the 2 threads have the same amount of work!!!
    m=round(Int64, nbkept/sqrt(2.))
    myRanges=Array{UnitRange{Int64}}(2)
    myRanges[1] = 2:m
    myRanges[2] = m+1:nbkept
    Threads.@threads  for i in 1:2
        res[1] = diameterEstimationByRange(D, node, keptI, 2:m)
        res[2] = diameterEstimationByRange(D, node, keptI, m+1:nbkept)
    end
    meanDiameter = res[1][1] + res[2][1]        # sum of means
    maxDiameter = max(res[1][2], res[2][2])     # max of max
    nbcouple = (nbkept-1) * nbkept
    meanDiameter = meanDiameter/nbcouple
    # compute costs
    costs=res[1][3] .+ res[2][3]
    costmin,imin = findmin(costs)
    #
    return meanDiameter, maxDiameter, imin
end

function splitNodeByDiameter(D::SemiMetric, node::TreeNode{KeyVector, RPTreeEvent}, pivot::Int64)
    #
    #
    # compute mean
    #
    nb=length(node.data)
    center=node.data[pivot]
#    center=node.data[1]
#    for i in 2:nb
#        center = center .+ node.data[i]
#    end
#    center = center ./ nb
    #
    # compute median to distances to center
    #
    distances = zeros(Float64,nb)
    vectors = collect(values(node.data))
    map!(vec -> evaluate(D, vec, center) , distances, vectors)
    medianDist=median(distances)
    #
    leftData = KeyVector()
    rightData = KeyVector()
    # it is documented that keys and values are in the same order !!
    # So the i-th distance correspond to the i-th values and to the i-th key
    vindex = collect(keys(node.data))
    for i in 1:length(vindex)
        if distances[i] <= medianDist
            # go to left
            leftData[vindex[i]] = vectors[i]
        else
            # go to right
            rightData[vindex[i]] = vectors[i]
        end
    end
    #
    # We can create children
    #
    leftNode = TreeNode{KeyVector, RPTreeEvent}(node, leftData)
    rightNode = TreeNode{KeyVector, RPTreeEvent}(node, rightData)
    #
    node.data = KeyVector()   # free .... all data have been pushed down
    #
    return leftNode,rightNode,medianDist
end # end of splitByDiameter

# generate random direction too split data in node

function generateRandomDirection2(node::TreeNode{KeyVector})
    item, state = iterate(node.data)
    dim = length(item[2])
    #
    direction = zeros(dim)
    randomNumbers = Vector{Float64}(undef, dim+1)
    #
    Random.rand!(randomNumbers)
    i = 1
    while i < dim
        s = sqrt(-2 * log(randomNumbers[i]))
        t = 2*randomNumbers[i+1]*pi
        g1 = s * cos(t)
        g2 = s * sin(t)
        direction[i] = g1
        if i+1 <= dim
            direction[i+1] = g2
            i = i+2
        else
            break
        end
    end    # end while
    return direction
end

function generateRandomDirection(node::TreeNode{KeyVector, RPTreeEvent})
    item, state = iterate(node.data)
    dim = length(item[2])
    #
    direction = zeros(dim)
    randomNumbers = Vector{Float64}(undef, dim+1)
    #
    Random.rand!(randomNumbers)
    halfsize = round(Int64,dim/2)
    @inbounds @simd for i in 1:halfsize
        j=2*i-1
        s = sqrt(-2 * log(randomNumbers[j]))
        t = 2*randomNumbers[j+1]*pi
        g1 = s * cos(t)
        g2 = s * sin(t)
        direction[j] = g1
        direction[j+1] = g2
    end
    if dim%2 != 0
        s = sqrt(-2 * log(randomNumbers[dim]))
        t = 2*randomNumbers[dim+1]*pi
        g1 = s * cos(t)
        direction[dim] = g1
    end
#
    return direction
end

function splitNodeByProjection(node::TreeNode{KeyVector, RPTreeEvent})
    # generate random normal direction
    direction = generateRandomDirection(node)
    # compute scalar product of each element in node with direction. we have the dot function
    nb=length(node.data)
    scalarProduct = Vector{Float64}(undef, nb)
    # possibly iterate to avoid the collect ?
    map!(x-> dot(x,direction), scalarProduct, collect(values(node.data)))
    medVal=median(scalarProduct)
    # dispatch
    leftData = KeyVector()
    rightData = KeyVector()
    indexes = collect(keys(node.data))
    #
    projEvent = RPTProjParams(direction,scalarProduct)
    #
    for i in 1:nb
        if scalarProduct[i] <= medVal
            # go to left
            leftData[indexes[i]] = node.data[indexes[i]]
        else
            # go to right
            rightData[indexes[i]] = node.data[indexes[i]]
        end
    end
    #
    # We can create children
    #
    leftNode = TreeNode{KeyVector, RPTreeEvent}(node, leftData)
    rightNode = TreeNode{KeyVector, RPTreeEvent}(node, rightData)
    #
    #
    node.data = KeyVector()   # free .... all data have been pushed down
    #
    return leftNode,rightNode,projEvent
end

# splitNodeDiamAndProjection
# Here we must choose the way we do the splitting

function splitNodeDiamAndProjection(rptarg::RPTreeArg , node::TreeNode{KeyVector})
    #
    D = rptarg.D
    threshold = rptarg.threshold
    split = splitNull
    private = nothing
    #
    if length(node.data) == 0
        exit(1)
    end
    #  do an estimation of diameter
    sampleSize = 512.
    fraction = min(1., sampleSize/length(node.data))
    #
#    medDiam, maxDiam, pivot  =  diameterEstimationThreaded(D, node, fraction)
#     medDiam, maxDiam, pivot  =  diameterEstimation2Threads(D, node, fraction)
    medDiam, maxDiam, pivot  = diameterEstimation(D, node, fraction)
    #
    # to store node properties
    #
    diameters=Vector{Float64}(undef, 2)
    diameters[1] = medDiam
    diameters[2] = maxDiam
    # depending on node sphericity we split one way or the other
    if maxDiam/medDiam >= threshold
        # as splitNodeByDiameter clears its data we must get at pivot before the split
        diamEvent = RPTDiamSplit(node.data[pivot])
        leftNode,rightNode, radius = splitNodeByDiameter(D, node, pivot)
        diamEvent.radius = radius
        split = splitDiam
        private = RPTreeEvent(split , diameters, diamEvent)
    elseif medDiam > 0
        leftNode,rightNode, projparams = splitNodeByProjection(node)
        split = splitProj
        private = RPTreeEvent(split , diameters, projparams)
        # else nothing , means all element are equal
    end
    #
    node.private = private
    #
    return leftNode,rightNode
end # end of splitNodeDiamAndProjection

# once we got in serial mode we stick to it

function splitNodeSerial(rptarg::RPTreeArg, node::TreeNode)
    #
    depth = rptarg.depth
    if rptreeDebugLevel > 1
        lock(debuglock)
        @printf stdout "\n splitNodeSerial process id : %d thread %d depth %d " myid() Threads.threadid()  node.depth
        unlock(debuglock)
    end
    #
    if node.depth < depth && length(node.data) > 0
        leftNode,rightNode = splitNodeDiamAndProjection(rptarg ,node)
        # possibly a lock on tree , but we have a lock on node in which we insert
        insertChildInNode(node, leftNode)
        insertChildInNode(node, rightNode)
        #
        splitNodeSerial(rptarg, leftNode)
        splitNodeSerial(rptarg, rightNode)
    end
    if rptreeDebugLevel > 1
        lock(debuglock)
        @printf stdout "\n splitNodeSerial process id : %d thread %d depth %d end " myid() Threads.threadid()  node.depth
        unlock(debuglock)
    end
    return node
end

# The following can replace splitNodeSerial , it just do one step by hand and call 2 threads for
# further splitting
#
function splitNodeThreaded(rptarg::RPTreeArg, node::TreeNode)
    depth = rptarg.depth
    if rptreeDebugLevel > 1
        lock(debuglock)
        @printf stdout "\n splitNodeTheaded process id : %d depth %d " myid()   node.depth
        unlock(debuglock)
    end
    #
    if node.depth < min(depth, 3) &&  Threads.nthreads() > 1
        leftNode,rightNode = splitNodeDiamAndProjection(rptarg ,node)
        # possibly a lock on tree , but we have a lock on node in which we insert
        insertChildInNode(node, leftNode)
        insertChildInNode(node, rightNode)
        # now a threaded loop
        nodesToSplit = Vector{TreeNode}(undef, 2)
        nodesToSplit[1] = leftNode
        nodesToSplit[2] = rightNode
        Threads.@threads for i in 1:2
            splitNodeThreaded(rptarg,nodesToSplit[i])
        end
    else
        splitNodeSerial(rptarg, node)
    end
    if rptreeDebugLevel > 0
        lock(debuglock)
        @printf stdout "\n splitNodeTheaded process id : %d depth %d end " myid()   node.depth
        unlock(debuglock)
    end
    return node
end

function splitNodeParallel(rptarg::RPTreeArg, node::TreeNode)
    #
    depth = rptarg.depth
    @debug "\n  splitNodeParallel process id : %d depth %d " myid()  node.depth
    #
    if node.depth < depth && length(node.data) > 0
        leftNode,rightNode = splitNodeDiamAndProjection(rptarg, node)
        #
        if nworkers() >= 2 && node.depth < 1
            resLeftSubTree = Future()
            resRightSubTree = Future()
            @sync begin
                resLeftSubTree =  @spawnat :any splitNodeThreaded(rptarg, leftNode)
                resRightSubTree = @spawnat :any splitNodeThreaded(rptarg, rightNode)
            end
            leftNode2 = fetch(resLeftSubTree)
            rightNode2 = fetch(resRightSubTree)
            insertChildInNode(node, leftNode2)
            insertChildInNode(node, rightNode2)
        else
            insertChildInNode(node, leftNode)
            insertChildInNode(node, rightNode)
            splitNodeSerial(rptarg, leftNode)
            splitNodeSerial(rptarg, rightNode)
        end
    end
    return node
end

#
# The driver method that grow a tree up to argument depth
#

"""

# function randomProjection(rptree::RPTree)

The driver method that grow a tree up to argument depth.

The splitting of nodes is parallelised:
- if ``nworkers() > 1`` the first node is split the left and right children are
    are affected to two independents tasks via @spawnat : any.

- Furthermore if ``Threads.nthreads() > 1`` the next generation of children
    is split in 2 independent threads, thus giving 4 independants sub trees generated
    in parallel.

Returns:
-------
* an Array{Vector{Float64},1} storing the center of leaves data vectors.

NOTA : The method cannot be run twice on the same rptree. A new RPTree
must be initialized.
"""
function randomProjection(rptree::RPTree)
    #
    node=rptree.treedata.root
    nbcpus = nworkers()
    @info "found nb workers " nbcpus
    #
    if nbcpus <= 1 && Threads.nthreads() <= 1
        @info "going serial"
        splitNodeSerial(rptree.argument, node)
    elseif nbcpus <= 1 && Threads.nthreads() > 1
        @info "going to threads"
        splitNodeThreaded(rptree.argument, node)
    else
        @info "going parallel"
        node = splitNodeParallel(rptree.argument, node)
        rptree.treedata.root = node
    end
    # now we must iterate through leaves of tree, compute centers and send that to clustering
    @info "randomProjection: collecting leaves of tree \n"
    leafCenters=Array{Vector{Float64},1}()
    leaf = getFirstLeftLeaf(rptree.treedata)
    while leaf !== nothing
#        @printf stdout "\n depth %d rank in Parent %d" leaf.depth leaf.rankInParent
        realdata = values(leaf.data)
        center = sum(realdata)
        center = center ./ length(realdata)
        push!(leafCenters,center)
        leaf=getNextLeafRight(rptree.treedata, leaf)
    end
    @printf stdout "\n randomProjection, collected nb leaves = %d \n" length(leafCenters)
    #
    return leafCenters
end

function getFirstLeftNodeAtDepth(rptree::RPTree, depth::Int64)
    if depth > rptree.argument.depth
        return NullNode
    end
    #
    node = rptree.treedata.root
    while node.depth < depth
        if length(node.children) > 0
            node = node.children[1]
        else
            return NullNode
        end
    end
    if node.depth == depth
        return node
    end
    #
    return NullNode
end

"""

    function fillSplittingInfo(rptree::RPTree)

fills `rptree.eventDict` with events describing how each splitted node was split
For leaves computes median and max diameter and associate an event with splitNull

"""
function fillSplittingInfo(rptree::RPTree)
    node = rptree.treedata.root
    sampleSize = 512
    #
    nbseen = 0
    nbproj = 0
    nbdiam = 0
    while node !== nothing
        if rptreeDebugLevel > 1
            @printf stdout " \n dump de noeud : %d, depth %d"  nbseen  node.depth
            dumpPos(node)
        end
        nbseen += 1
        if length(node.children) > 0
            rptree.eventDict[node] = node.private
        else
            # we have a leaf
            fraction = min(1., sampleSize/length(node.data))
            res = diameterEstimation(rptree.argument.D, node, fraction)
            diameters=[res[1] , res[2]]
            # the following just store diameters but do not register an event so that analyzeSplittingInfo
            # will not be polluted
            event = RPTreeEvent(diameters)
            node.private = event
            rptree.eventDict[node] = event
        end
        node = getDepthFirstNextRight(rptree.treedata, node)
    end
    @printf stdout "\n nbstored = %d \n" length(rptree.eventDict)
end     # end of fillSplittingInfo

"""
# function analyzeSplittingInfo(rptree::RPTree)

This function is used a posteriori to get statistics on how nodes were split.
    Presently it computes statistics on leaves diameter.
    It returns :
- leaves = Array{TreeNode{KeyVector,RPTreeEvent}}()
- leafDiameters : a Array{Float64,1} containing for each leaf its final diameter.

"""
function analyzeSplittingInfo(rptree::RPTree)
    #
    @info "in analyzeSplittingInfo"
    #
    leafDiameters = Array{Float64,1}()
    leaves = Vector{TreeNode{KeyVector,RPTreeEvent}}()
    nbsplitDiam = 0
    nbsplitProj = 0
    sampleSize = 512
    nbleaf = 0
    #
    if length(rptree.eventDict) == 0
        fillSplittingInfo(rptree)
    end
    #
    res_iter = iterate(rptree.eventDict)
    while res_iter !== nothing
        item = res_iter[1]
        #        item[1] is a node, item[2] a RPTreeEvent
        diams  = item[2].diameters
        if item[2].split == splitDiam
            nbsplitDiam = nbsplitDiam+1
        elseif item[2].split == splitProj
            nbsplitProj = nbsplitProj+1
        elseif item[2].split == splitNull
            # we have a leaf
            nbleaf = nbleaf + 1
            leaf = item[1]
            # diameters data have been stored in by fillSplittingInfo
            medDiam = leaf.private.diameters[1]
            #
            push!(leafDiameters, medDiam)
            push!(leaves, leaf)
        end
        res_iter = iterate(rptree.eventDict, res_iter[2])
    end
    meanDiam = mean(leafDiameters)
    @printf stdout "\n meanDiameter at max depth = %f" meanDiam
    @printf stdout "\n nb split diam = %d" nbsplitDiam
    @printf stdout "\n nb split proj = %d\n" nbsplitProj
    #
    return leaves, leafDiameters
end    # end of analyzeSplittingInfo

"""
    `function fillItemLeafDict(rptree::RPTree)`

    A function to fill a dictionary giving the leaf in which
    is initial item was dipatched
    return a `Dict{Int64,TreeNode}`
    """
function fillItemLeafDict(rptree::RPTree)
    itemLeafDict = Dict{Int64, TreeNode{KeyVector, RPTreeEvent}}()
    # must check that eventDict is filled (or fillSplittingInfo has been called)
    if length(rptree.eventDict) == 0
        fillSplittingInfo(rptree)
    end
    # loop from First left Leaf
    leaf = getFirstLeftLeaf(rptree)
    ileaf = 1
    while leaf !== nothing
        # transfer dict from leaf.data to
        idx = collect(keys(leaf.data))
        for i in idx
            itemLeafDict[i] = leaf
        end
        #
        leaf = getNextLeafRight(rptree,leaf)
    end
    return itemLeafDict
end