text:
# untyped generic code: perfect to inject our own types :)
function magic(y, func)
    l = length(y)
    k = ceil(Int, log2(l))
    for j=1:k, i=2^j:2^j:min(l, 2^k)
        y[i] = func(y[i-2^(j-1)], y[i])
    end
    for j=(k-1):-1:1, i=3*2^(j-1):2^j:min(l, 2^k)
        y[i] = func(y[i-2^(j-1)], y[i])
    end
    y
end
import Base: getindex, setindex!, length, size

mutable struct AccessArray <: AbstractArray{Nothing, 1}
    length::Int
    read::Vector
    history::Vector
end

function AccessArray(length, read = [], history = [])
    AccessArray(length, read, history)
end

length(A::AccessArray) = A.length
size(A::AccessArray) = (A.length,)

function getindex(A::AccessArray, i)
    push!(A.read, i)
    return
end

function setindex!(A::AccessArray, x, i)
    push!(A.history, (A.read, [i]))
    A.read = []
end

import Base.+
+(a::Nothing, b::Nothing)=a
A = magic(AccessArray(8), +)

using Makie, GeometryTypes
function render(A::AccessArray)
    olast = depth = 0
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        olast = maximum(y[2])
    end
    maxdepth = depth
    olast = depth = 0
    C = []
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        push!(C, ((y...,), A.length, maxdepth, depth))
        olast = maximum(y[2])
    end
    msize = 0.1
    outsize = 0.15
    x1 = Point2f0.(first.(first.(first.(C))), last.(C) .+ outsize .+ 0.05)
    x2 = Point2f0.(last.(first.(first.(C))), last.(C) .+ outsize .+ 0.05)
    x3 = Point2f0.(first.(last.(first.(C))), last.(C) .+ 1)
    connections = Point2f0[]

    yoff = Point2f0(0, msize / 2)
    ooff = Point2f0(0, outsize / 2 + 0.05)
    for i = 1:length(x3)
        push!(connections, x3[i] .- ooff, x1[i] .+ yoff, x3[i] .- ooff, x2[i] .+ yoff)
    end
    node_theme = Theme(
        markersize = msize, strokewidth = 3,
        strokecolor = :black, color = (:white, 0.0),
        axis = (
            ticks = (ranges = (1:8, 1:5),),
            names = (axisnames = ("Array Index", "Depth"),),
            frame = (axis_position = :none,)
        )
    )
    s = scatter(node_theme, x1)
    scatter!(node_theme, x2)
    scatter!(x3, color = :white, markersize = 0.2, strokewidth = 4, strokecolor = :black)
    scatter!(x3, color = :red, marker = '+', markersize = outsize)
    linesegments!(connections, color = :red)
    s
end
render(A)


text:
using VisualCircuits, Compose

width, height = 7, 4
nand1 = Gate(4.5, 1, :NAND)
nand2 = Gate(4.5, 3, :NAND)
nand3 = Gate(2.5, nand1[:in1][2], :NAND)
nand4 = Gate(2.5, nand2[:in2][2], :NAND)

circuit = compose(
    # Range of coordinate system
    context(units=UnitBox(0, 0, width, height)),

    # Gates
    render(nand1), render(nand2), render(nand3), render(nand4),

    # Connect to edges (point, x_value)
    connect(nand1[:out], width-1),
    connect(nand2[:out], width-1),
    connect(nand3[:in1], 1),
    connect(nand4[:in2], 1),
    connect((1.5, 2), 1),

    # Connect Gates
    connect(nand1[:in1], nand3[:out]),
    connect(nand2[:in2], nand4[:out]),

    # crossing lines
    connect(
        nand1[:in2],
        (3.5, nand1[:in2][2]),
        (3.5, 1.5), (5.5, 2.5),
        (5.5, nand2[:out][2]),
        nand2[:out]
    ),
    connect(
        nand2[:in1],
        (3.5, nand2[:in1][2]),
        (3.5, 2.5), (5.5, 1.5),
        (5.5, nand1[:out][2]),
        nand1[:out]
    ),
    connect(
        nand3[:in2],
        (1.5, nand3[:in2][2]),
        (1.5, nand4[:in1][2]),
        nand4[:in1]
    ),

    # connections/dots
    dots(
        (5.5, nand1[:out][2]),
        (5.5, nand2[:out][2]),
        (1.5, 2)
    ),

    # render text
    render(
        (0.5, nand3[:in1][2], "S"),
        (0.5, nand4[:in2][2], "R"),
        (0.5, 2, "C"),
        (width-0.5, nand1[:out][2], "Q"),
        (width-0.5, nand2[:out][2], "~Q")
    ),

    # compose options
    stroke("black"), linewidth(2mm)
)

# draw(PDF("Clocked_RS_FlipFlop.pdf", autosize(width, height)...), circuit)
draw(SVG("Clocked_RS_FlipFlop.svg", autosize(width, height)...), circuit)

ext: jl


text:
##################
# L2GTransformer #
##################
#
# maps points of reference geometries to global world
# and is e.g. used by FEBasisEvaluator or CellFinder
#
# needs call of update_trafo! on entry of a new cell (to update the trafo matrix and the determinant det)
#
# eval_trafo! maps local xref on cell to global x (e.g. for evaluation of data functions)
# mapderiv! gives the derivative of the mapping (for computation of derivatives of basis functions)

"""
    L2GTransformer

Transforms reference coordinates to global coordinates
"""
mutable struct L2GTransformer{Tv <: Real, Ti <: Integer, EG <: AbstractElementGeometry, CS <: AbstractCoordinateSystem}
    citem::Int
    nonlinear::Bool # so that users know if derivatives of map change in every quadrature point of cell or not
    Coords::Array{Tv,2}
    Nodes::Adjacency{Ti}
    ItemVolumes::Array{Tv,1}
    A::Matrix{Tv}
    b::Vector{Tv}
    C::Matrix{Tv} # cache for subcalculations that stay the same for each x (like adjugates)
    det::Tv
end

function L2GTransformer(EG::Type{<:AbstractElementGeometry}, grid::ExtendableGrid{Tv,Ti}, AT::Type{<:AssemblyType} = ON_CELLS)  where {Tv, Ti}
    A = zeros(Tv,size(grid[Coordinates],1),dim_element(EG))
    b = zeros(Tv,size(grid[Coordinates],1))
    return L2GTransformer{Tv,Ti,EG,grid[CoordinateSystem]}(0,false,grid[Coordinates],grid[GridComponentNodes4AssemblyType(AT)],grid[GridComponentVolumes4AssemblyType(AT)],A,b,zeros(Tv,0,0),0)
end

function L2GTransformer(EG::Union{Type{<:Tetrahedron3D}, Type{<:Parallelepiped3D}}, grid::ExtendableGrid{Tv,Ti}, AT::Type{<:AssemblyType} = ON_CELLS)  where {Tv, Ti}
    A = zeros(Tv,size(grid[Coordinates],1),dim_element(EG))
    b = zeros(Tv,size(grid[Coordinates],1))
    return L2GTransformer{Tv,Ti,EG,grid[CoordinateSystem]}(0,false,grid[Coordinates],grid[GridComponentNodes4AssemblyType(AT)],grid[GridComponentVolumes4AssemblyType(AT)],A,b,zeros(Tv,3,3),0)
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian1D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian2D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian3D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.b[3] = T.Coords[3,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[3,1] = T.Coords[3,T.Nodes[2,item]] - T.b[3]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Triangle2D,Cartesian2D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[3,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[3,item]] - T.b[2]
        T.det = 2*T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Parallelogram2D,Cartesian2D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[4,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[4,item]] - T.b[2]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Triangle2D,Cartesian3D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.b[3] = T.Coords[3,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[3,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[3,item]] - T.b[2]
        T.A[3,1] = T.Coords[3,T.Nodes[2,item]] - T.b[3]
        T.A[3,2] = T.Coords[3,T.Nodes[3,item]] - T.b[3]
        T.det = 2*T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Parallelogram2D,Cartesian3D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.b[3] = T.Coords[3,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[4,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[4,item]] - T.b[2]
        T.A[3,1] = T.Coords[3,T.Nodes[2,item]] - T.b[3]
        T.A[3,2] = T.Coords[3,T.Nodes[4,item]] - T.b[3]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Tetrahedron3D,Cartesian3D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.b[3] = T.Coords[3,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[3,item]] - T.b[1]
        T.A[1,3] = T.Coords[1,T.Nodes[4,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[3,item]] - T.b[2]
        T.A[2,3] = T.Coords[2,T.Nodes[4,item]] - T.b[2]
        T.A[3,1] = T.Coords[3,T.Nodes[2,item]] - T.b[3]
        T.A[3,2] = T.Coords[3,T.Nodes[3,item]] - T.b[3]
        T.A[3,3] = T.Coords[3,T.Nodes[4,item]] - T.b[3]
        T.det = 6*T.ItemVolumes[item]
    end
    return nothing
end

function update_trafo!(T::L2GTransformer{<:Real,<:Integer,<:Parallelepiped3D,Cartesian3D}, item::Int)
    if T.citem != item
        T.citem = item
        T.b[1] = T.Coords[1,T.Nodes[1,item]]
        T.b[2] = T.Coords[2,T.Nodes[1,item]]
        T.b[3] = T.Coords[3,T.Nodes[1,item]]
        T.A[1,1] = T.Coords[1,T.Nodes[2,item]] - T.b[1]
        T.A[1,2] = T.Coords[1,T.Nodes[4,item]] - T.b[1]
        T.A[1,3] = T.Coords[1,T.Nodes[5,item]] - T.b[1]
        T.A[2,1] = T.Coords[2,T.Nodes[2,item]] - T.b[2]
        T.A[2,2] = T.Coords[2,T.Nodes[4,item]] - T.b[2]
        T.A[2,3] = T.Coords[2,T.Nodes[5,item]] - T.b[2]
        T.A[3,1] = T.Coords[3,T.Nodes[2,item]] - T.b[3]
        T.A[3,2] = T.Coords[3,T.Nodes[4,item]] - T.b[3]
        T.A[3,3] = T.Coords[3,T.Nodes[5,item]] - T.b[3]
        T.det = T.ItemVolumes[item]
    end
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Union{Triangle2D, Parallelogram2D},Cartesian2D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.A[1,2]*xref[2] + T.b[1]
    x[2] = T.A[2,1]*xref[1] + T.A[2,2]*xref[2] + T.b[2]
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Union{Triangle2D, Parallelogram2D},Cartesian3D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.A[1,2]*xref[2] + T.b[1]
    x[2] = T.A[2,1]*xref[1] + T.A[2,2]*xref[2] + T.b[2]
    x[3] = T.A[3,1]*xref[1] + T.A[3,2]*xref[2] + T.b[3]
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Union{Tetrahedron3D, Parallelepiped3D},Cartesian3D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.A[1,2]*xref[2] + T.A[1,3]*xref[3] + T.b[1]
    x[2] = T.A[2,1]*xref[1] + T.A[2,2]*xref[2] + T.A[2,3]*xref[3] + T.b[2]
    x[3] = T.A[3,1]*xref[1] + T.A[3,2]*xref[2] + T.A[3,3]*xref[3] + T.b[3]
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian1D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.b[1]
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian2D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.b[1]
    x[2] = T.A[2,1]*xref[1] + T.b[2]
    return nothing
end

function eval_trafo!(x::Vector, T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian3D}, xref)
    x[1] = T.A[1,1]*xref[1] + T.b[1]
    x[2] = T.A[2,1]*xref[1] + T.b[2]
    x[3] = T.A[3,1]*xref[1] + T.b[3]
    return nothing
end

# EDGE1D/CARTESIAN1D map derivative
# x = a*xref + b
# Dxref/dx = a^{-1} = |E|^{-1}
function mapderiv!(M::Matrix, T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian1D}, xref)
    # transposed inverse of A
    M[1,1] = 1.0/T.det
    return nothing
end

# EDGE1D/CARTESIAN2D (tangential) map derivative
# x = A*xref + b
# Dxref/dx = A*tangent^{-1} = |E|^{-1}
function mapderiv!(M::Matrix, T::L2GTransformer{<:Real,<:Integer,<:Edge1D,Cartesian2D}, xref)
    # transposed inverse of A
    M[1,1] = 1.0/T.det
    return nothing
end
# TRIANGLE2D/CARTESIAN2D map derivative
# x = A*xref + b
# Dxref/dx = A^{-T}
function mapderiv!(M::AbstractMatrix, T::L2GTransformer{<:Real,<:Integer,<:Triangle2D,Cartesian2D}, xref)
    # transposed inverse of A
    M[2,2] = T.A[1,1]/T.det
    M[2,1] = -T.A[1,2]/T.det
    M[1,2] = -T.A[2,1]/T.det
    M[1,1] = T.A[2,2]/T.det
    return nothing
end
# similar for parallelogram
function mapderiv!(M::Matrix, T::L2GTransformer{<:Real,<:Integer,<:Parallelogram2D,Cartesian2D}, xref)
    # transposed inverse of A
    M[2,2] = T.A[1,1]/T.det
    M[2,1] = -T.A[1,2]/T.det
    M[1,2] = -T.A[2,1]/T.det
    M[1,1] = T.A[2,2]/T.det
    return nothing
end

function mapderiv!(M::Matrix, T::L2GTransformer{<:Real,<:Integer,<:Tetrahedron3D,Cartesian3D}, xref)
    # adjugate matrix = T.determinant of subblocks (for faster map_deriv!)
    T.C[1,1] =   T.A[2,2]*T.A[3,3] - T.A[2,3] * T.A[3,2]
    T.C[1,2] = -(T.A[2,1]*T.A[3,3] - T.A[2,3] * T.A[3,1])
    T.C[1,3] =   T.A[2,1]*T.A[3,2] - T.A[2,2] * T.A[3,1]
    T.C[2,1] = -(T.A[1,2]*T.A[3,3] - T.A[1,3] * T.A[3,2])
    T.C[2,2] =   T.A[1,1]*T.A[3,3] - T.A[1,3] * T.A[3,1]
    T.C[2,3] = -(T.A[1,1]*T.A[3,2] - T.A[1,2] * T.A[3,1])
    T.C[3,1] =   T.A[1,2]*T.A[2,3] - T.A[1,3] * T.A[2,2]
    T.C[3,2] = -(T.A[1,1]*T.A[2,3] - T.A[1,3] * T.A[2,1])
    T.C[3,3] =   T.A[1,1]*T.A[2,2] - T.A[1,2] * T.A[2,1]

    # transposed inverse of A
    for j = 1 : 3, k = 1 : 3
        M[j,k] = T.C[j,k] / T.det
    end
    return nothing
end

function mapderiv!(M::Matrix, T::L2GTransformer{<:Real,<:Integer,<:Parallelepiped3D,Cartesian3D}, xref)

    # adjugate matrix = T.determinant of subblocks (for faster map_deriv!)
    T.C[1,1] =   T.A[2,2]*T.A[3,3] - T.A[2,3] * T.A[3,2]
    T.C[1,2] = -(T.A[2,1]*T.A[3,3] - T.A[2,3] * T.A[3,1])
    T.C[1,3] =   T.A[2,1]*T.A[3,2] - T.A[2,2] * T.A[3,1]
    T.C[2,1] = -(T.A[1,2]*T.A[3,3] - T.A[1,3] * T.A[3,2])
    T.C[2,2] =   T.A[1,1]*T.A[3,3] - T.A[1,3] * T.A[3,1]
    T.C[2,3] = -(T.A[1,1]*T.A[3,2] - T.A[1,2] * T.A[3,1])
    T.C[3,1] =   T.A[1,2]*T.A[2,3] - T.A[1,3] * T.A[2,2]
    T.C[3,2] = -(T.A[1,1]*T.A[2,3] - T.A[1,3] * T.A[2,1])
    T.C[3,3] =   T.A[1,1]*T.A[2,2] - T.A[1,2] * T.A[2,1]

    # transposed inverse of A
    for j = 1 : 3, k = 1 : 3
        M[j,k] = T.C[j,k] / T.det
    end
    return nothing
end


text:
# @testset "fps and every" begin
    #warmup
    A = fps(10; duration = 1)
    pA = previous(A)
    C = Signal(-, A, pA)
    nupdates = count(C)
    sum_updates = foldp(+, 0, C)
    avg_delay = Signal(/, sum_updates, nupdates)
    Signals.run_till_now()
    switch = Signal(false)
    A = fpswhen(switch, 10; duration = 1)

    #fps
    for tv in [true, false]
        global A,pA,C,nupdates,sum_updates,avg_delay,switch
        Signals.async_mode(tv)
        A = fps(10; duration = 1)
        pA = previous(A)
        C = Signal(-, A, pA)
        nupdates = count(C)
        sum_updates = foldp(+, 0, C)
        avg_delay = Signal(/, sum_updates, nupdates)
        sleep(1)
        @test abs(avg_delay() - 1/10) < 0.05 # relax test to pass appveyor build test

        #every
        A = every(1/10; duration = 1.5)
        pA = previous(A)
        C = Signal(-, A, pA)
        nupdates = count(C)
        sum_updates = foldp(+, 0, C)
        avg_delay = Signal(/, sum_updates, nupdates)
        sleep(1)
        @test abs(avg_delay() - 1/10) < 0.05

        #fpswhen
        switch = Signal(false)
        A = fpswhen(switch, 10; duration = 1.2)
        C = count(A)
        sleep(1)
        @test C() <= 1

        switch(true)
        sleep(1.5)
        @test C() >= 9 # relax test to pass appveyor build test
    end
# end

@testset "throttle and debounce" begin
    for tv in [true, false]
        Signals.async_mode(tv)
        #throttle
        A = every(1/100; duration = 1)
        C = throttle(x->x, A; maxfps = 10)
        nupdates = count(C)
        sleep(1)
        @test nupdates() <= 10

        #debounce
        A = Signal(time())
        B = debounce(A; delay = 1 ) do a
            time() - a
        end

        A(time())
        sleep(1.2)
        @test B() >= 1 # relax the test for appveyor windows
    end
end

@testset "for_signal" begin
    for tv in [true, false]
        Signals.async_mode(tv)
        #for_signal
        range = Signal(1:5)
        A = Signal(2)
        B = for_signal(range, A; fps = 30) do i, a
            # println(a^i)
            i+a
        end
        C = foldp(+, 0, B)
        D = count(B)
        @test D() == 1
        @test C() == 3

        D(0); state(D, 0)
        C(0)
        A(1)
        sleep(1)
        @test D() == 5
        @test C() == A()*5+6*5/2
    end
end

@testset "buffer" begin
    #warmup
    A = fps(30; duration = 1)
    B = buffer(A; buf_size = Inf, timespan = 0.5, type_stable = true)
    sleep(1)
    typeof(B()) == Vector{typeof(A())}
    abs(B()[end] - B()[1] - 0.5) < 0.15
    Signals.run_till_now()

    for tv in [true, false]
        Signals.async_mode(tv)
        A = Signal(2; strict_push = true)
        B = buffer(A; buf_size = 3, timespan = Inf)
        @test B() == [2]
        A(1.2)
        B() == [2, 1.2]
        @test typeof(B()) == Vector{Any}
        A(10); A(11); yield()
        @test B() == [2, 1.2, 10]
        A(12); A(13); yield()
        @test B() == [11, 12, 13]

        A = fps(30; duration = 1)
        B = buffer(A; buf_size = Inf, timespan = 0.5, type_stable = true)
        sleep(2)
        @test typeof(B()) == Vector{typeof(A())}
        @test abs(B()[end] - B()[1] - 0.5) < 0.15
    end
end


text:
## recurrent networks
## LSTM implementation is taken from:
## https://github.com/FluxML/Flux.jl/blob/7104fd933200833d3d73e6dec2f9572e3436f4dc/src/layers/recurrent.jl

################################################################################
#                                 Vanilla RNN                                  #
################################################################################

mutable struct RNNCell
  W_ih::AbstractMatrix
  W_hh::AbstractMatrix
  b_ih::AbstractVector
  b_hh::AbstractVector
end

function RNNCell(inp::Integer, hid::Integer)
    k_sqrt = sqrt(1 / hid)
    d = Uniform(-k_sqrt, k_sqrt)
    return RNNCell(rand(d, hid, inp), rand(d, hid, hid), rand(d, hid), rand(hid))
end

Base.show(io::IO, m::RNNCell) = print(io, "RNNCell($(size(m.W_ih, 2)) => $(size(m.W_ih, 1)))")

# using proxy function to be able to compute numeric gradient in tests
function rnn_forward(W_ih, W_hh, b_ih, b_hh, x, h)
    inp_v = W_ih * x .+ b_ih
    hid_v = W_hh * h .+ b_hh
    h_ = tanh.(inp_v .+ hid_v)
    return h_
end

# input should be of size (inp_size, batch)
function forward(m::RNNCell, x::AbstractMatrix, h::AbstractMatrix)
    rnn_forward(m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h)
end

(m::RNNCell)(x::AbstractMatrix, h::AbstractMatrix) = forward(m, x, h)

mutable struct RNN
    cell::RNNCell
end
RNN(inp_size, hid_size) = RNN(RNNCell(inp_size, hid_size))
RNN(inp_hid::Pair{Int}) = RNN(RNNCell(inp_hid[1], inp_hid[2]))

function init_hidden(m::RNN, batch_size::Integer)
    hid_size = length(m.cell.b_hh)
    return zeros(hid_size, batch_size)
end

# x_seq should be of size (inp_size, batch, seq_len)
# h should be of size (hid_size, batch)
function forward(m::RNN, x_seq::AbstractArray{T, 3}, h::AbstractArray{T, 2}) where T
    # device = Yota.guess_device([h])
    inp_size, batch, seq_len = size(x_seq)
    hid_size = length(m.cell.b_hh)
    h_all_sz = (size(h)..., 1)
    cell = m.cell
    # 1st element
    x = x_seq[:, :, 1]
    h = forward(cell, x, h)
    h_all = reshape(h, h_all_sz)
    # all other elements
    for i=2:size(x_seq, 3)
        x = x_seq[:, :, i]
        h = forward(cell, x, h)
        h_all = cat(h_all, reshape(h, h_all_sz); dims=3)
    end
    return h_all, h
end

(m::RNN)(x::AbstractArray{T,3}, h::AbstractArray{T,2}) where T = forward(m, x, h)

################################################################################
#                                  LSTM                                        #
################################################################################

mutable struct LSTMCell
  W_ih::AbstractMatrix
  W_hh::AbstractMatrix
  b_ih::AbstractVector
  b_hh::AbstractVector
end

function LSTMCell(inp::Integer, hid::Integer)
    k_sqrt = sqrt(1 / hid)
    d = Uniform(-k_sqrt, k_sqrt)
    return LSTMCell(rand(d, 4*hid, inp), rand(d, 4*hid, hid), rand(d, 4*hid), rand(4*hid))
end
LSTMCell(inp_hid::Pair{Int}) = LSTMCell(inp_hid[1], inp_hid[2])

Base.show(io::IO, m::LSTMCell) = print(io, "LSTMCell($(size(m.W_ih, 2)) => $(size(m.W_ih, 1)))")

slice(len, idx) = (1:len) .+ len*(idx-1)

function lstm_forward(W_ih, W_hh, b_ih, b_hh, x, h, c)
    hid_len = size(h, 1)
    σ = tanh
    # weight breakdown in PyTorch is described here:
    # https://github.com/pytorch/pytorch/blob/0c48092b2270d56cdab327bd1ff0ca89f5b7d569/torch/nn/modules/rnn.py#L479-L487
    # however, we don't split weights into W_ii, W_if, etc., but instead first calculate
    # all linear transformations as `y` and only after that take slices of this variable
    # to calculate i, f, g and o
    y = W_ih*x .+ b_ih .+ W_hh*h .+ b_hh
    i = σ.(y[slice(hid_len, 1), :])
    f = σ.(y[slice(hid_len, 2), :])
    g = tanh.(y[slice(hid_len, 3), :])
    o = σ.(y[slice(hid_len, 4), :])
    c_ = f .* c .+ i .* g
    h_ = o .* tanh.(c_)
    return h_, c_
end

# x should be of size (inp_size, batch_size)
# h and c should be of size (hid_size, batch_size)
function forward(m::LSTMCell, x::AbstractMatrix, h::AbstractMatrix, c::AbstractMatrix)
    return lstm_forward(m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h, c)
end

(m::LSTMCell)(x::AbstractMatrix, h::AbstractMatrix, c::AbstractMatrix) = forward(m, x, h, c)

mutable struct LSTM
    cell::LSTMCell
end
LSTM(inp_size, hid_size) = LSTM(LSTMCell(inp_size, hid_size))
LSTM(inp_hid::Pair{Int}) = LSTM(LSTMCell(inp_hid[1], inp_hid[2]))

function init_hidden(m::LSTM, batch_size::Integer)
    hid_size = length(m.cell.b_hh) ÷ 4
    return zeros(hid_size, batch_size), zeros(hid_size, batch_size)
end

# x_seq should be of size (inp_size, batch, seq_len)
# h and c should be of size (hid_size, batch)
function forward(m::LSTM, x_seq::AbstractArray{T, 3}, h::AbstractArray{T, 2}, c::AbstractArray{T, 2}) where T
    # device = Yota.guess_device([h])
    inp_size, batch, seq_len = size(x_seq)
    hid_size = size(h, 1)
    h_all_sz = (size(h)..., 1)
    cell = m.cell
    # 1st element
    x = x_seq[:, :, 1]
    h, c = forward(cell, x, h, c)
    h_all = reshape(h, h_all_sz)
    # all other elements
    for i=2:size(x_seq, 3)
        x = x_seq[:, :, i]
        h, c = forward(cell, x, h, c)
        h_all = cat(h_all, reshape(h, h_all_sz); dims=3)
    end
    return h_all, h, c
end

(m::LSTM)(x::AbstractArray{T,3}, h::AbstractArray{T,2}, c::AbstractArray{T,2}) where T = forward(m, x, h, c)

################################################################################
#                                  GRU                                         #
################################################################################

mutable struct GRUCell
  W_ih::AbstractMatrix
  W_hh::AbstractMatrix
  b_ih::AbstractVector
  b_hh::AbstractVector
end

function GRUCell(inp::Integer, hid::Integer)
    k_sqrt = sqrt(1 / hid)
    d = Uniform(-k_sqrt, k_sqrt)
    return GRUCell(rand(d, 3*hid, inp), rand(d, 3*hid, hid), rand(d, 3*hid), rand(3*hid))
end
GRUCell(inp_hid::Pair{Int}) = GRUCell(inp_hid[1], inp_hid[2])

Base.show(io::IO, m::GRUCell) = print(io, "GRUCell($(size(m.W_ih, 2)) => $(size(m.W_ih, 1)))")

function gru_forward(W_ih, W_hh, b_ih, b_hh, x, h)
    hid_len = size(h, 1)
    σ = tanh
    y_i = W_ih*x .+ b_ih
    y_h = W_hh*h .+ b_hh
    # parameters are slices of precalculated y_i and y_h
    s1 = slice(hid_len, 1); s2 = slice(hid_len, 2); s3 = slice(hid_len, 3)
    r = σ.(y_i[s1, :] .+ y_h[s1, :])
    z = σ.(y_i[s2, :] .+ y_h[s2, :])
    n = tanh.(y_i[s3, :] .+ r .* y_h[s3, :])
    h_ = (1 .- z) .* n .+ z .* h
    return h_
end

# x should be of size (inp_size, batch_size)  -- check
# h and c should be of size (hid_size, batch_size)  -- check
function forward(m::GRUCell, x::AbstractMatrix, h::AbstractMatrix)
    return gru_forward(m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h)
end

(m::GRUCell)(x::AbstractMatrix, h::AbstractMatrix) = forward(m, x, h)

mutable struct GRU
    cell::GRUCell
end
GRU(inp_size, hid_size) = GRU(GRUCell(inp_size, hid_size))
GRU(inp_hid::Pair{Int}) = GRU(GRUCell(inp_hid[1], inp_hid[2]))

function init_hidden(m::GRU, batch_size::Integer)
    hid_size = length(m.cell.b_hh) ÷ 3
    return zeros(hid_size, batch_size)
end

# x_seq should be of size (inp_size, batch, seq_len)
# h should be of size (hid_size, batch)
function forward(m::GRU, x_seq::AbstractArray{T, 3}, h::AbstractArray{T, 2}) where T
    # device = Yota.guess_device([h])
    inp_size, batch, seq_len = size(x_seq)
    hid_size = size(h, 1)
    h_all_sz = (size(h)..., 1)
    cell = m.cell
    # 1st element
    x = x_seq[:, :, 1]
    h = forward(cell, x, h)
    h_all = reshape(h, h_all_sz)
    for i=2:size(x_seq, 3)
        x = x_seq[:, :, i]
        h = forward(cell, x, h)
        h_all = cat(h_all, reshape(h, h_all_sz); dims=3)
    end
    return h_all, h
end

(m::GRU)(x::AbstractArray{T,3}, h::AbstractArray{T,2}) where T = forward(m, x, h)

## SAVING THIS CODE FOR FUTURE
## PyTorch provides optionally multilayer and biderectional RNN implementation.
## There are many ways to do the same thing, but to align the code with PyTorch
## I need to clearly understand their implementation, which is out of scope for me now.
## Thus I'm just saving this in-progress piece of code for future and stick
## with simple implementation for now.
## It's worth to note that at least multiple layers can be achieved on user level
## by stacking several RNNs

## mutable struct RNN
##     num_layers::Int
##     num_directions::Int
##     cells::RNNCell
## end

## RNN(inp::Integer, hid::Integer; num_layers::Int=1, biderectional::Bool=false) =
##     RNN(num_layers, biderectional ? 2 : 1, RNNCell(inp, hid))
## RNN(inp_hid::Pair{Int}; num_layers::Int=1, biderectional::Bool=false) =
##     RNN(inp_hid[1], inp_hid[2]; num_layers=num_layers, biderectional=biderectional)

## Base.show(io::IO, m::RNN) = print(io, "RNN($(m.cell)")

## function init_hidden(m::RNN, batch_size::Integer)
##     hid_size = length(m.cell.b_hh)
##     return zeros(hid_size, batch_size, m.num_layers * m.num_directions)
## end

## # x_seq should be of size (inp_size, batch_size, seq_len)
## # H should be of size (hid_size, batch_size, num_layers*num_directions)
## function forward(m::RNN, x_seq::AbstractArray{T, 3}, H::AbstractArray{T, 3}) where T
##     inp_size, batch_size, seq_len = size(x_seq)
##     hid_size = length(m.cell.b_hh)
##     # h_all = zeros(hid_size, batch_size, m.num_layers * m.num_directions)
##     # h_all_sz = (size(h)..., 1)
##     cell = m.cell
##     for l=1:m.num_layers
##         h = H[:, :, 2l - 1]           # l-th layer, forward direction
##         for i=1:size(x_seq, 3)
##             x = x_seq[:, :, i]
##             h = forward(cell, x, h)
##             # h_all = cat(h_all, reshape(h, h_all_sz); dims=3)
##         end
##         h = H[:, :, 2l]  # l-th layer, reverse direction
##         for i=size(x_seq, 3):1
##             x = x_seq[:, :, i]
##             h = forward(cell, x, h)
##             # h_all = cat(h_all, reshape(h, h_all_sz); dims=3)
##         end
##         # TODO: add another direction
##     end
##     # return h_all, h
##     return h
## end

## (m::RNN)(x::AbstractArray{T,3}, h::AbstractArray{T,2}) where T = forward(m, x, h)


text:
println("                 *Notices Bulge*")
println("__        ___           _    _        _   _     _      ")
println("\\ \\      / / |__   __ _| |_ ( ) ___  | |_| |__ (_) ___ ")
println(" \\ \\ /\\ / /| '_ \\ / _\\`| __|// / __| | __| '_ \\| |/ __|")
println("  \\ V  V / | | | | (_| | |_    \\__ \\ | |_| | | | |\\__ \\")
println("   \\_/\\_/  |_| |_|\\__,_|\\__|   |___/ \\___|_| |_|_|/___/")

max_stars_repo_name: doczkal/OwO

max_issues_repo_name: doczkal/OwO

max_forks_repo_name: doczkal/OwO


text:
using Compat

## FIXME: broken from bitrot
## This test used to work. Something changed in a library, or Julia, or something. Now it causes a segfault
#@test( mfactor( BigInt(2)^201-1)  == @compat Dict{BigInt,Int}(7=>1,761838257287=>1,87449423397425857942678833145441=>1,1609=>1,22111=>1,193707721=>1))
#@test( mfactor( @bigint(2^201-1))  == @compat Dict{BigInt,Int}(7=>1,761838257287=>1,87449423397425857942678833145441=>1,1609=>1,22111=>1,193707721=>1))

@test mfactor(10) ==  @compat Dict( 2 => 1, 5 => 1 )

@test mfactor(2^20 - 1) == @compat Dict(31=>1,3=>1,11=>1,5=>2,41=>1)


text:
using quantumoptics

N = 200
xmin = -10.
xmax = 10.
omega = 2.

x0 = 1.
p0 = 0.
sigma0 = 1.

T = [0.:0.2:2.;]

basis_position = quantumoptics.particle.PositionBasis(xmin, xmax, N)
basis_momentum = quantumoptics.particle.MomentumBasis(basis_position)
xsamplepoints = quantumoptics.particle.samplepoints(basis_position)
psamplepoints = quantumoptics.particle.samplepoints(basis_momentum)

Tpx = quantumoptics.particle.FFTOperator(basis_momentum, basis_position)
Txp = quantumoptics.particle.FFTOperator(basis_position, basis_momentum)

# Position Basis
#p_pos = SparseOperator(quantumoptics.particle.momentumoperator(basis_position))
p_pos = Txp*quantumoptics.particle.momentumoperator(basis_momentum)*Tpx
x_pos = quantumoptics.particle.positionoperator(basis_position)
V_pos = Operator(basis_position, diagm(omega^2*xsamplepoints.^2))
H_pos = p_pos*p_pos/2 + V_pos

psi0_pos = quantumoptics.particle.gaussianstate(basis_position, x0, p0, sigma0)
@time tout, psi_pos_t = timeevolution.schroedinger(T, psi0_pos, H_pos)
@time tout, psi_pos_t = timeevolution.schroedinger(T, psi0_pos, H_pos)

# Momentum Basis
p_mom = SparseOperator(quantumoptics.particle.momentumoperator(basis_momentum))
x_mom = SparseOperator(quantumoptics.particle.positionoperator(basis_momentum))
#x_mom = SparseOperator(Tpx*quantumoptics.particle.positionoperator(basis_position)*Txp)
V_mom = SparseOperator(Tpx*Operator(basis_position, diagm(omega^2*xsamplepoints.^2))*Txp)
H_mom = p_mom*p_mom/2 + V_mom

psi0_mom = quantumoptics.particle.gaussianstate(basis_momentum, x0, p0, sigma0)
@time tout, psi_mom_t = timeevolution.schroedinger(T, psi0_mom, H_mom)
@time tout, psi_mom_t = timeevolution.schroedinger(T, psi0_mom, H_mom)

# Mixed evolution
H_mix = LazyProduct(Txp, p_mom*p_mom/2, Tpx) + SparseOperator(V_pos)
# println(expect(H_mix, psi0_pos))
# println(expect(H_pos, psi0_pos))
# println(expect(H_mix, psi0_pos + 0.1im*H_mix*psi0_pos))
# println(expect(H_pos, psi0_pos + 0.1im*H_pos*psi0_pos))
@time tout, psi_mix_t = timeevolution.schroedinger(T, psi0_pos, H_mix)
@time tout, psi_mix_t = timeevolution.schroedinger(T, psi0_pos, H_mix)

for i=1:length(T)
    println(expect(p_mom, psi_mom_t[i]), ", ", expect(p_pos, psi_pos_t[i]))
    println(expect(x_mom, psi_mom_t[i]), ", ", expect(x_pos, psi_pos_t[i]))
    println(norm(Txp*psi_mom_t[i] - psi_pos_t[i]))
    println(norm(psi_mix_t[i] - psi_pos_t[i]))
    println("================================================")
end


text:
module  LevelledTrees

import FMMTrees

struct Data{T}
    sector::Int
    values::Vector{T}
end

struct HNode{D}
    node::FMMTrees.PointerBasedTrees.Node{D}
    height::Int
end

# FMMTrees.PointerBasedTrees.data(n::HNode) = data(n.node)

struct LevelledTree{D,P,T} <: FMMTrees.PointerBasedTrees.APBTree
    nodes::Vector{HNode{D}}
    root::Int
    center::P
    halfsize::T
    levels::Vector{Int}
end

FMMTrees.root(tree::LevelledTree) = tree.root
FMMTrees.data(tree::LevelledTree, node=FMMTrees.root(tree)) = tree.nodes[node].node.data

FMMTrees.PointerBasedTrees.nextsibling(tree::LevelledTree, node_idx) = tree.nodes[node_idx].node.next_sibling
FMMTrees.PointerBasedTrees.parent(tree::LevelledTree, node_idx) = tree.nodes[node_idx].node.parent
FMMTrees.PointerBasedTrees.firstchild(tree::LevelledTree, node_idx) = tree.nodes[node_idx].node.first_child
height(tree::LevelledTree, node_idx) = tree.nodes[node_idx].height

_setfirstchild!(node::HNode, child) = HNode(FMMTrees.PointerBasedTrees.Node(node.node.data, node.node.num_children, node.node.next_sibling, node.node.parent, child), node.height)
FMMTrees.PointerBasedTrees.setfirstchild!(tree::LevelledTree, node, child) = tree.nodes[node] = _setfirstchild!(tree.nodes[node], child)

_setnextsibling!(node::HNode, next) = HNode(FMMTrees.PointerBasedTrees.Node(node.node.data, node.node.num_children, next, node.node.parent, node.node.first_child), node.height)
FMMTrees.PointerBasedTrees.setnextsibling!(tree::LevelledTree, node, next) = tree.nodes[node] = _setnextsibling!(tree.nodes[node], next)

_setheight!(node::HNode, height) = HNode(node.node, height)
setheight!(tree::LevelledTree, node, height) = tree.nodes[node] = _setheight!(tree.nodes[node], height)

function FMMTrees.insert!(tree::LevelledTree, data; parent, next, prev)
    push!(tree.nodes, HNode(FMMTrees.PointerBasedTrees.Node(data, 0, next, parent, 0), 0))
    fs = FMMTrees.PointerBasedTrees.firstchild(tree, parent)
    if fs < 1 || fs == next
        FMMTrees.PointerBasedTrees.setfirstchild!(tree, parent, length(tree.nodes))
    end
    if !(prev < 1)
        FMMTrees.PointerBasedTrees.setnextsibling!(tree, prev, length(tree.nodes))
    end

    # Walk up the tree and add one to the height tracker along the path
    id = length(tree.nodes)
    h = 1
    while true
        setheight!(tree, id, max(height(tree,id),h))
        id = FMMTrees.PointerBasedTrees.parent(tree, id)
        id < 1 && break
        h += 1
    end

    return length(tree.nodes)
end

function sector_center_size(pt, ct, hs)
    hs = hs / 2
    bl = pt .> ct
    ct = ifelse.(bl, ct.+hs, ct.-hs)
    sc = sum(b ? 2^(i-1) : 0 for (i,b) in enumerate(bl))
    return sc, ct, hs
end

function center_size(sector, center, size)
    size = size/2
    bl = [2^(i-1) & sector != 0 for i in eachindex(center)]
    center = ifelse.(bl, center.+size, center.-size)
    return center, size
end

function contains(pt, ct, hs)
    maximum(abs.(pt - ct)) <= hs
end

struct Router{T,P}
    smallest_box_size::T
    target_point::P
end

const hilbert_states = [
    [1, 2, 3, 2, 4, 5, 3, 5],
    [2, 6, 0, 7, 8, 8, 0, 7],
    [0, 9,10, 9, 1, 1,11,11],
    [6, 0, 6,11, 9, 0, 9, 8],
    [11,11, 0, 7, 5, 9, 0, 7],
    [4, 4, 8, 8, 0, 6,10, 6],
    [5, 7, 5, 3, 1, 1,11,11],
    [6, 1, 6,10, 9, 4, 9,10],
    [10, 3, 1, 1,10, 3, 5, 9],
    [4, 4, 8, 8, 2, 7, 2, 3],
    [7, 2,11, 2, 7, 5, 8, 5],
    [10, 3, 2, 6,10, 3, 4, 4]]

const hilbert_positions = [
    [0,1,3,2,7,6,4,5],
    [0,7,1,6,3,4,2,5],
    [0,3,7,4,1,2,6,5],
    [2,3,1,0,5,4,6,7],
    [4,3,5,2,7,0,6,1],
    [6,5,1,2,7,4,0,3],
    [4,7,3,0,5,6,2,1],
    [6,7,5,4,1,0,2,3],
    [2,5,3,4,1,6,0,7],
    [2,1,5,6,3,0,4,7],
    [4,5,7,6,3,2,0,1],
    [6,1,7,0,5,2,4,3]]

function FMMTrees.route!(tree::LevelledTree, state, router)

    point = router.target_point
    smallest_box_size = router.smallest_box_size

    node_idx, center, size, sfc_state, depth = state
    size <= smallest_box_size && return state
    target_sector, target_center, target_size = sector_center_size(point, center, size)
    target_pos = hilbert_positions[sfc_state][target_sector+1] + 1
    target_sfc_state = hilbert_states[sfc_state][target_sector+1] + 1
    prev_child, next_child = 0, 0
    for child in FMMTrees.children(tree, node_idx)
        child_sector = FMMTrees.data(tree,child).sector
        child_pos = hilbert_positions[sfc_state][child_sector+1]+1
        target_pos < child_pos  && (next_child = child; break)
        if child_sector == target_sector
            @assert contains(point, target_center, target_size)
            return (child, target_center, target_size, target_sfc_state, depth+1)
        end
        prev_child = child
    end
    data = Data(target_sector, Int[])
    new_node_idx = FMMTrees.insert!(tree, data, next=next_child, prev=prev_child, parent=node_idx)

    if prev_child < 1
        prev_node_idx = findprevnode(tree, router, new_node_idx)
        prev_node_idx < 1 || FMMTrees.PointerBasedTrees.setnextsibling!(tree, prev_node_idx, new_node_idx)
        if prev_node_idx < 1
            depth+1 > length(tree.levels) && resize!(tree.levels,depth+1)
            tree.levels[depth+1] = new_node_idx
        end
    end

    if next_child < 1
        next_node_idx = findnextnode(tree, router, new_node_idx)
        FMMTrees.PointerBasedTrees.setnextsibling!(tree, new_node_idx, next_node_idx)
    end

    return new_node_idx, target_center, target_size, target_sfc_state, depth+1
end

function findprevnode(tree::LevelledTree, target, new_node)

    prev_node = 0
    prev_branch = 0
    node = new_node
    parent = FMMTrees.PointerBasedTrees.parent(tree, node)
    height = 1
    while true
        if parent < 1
            break
        end
        for chd in FMMTrees.children(tree, parent)
            if chd == node
                break
            end
            if LevelledTrees.height(tree,chd) >= height
                prev_branch = chd
            end
        end
        if !(prev_branch < 1)
            break
        end
        node = parent
        parent = FMMTrees.PointerBasedTrees.parent(tree, parent)
        height += 1
    end

    parent < 1 && return prev_node

    while height > 1
        for chd in FMMTrees.children(tree, prev_branch)
            if LevelledTrees.height(tree,chd) >= height-1
                prev_node = chd
            end
        end
        prev_branch = prev_node
        height -= 1
    end

    return prev_node
end

function findnextnode(tree::LevelledTree, target, new_node)

    prev_node = 0
    prev_branch = 0
    node = new_node
    parent = FMMTrees.PointerBasedTrees.parent(tree, node)
    height = 1
    while true
        if parent < 1
            break
        end
        pay_attention = false
        for chd in FMMTrees.children(tree, parent)
            if chd == node
                pay_attention = true
                continue
            end
            pay_attention || continue
            if LevelledTrees.height(tree,chd) >= height
                prev_branch = chd
                break
            end
        end
        if !(prev_branch < 1)
            break
        end
        node = parent
        parent = FMMTrees.PointerBasedTrees.parent(tree, parent)
        height += 1
    end

    parent < 1 && return prev_node

    while height > 1
        for chd in FMMTrees.children(tree, prev_branch)
            if LevelledTrees.height(tree,chd) >= height-1
                prev_node = chd
                break
            end
        end
        prev_branch = prev_node
        height -= 1
    end

    return prev_node
end

struct LevelIterator
    tree
    level::Int
end

function Base.iterate(itr::LevelIterator, node=itr.tree.levels[itr.level])
    node < 1 && return nothing
    return node, FMMTrees.PointerBasedTrees.nextsibling(itr.tree,node)
end

Base.IteratorSize(::LevelIterator) = Base.SizeUnknown()

numlevels(t::LevelledTree) = length(t.levels)

end # module LevelledTrees


text:
include("Core.jl")
include("Metropolis.jl")


text:
using LightGraphs, MetaGraphs, SumOfSquares

export grid_graph, Partition, props

struct Partition
    graph
    get_vertex
end

linearize_index(idx,rs) = idx[1] + (length(idx) > 1 ? sum((idx[i] - 1) * prod(rs[1:i-1]) for i in 2:length(idx)) : 0)

function invert_index(idx, rs)
    inv_idx = similar(rs)
    for i in length(rs):-1:2
        fac = prod(rs[1:i-1])
        n = div(idx - 1, fac)
        inv_idx[i] = n + 1
        idx -= n*fac
    end
    inv_idx[1] = idx
    return inv_idx
end

function grid_graph(x, lb, ub, n; inf_top = zeros(Int64, length(ub)), inf_floor = zeros(Int64, length(lb)))
    degen_dim = []
    nondegen_dim = []
    for k in 1:length(n)
        if n[k] == 1
            push!(degen_dim, k)
        else
            push!(nondegen_dim, k)
        end
    end
    n_eff = n .- inf_top .- inf_floor
    Δx = zeros(length(n))
    for k in nondegen_dim
        if n_eff[k] == 0
            lb[k] = (lb[k] + ub[k])/2
            ub[k] = lb[k]
            Δx[k] = 1.0
        else
            Δx[k] = (ub[k] - lb[k]) / n_eff[k]
        end
    end
    nv = prod(n)
    mg = MetaGraph(SimpleGraph(nv))
    for i in 1:nv
        idx = invert_index(i, n)
        subset = FullSpace()
        for k in nondegen_dim
            if inf_top[k] == 1 && idx[k] == n[k]
                subset = intersect(subset, @set(x[k] >= ub[k]))
            elseif inf_floor[k] == 1 && idx[k] == 1
                subset = intersect(subset, @set(x[k] <= lb[k]))
            else
                subset = intersect(subset, @set(x[k] >= lb[k] + (idx[k]-inf_floor[k]-1)*Δx[k] && x[k] <= lb[k] + (idx[k]-inf_floor[k])*Δx[k]))
            end
        end
        for k in degen_dim
            if inf_top[k] == 0
                subset = intersect(subset, @set(x[k] <= ub[k]))
            end
            if inf_floor[k] == 0
                subset = intersect(subset, @set(x[k] >= lb[k]))
            end
        end
        set_prop!(mg, i, :cell, subset)
        for k in nondegen_dim
            if idx[k] < n[k]
                idx[k] += 1
                j = linearize_index(idx, n)
                add_edge!(mg, i, j)
                set_prop!(mg, Edge(i,j), :interface, (k, lb[k] + (idx[k]-1-inf_floor[k])*Δx[k]))
                idx[k] -= 1
            end
        end
    end
    get_vertex = function (x)
                    idx = ones(Int64,length(n))
                    for k in nondegen_dim
                        if inf_floor[k] == 1
                            idx[k] = x[k] <= lb[k] ? 1 : min(ceil(Int64, (x[k] - lb[k])/Δx[k]) + 1, n[k])
                        else
                            idx[k] = min(floor(Int64, (x[k] - lb[k])/Δx[k]) + 1, n[k])
                        end
                    end
                    return linearize_index(idx, n)
                 end
    return mg, get_vertex
end

function grid_graph(x, x_ranges)
    n = length.(x_ranges) .- 1
    nv = prod(n)

    mg = MetaGraph(SimpleGraph(nv))
    for i in 1:nv
        idx = invert_index(i, n)
        subset = FullSpace()
        for k in 1:length(idx)
            if !isinf(x_ranges[k][idx[k]])
                subset = intersect(subset, @set(x[k] >= x_ranges[k][idx[k]]))
            end
            if !isinf(x_ranges[k][idx[k]+1])
                subset = intersect(subset, @set(x[k] <= x_ranges[k][idx[k]+1]))
            end
        end
        set_prop!(mg, i, :cell, subset)
        for k in 1:length(idx)
            if idx[k] < n[k]
                idx[k] += 1
                j = linearize_index(idx, n)
                add_edge!(mg, i, j)
                set_prop!(mg, Edge(i,j), :interface, (k, x_ranges[k][idx[k]]))
                idx[k] -= 1
            end
        end
    end
    get_vertex = function (x)
                    idx = zeros(Int64, length(n))
                    for k in 1:length(n)
                        j = findfirst(m -> m >= x[k], x_ranges[k]) - 1
                        if isnothing(j)
                            @error "state outside domain"
                        end
                        idx[k] = max(j, 1)
                    end
                    return linearize_index(idx, n)
                 end
    return mg, get_vertex
end

function control_gmp(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, p::Partition, solver; moments = Dict())
    @assert typeof(CP.Objective) == LagrangeMayer

    MP = CP.MP
    nₜ = length(trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:nₜ]...]

    model = SOSModel(solver)

    @variable(model, w[vertices(p.graph), k in 1:nₜ], Poly(monomials([MP.x..., CP.t], 0:order)))

    for v in vertices(p.graph), k in 1:nₜ
        X = intersect(props(p.graph, v)[:cell], @set(CP.t >= 0 && 1-CP.t >= 0), CP.U)
        @constraint(model, extended_inf_generator(MP, w[v,k], CP.t; scale = Δt[k]) + Δt[k]*CP.Objective.l >= 0, domain = X)
    end

    for v in vertices(p.graph), k in 2:nₜ
        @constraint(model, subs(w[v,k], CP.t => 0) - subs(w[v,k-1], CP.t => 1) >= 0, domain=props(p.graph, v)[:cell])
    end

    for e in edges(p.graph), k in 1:nₜ
        i, val = props(p.graph, e)[:interface]
        if all(subs(MP.σ, MP.x[i] => val) .== 0)
            X = FullSpace()
            for p in props(p.graph, e.dst)[:cell].p
                if !(variables(p) == [MP.x[i]])
                    p = subs(p, MP.x[i] => val)
                    X = intersect(X, @set(p >= 0))
                end
            end
            X = intersect(X, CP.U, @set(CP.t >= 0 && 1-CP.t >= 0))
            @constraint(model, subs((w[e.dst,k] - w[e.src,k])*MP.f[i], MP.x[i] => val) >= 0, domain=X)
        else
            @constraint(model, subs(w[e.src,k] - w[e.dst,k], MP.x[i] => val) == 0)
        end
    end
    for v in vertices(p.graph)
        @constraint(model, CP.Objective.m - subs(w[v,nₜ], CP.t => 1) >= 0, domain = props(p.graph, v)[:cell])
    end
    @objective(model, Max, (isempty(moments) ? set_objective(w, x0, p) : set_objective(w, CP.t, moments, p)))
    return model
end

function optimal_control(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, p::Partition, solver; moments=Dict(), value_function = false)
    gmp = control_gmp(CP, x0, order, trange, p, solver, moments=moments)
    optimize!(gmp)
    lb, stat, time = extract_primal_solution(gmp)
    if value_function
        return Bounds(:control_problem, order, lb, InfoData(stat, time)), value_function_approximation(gmp, p, CP.t, trange)
    else
        return Bounds(:control_problem, order, lb, InfoData(stat, time))
    end
end

function value_function_approximation(gmp, p::Partition, t, trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:length(trange)]...]
    V_pieces = Dict((k, i) => subs(value(gmp[:w][k,i]), t => (t - (i > 1 ? trange[i] : 0))/Δt[i]) for k in vertices(p.graph), i in 1:length(trange))
    V = function (x,t)
            if t > trange[end] || t < 0
                @warn string("t = ", t, " outside the domain of V")
                i = t < 0 ? 0 : length(trange) - 1
            else
                i = findlast(ti -> ti <= t, trange[1:end-1])
            end
            k = p.get_vertex(x)
            return isnothing(i) ? (V_pieces[k,1], V_pieces[k,1](x...,t)) : (V_pieces[k,i+1], V_pieces[k,i+1](x...,t))
        end
    return V
end

function inf_horizon_control_gmp(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, p::Partition, solver; moments=Dict())
    @assert typeof(CP.Objective) == LagrangeMayer
    @assert CP.Objective.m == 0

    MP = CP.MP
    nₜ = length(trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:nₜ]...]
    ρ = CP.discount_factor
    model = SOSModel(solver)

    @variable(model, w[vertices(p.graph), k in 1:nₜ], Poly(monomials([MP.x..., CP.t], 0:order)))
    @variable(model, w∞[vertices(p.graph)], Poly(monomials(MP.x, 0:order)))

    for v in vertices(p.graph), k in 1:nₜ
        X = intersect(props(p.graph, v)[:cell], @set(CP.t >= 0 && 1-CP.t >= 0), CP.U)
        @constraint(model, extended_inf_generator(MP, w[v,k], CP.t; scale = Δt[k]) - ρ*w[v,k] + Δt[k]*CP.Objective.l >= 0, domain = X)
    end

    for v in vertices(p.graph), k in 2:nₜ
        @constraint(model, subs(w[v,k], CP.t => 0) - subs(w[v,k-1], CP.t => 1) >= 0, domain=props(p.graph, v)[:cell])
    end

    for e in edges(p.graph), k in 1:nₜ
        i, val = props(p.graph, e)[:interface]
        if all(subs(MP.σ, MP.x[i] => val) .== 0)
            X = FullSpace()
            for p in props(p.graph, e.dst)[:cell].p
                if !(variables(p) == [MP.x[i]])
                    p = subs(p, MP.x[i] => val)
                    X = intersect(X, @set(p >= 0))
                end
            end
            X = intersect(X, CP.U, @set(CP.t >= 0 && 1-CP.t >= 0))
            @constraint(model, subs((w[e.dst,k]-w[e.src,k])*MP.f[i], MP.x[i] => val) >= 0, domain=X)
        else
            @constraint(model, subs(w[e.src,k]-w[e.dst,k], MP.x[i] => val) == 0)
        end
    end

    for v in vertices(p.graph)
        X = intersect(props(p.graph, v)[:cell], @set(CP.t >= 0), CP.U)
        @constraint(model, extended_inf_generator(MP, w∞[v], CP.t) - ρ*w∞[v] + CP.Objective.l >= 0, domain = X)
        @constraint(model, w∞[v] - subs(w[v,nₜ], CP.t => 1) >= 0, domain=props(p.graph, v)[:cell])
    end

    for e in edges(p.graph), k in 1:nₜ
        i, val = props(p.graph, e)[:interface]
        if all(subs(MP.σ, MP.x[i] => val) .== 0)
            X = FullSpace()
            for p in props(p.graph, e.dst)[:cell].p
                if !(variables(p) == [MP.x[i]])
                    p = subs(p, MP.x[i] => val)
                    X = intersect(X, @set(p >= 0))
                end
            end
            X = intersect(X, CP.U, @set(CP.t >= 0))
            @constraint(model, subs((w∞[e.dst]-w∞[e.src])*MP.f[i], MP.x[i] => val) >= 0, domain=X)
        else
            @constraint(model, subs(w∞[e.src]-w∞[e.dst], MP.x[i] => val) == 0)
        end
    end
    @objective(model, Max, isempty(moments) ? set_objective(w, x0, p) : set_objective(w, CP.t, moments, p))
    return model
end

function inf_horizon_control(CP::ControlProcess, x0::AbstractVector, order::Int, trange::AbstractVector, p::Partition, solver; moments=Dict(), value_function = false)
    gmp = inf_horizon_control_gmp(CP, x0, order, trange, p, solver; moments=moments)
    optimize!(gmp)
    lb, stat, time = extract_primal_solution(gmp)
    if value_function
        return Bounds(:control_problem, order, lb, InfoData(stat, time)), value_function_approximation_inf_horizon(gmp, p, CP.t, trange)
    else
        return Bounds(:control_problem, order, lb, InfoData(stat, time))
    end
end

function value_function_approximation_inf_horizon(gmp, p::Partition, t, trange)
    nt = length(trange)
    Δt = [trange[1], [trange[i] - trange[i-1] for i in 2:nt]...]
    V_pieces = Dict((k, i) => subs(value(gmp[:w][k,i]), t => (t - (i > 1 ? trange[i] : 0))/Δt[i]) for k in vertices(p.graph), i in 1:length(trange))
    merge!(V_pieces, Dict((k, nt+1) => subs(value(gmp[:w∞][k]), t => t - trange[end])+0.0*t for k in vertices(p.graph)))
    V = function (x,t)
            k = p.get_vertex(x)
            if t < 0
                @warn string("t = ", t, " outside the domain of V")
                i = t < 0 ? 0 : length(trange) - 1
            else
                i = findlast(ti -> ti <= t, trange)
            end
            return isnothing(i) ? (V_pieces[k,1], V_pieces[k,1](x...,t)) : (V_pieces[k,i+1], V_pieces[k,i+1](x...,t))
        end
    return V
end

function set_objective(w, t, moments::Dict, p)
    obj = AffExpr(0.0)
    for v in vertices(p.graph)
        aux = subs(w[v,1], t => 0)
        for i in 1:length(aux.a)
            obj += aux.a[i]*moments[v, aux.x[i]]
        end
    end
    return obj
end

function set_objective(w, x0::Vector, p)
    return w[p.get_vertex(x0),1](x0..., 0)
end


text:
# Profile data processing. Skip plotting for now.

using RosDataProcess
using ProfileView
Profile.init(100000000, 0.001)

function profile_fun()
  println("Entering profiling code")

  bag_directory = "$(homedir())/bagfiles/decentralized_exploration_journal/"
  bags = load_directory(bag_directory)

  # show topics
  println("Topic names:")
  foreach(println, get_topic_names(bags[1]))

  params = Dict(:num_robots => "num_robots",
                :num_decentralized =>
                "kinematic_exploration/num_decentralized_planning_rounds")

  ranges = Dict(key => get_range(value, bags) for (key, value) in params)
  for num_robots = ranges[:num_robots][1:1]
    for (ii, num_decentralized) = enumerate(ranges[:num_decentralized][1:1])
      fixed = Dict(params[:num_robots] => num_robots,
                   params[:num_decentralized] => num_decentralized)
      trials = filter_bags(fixed, bags)

      entropy = read_series("/kinematic_exploration/iteration",
                            "/kinematic_exploration/entropy_reduction", trials;
                            access_data=x->x.data,
                            access_time=x-> num_robots * x.data,
                            intersect=true)
    end
  end

  foreach(close, bags)

  nothing
end

profile_fun()

@profile foreach(1:3) do _
  @time profile_fun()
end

ProfileView.view()

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


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

DOCUMENTER_DEBUG=true

using Documenter, RandomProjectionTree

makedocs(
    format = Documenter.HTML(prettyurls = false),
    sitename = "RandomProjectionTree",
    pages = Any[
        "Introduction" => "INTRO.md",
        "RPTree.jl " => "index.md",
    ]
)


text:
# julia --handle-signals=no repo_test.jl

module adapters_jdbc_repo_test

using Test # @test
using Octo.Adapters.JDBC # Repo Schema Raw
import JDBC # pathof(JDBC) JDBC.usedriver JDBC.init JDBC.destroy

Repo.debug_sql()

jdbc_test_dir = normpath(dirname(pathof(JDBC)), "..", "test")
derby_jar = normpath(jdbc_test_dir, "derby.jar")
toursdb_jar = normpath(jdbc_test_dir, "toursdb.jar")

JDBC.usedriver(derby_jar)
JDBC.init()

Repo.connect(
    adapter    = Octo.Adapters.JDBC,
    database   = Octo.DBMS.SQL,
    connection = (url="jdbc:derby:jar:($toursdb_jar)toursdb",),
)

struct Airline
end
Schema.model(Airline, table_name="airlines")

df = Repo.query(Airline)
@test size(df) == (2,)
@test df[1].AIRLINE == "AA"

flights = Repo.query([SELECT * FROM :flights])
@test size(flights) == (542,)
@test flights[1].FLIGHT_ID == "AA1111"

Repo.disconnect()

if isdir("tmptest")
    rm("tmptest", recursive=true)
end
@assert !isdir("tmptest")

Repo.connect(
    adapter    = Octo.Adapters.JDBC,
    database   = Octo.DBMS.SQL,
    connection = (url="jdbc:derby:tmptest", create=true),
)

struct FirstTable
end
Schema.model(FirstTable, table_name="firsttable")

Repo.execute(Raw("CREATE TABLE FIRSTTABLE (ID INT PRIMARY KEY, NAME VARCHAR(12))"))
multiple_changes = [
    (ID=10, NAME="TEN"),
    (ID=20, NAME="TWENTY"),
]
Repo.insert!(FirstTable, multiple_changes)
df = Repo.query(FirstTable)
@test size(df) == (2,)

Repo.disconnect()
rm("tmptest", recursive=true)

JDBC.destroy()

end # module adapters_jdbc_repo_test


text:
isdefined(Base, :__precompile__) && __precompile__()
module UnicodePlots
using Compat

export Canvas, BrailleCanvas, BarplotCanvas, DensityCanvas
export setPixel!, setPoint!, drawLine!, addRow!
export printRow, nrows, ncols

export Plot
export title, xlabel, ylabel
export annotate!, title!, xlabel!, ylabel!

export barplot, lineplot, scatterplot, stairs, histogram, densityplot
export barplot!, lineplot!, scatterplot!, stairs!, densityplot!

export spy

include("common.jl")
include("canvas.jl")
include("braillecanvas.jl")
include("densitycanvas.jl")
include("plot.jl")
include("barplot.jl")
include("histogram.jl")
include("scatterplot.jl")
include("densityplot.jl")
include("spy.jl")

end


text:
import Base.ht_keyindex

# Cf. https://github.com/JuliaStats/StatsBase.jl/issues/135
immutable UnitWeights <: AbstractVector{Int}
end
Base.getindex(w::UnitWeights, ::Integer...) = 1
Base.getindex(w::UnitWeights, ::AbstractVector) = w

# @pure only exists in Julia 0.5
if isdefined(Base, Symbol("@pure"))
    import Base.@pure
else
    macro pure(x) esc(x) end
end

# About the type inference limitation which prompts this workaround, see
# https://github.com/JuliaLang/julia/issues/10880
@pure eltypes(T) = Tuple{map(eltype, T.parameters)...}

# Internal function needed for now so that n is inferred
function _freqtable{n,T<:Real}(x::NTuple{n},
                               weights::AbstractVector{T} = UnitWeights(),
                               subset::Union{Void, AbstractVector{Int}, AbstractVector{Bool}} = nothing)
    if !isa(subset, Void)
        x = map(y -> y[subset], x)
        weights = weights[subset]
    end

    l = map(length, x)
    vtypes = eltypes(typeof(x))

    for i in 1:n
        if l[1] != l[i]
            error("arguments are not of the same length: $l")
        end
    end

    if !isa(weights, UnitWeights) && length(weights) != l[1]
        error("'weights' (length $(length(weights))) must be of the same length as vectors (length $(l[1]))")
    end

    d = Dict{vtypes, eltype(weights)}()

    for (i, el) in enumerate(zip(x...))
        index = ht_keyindex(d, el)

        if index > 0
            @inbounds d.vals[index] += weights[i]
        else
            @inbounds d[el] = weights[i]
        end
    end

    k = collect(keys(d))

    dimnames = Vector{Any}(n)
    for i in 1:n
        s = Set{vtypes.parameters[i]}()
        for j in 1:length(k)
            push!(s, k[j][i])
        end

        dimnames[i] = unique(s)
        elty = eltype(dimnames[i])
        if method_exists(isless, (elty, elty))
            sort!(dimnames[i])
        end
    end

    a = zeros(eltype(weights), map(length, dimnames)...)
    na = NamedArray(a, tuple(dimnames...), ntuple(i -> "Dim$i", n))

    for (k, v) in d
        na[k...] = v
    end

    na
end

freqtable{T<:Real}(x::AbstractVector...;
                   weights::AbstractVector{T} = UnitWeights(),
                   subset::Union{Void, AbstractVector{Int}, AbstractVector{Bool}} = nothing) =
    _freqtable(x, weights, subset)

# Internal function needed for now so that n is inferred
function _freqtable{n}(x::NTuple{n, PooledDataVector}, usena = false)
	len = map(length, x)
	lev = map(levels, x)

	for i in 1:n
	    if len[1] != len[i]
	        error(string("arguments are not of the same length: ", tuple(len...)))
	    end
	end

	if usena
        dims = map(l -> length(l) + 1, lev)
	    sizes = cumprod([dims...])
	    a = zeros(Int, dims)

	    for i in 1:len[1]
	        el = Int(x[1].refs[i])

            if el == 0
	            el = dims[1]
	        end

	        for j in 2:n
	            val = Int(x[j].refs[i])

	            if val == zero(val)
	                val = dims[j]
	            end

	            el += Int((val - 1) * sizes[j - 1])
	        end

	        a[el] += 1
	    end

	    NamedArray(a, map(l -> [l; "NA"], lev), ntuple(i -> "Dim$i", n))
	else
        dims = map(length, lev)
	    sizes = cumprod([dims...])
	    a = zeros(Int, dims)

	    for i in 1:len[1]
	        pos = (x[1].refs[i] != zero(UInt))
	        el = Int(x[1].refs[i])

	        for j in 2:n
	            val = x[j].refs[i]

	            if val == zero(val)
	                pos = false
	                break
	            end

	            el += Int((val - 1) * sizes[j - 1])
	        end

	        if pos
	            @inbounds a[el] += 1
	        end
	    end

	    NamedArray(a, lev, ntuple(i -> "Dim$i", n))
	end
end

freqtable(x::PooledDataVector...; usena::Bool = false) = _freqtable(x, usena)

function freqtable(d::DataFrame, x::Symbol...; args...)
    a = freqtable([d[y] for y in x]...; args...)
    setdimnames!(a, x)
    a
end


text:
## Counting invalidations

using SnoopCompile
umis(list) = unique(filter(x->isa(x, Core.MethodInstance), list))
# length(umis(@snoopr using <SomePkg>))

## Counting MethodInstances

using MethodAnalysis
using MethodAnalysis: equal
function getmis(parent...)  # use either 0 or 1 argument
    mis = Set{Core.MethodInstance}()
    visit(parent...) do mi
        isa(mi, Core.MethodInstance) && push!(mis, mi)
    end
    return collect(mis)
end
mis = getmis()

# counting by module

const modcounts = Dict{Module,Int}()
for mi in mis
    modcounts[mi.def.module] = get(modcounts, mi.def.module, 0) + 1
end

## When invalidations do and do not happen: the applyf example

f(::Int) = 1
f(::Bool) = 2
function applyf(container)
    x1 = f(container[1])
    x2 = f(container[2])
    return x1 + x2
end
c = Any[1, true]
applyf(c)
w2 = worlds(only(getmis(applyf)))
code2 = only(code_typed(applyf, (Vector{Any},)))
f(::String) = 3
applyf(c)
w3 = worlds(only(getmis(applyf)))
code3 = only(code_typed(applyf, (Vector{Any},)))
f(::AbstractVector) = 4   # if we replaced this with Vector{Int} we'd not have the unnecessary(?) invalidation
applyf(c)
w4 = worlds(only(getmis(applyf)))
code4 = only(code_typed(applyf, (Vector{Any},)))
f(::Missing) = 5
applyf(c)
w5 = worlds(only(getmis(applyf)))
code5 = only(code_typed(applyf, (Vector{Any},)))
f(::Nothing) = 6
applyf(c)
w6 = worlds(only(getmis(applyf)))
code6 = only(code_typed(applyf, (Vector{Any},)))

# One apparently-unnecessary invalidation?
@show equal(code3, code4)
@show w3 == w4

## Determining the source of invalidations

# Mechanism 1: MethodInstances with TypeVar parameters
function hastv(typ)
    isa(typ, UnionAll) && return true
    if isa(typ, DataType)
        for p in typ.parameters
            hastv(p) && return true
        end
    end
    return false
end
mitv = filter(mi->hastv(mi.specTypes), mis)
mitvi = Set(with_all_backedges(mitv))


text:
############################################################################################
"""
$(TYPEDEF)

Default configuration for stepsize adaption.

# Fields
$(TYPEDFIELDS)
"""
struct ConfigStepnumber{A<:BaytesCore.UpdateBool}
    "Step Number adaption"
    stepnumberadaption::A
    "Initial number of Steps"
    steps::Int64
    "Maximal number of steps"
    maxsteps::Int64
    "Desired Integration time"
    ∫dt::Float64
    function ConfigStepnumber(;
        stepnumberadaption=BaytesCore.UpdateTrue(),
        steps=10,
        maxsteps=100,
        ∫dt=1.0,
        )
        @argcheck steps <= maxsteps "Upper bound for L is maxsteps."
        @argcheck 0.0 < ∫dt "Integration time needs to be positive."
        return new{typeof(stepnumberadaption)}(
            stepnumberadaption, steps, maxsteps,∫dt
            )
    end
end

############################################################################################
"""
$(TYPEDEF)
Default Configuration for HMC sampler.

# Fields
$(TYPEDFIELDS)
"""
struct ConfigHMC{
    K<:KineticEnergy,S<:ConfigStepnumber
} <: AbstractConfiguration
    "Target Acceptance Rate"
    δ::Float64
    "Default size for tuning iterations in each cycle."
    window::ConfigTuningWindow
    "Kinetic Energy used in Hamiltonian: GaussianKineticEnergy"
    energy::K
    "Step Number adaption"
    stepnumber::S
    function ConfigHMC(
        δ::Float64,
        window::ConfigTuningWindow,
        energy::K,
        stepnumber::S,
    ) where {K<:KineticEnergy,S<:ConfigStepnumber}
        @argcheck 0.0 < δ <= 1.0 "Acceptance rate not bounded between 0 and 1."
        return new{K,S}(
            δ,
            window,
            energy,
            stepnumber
        )
    end
end

"""
$(SIGNATURES)
Initialize HMC custom configurations.

# Examples
```julia
```

"""
function init(
    ::Type{AbstractConfiguration},
    mcmc::Type{HMC},
    objective::Objective,
    proposalconfig::ConfigProposal;
    ## Target Acceptance Rate
    δ=0.80,
    ## MCMC Phase tune variables based on standard HMC best practice target acceptance rate of 0.80
    window= ConfigTuningWindow(
        [1, 5, 1],
        [75, 50, 25],
        [Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()]
        ),
    ## Kinetic energy in Hamiltonian
    energy=GaussianKineticEnergy(
        init(objective.model.info.reconstruct.default.output, proposalconfig.metric, length(objective.tagged))...,
    ),
    stepnumber = ConfigStepnumber()
)
    return ConfigHMC(
        δ,
        window,
        energy,
        stepnumber
    )
end

############################################################################################
#export
export ConfigHMC, ConfigStepnumber


text:
using Helpme
using Base.Test

for (error, fnstr, suggestion) in Helpme.testbase
	results = Helpme.search(error, fnstr)
	@test suggestion in map(x->x.suggestion, results)
end


text:
###############################################################################
#
#   fmpq_abs_series.jl : Power series over flint fmpq rationals (using fmpq_poly)
#
###############################################################################

export fmpq_abs_series, FmpqAbsSeriesRing, tan, tanh, sin, sinh, asin, asinh,
       atan, atanh, sqrt, log

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fmpq_abs_series)
   prec = length(a) - 1
   prec < 0 && throw(DomainError())
   z = parent(a)()
   z.prec = prec
   z.parent = parent(a)
   return z
end

elem_type(::Type{FmpqAbsSeriesRing}) = fmpq_abs_series

parent_type(::Type{fmpq_abs_series}) = FmpqAbsSeriesRing

base_ring(R::FmpqAbsSeriesRing) = R.base_ring

var(a::FmpqAbsSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FmpqAbsSeriesRing) = R.prec_max

function normalise(a::fmpq_abs_series, len::Int)
   if len > 0
      c = fmpq()
      ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq_abs_series}, Int), c, a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
            (Ref{fmpq}, Ref{fmpq_abs_series}, Int), c, a, len - 1)
      end
   end

   return len
end

function coeff(x::fmpq_abs_series, n::Int)
   if n < 0
      return fmpq(0)
   end
   z = fmpq()
   ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq_abs_series}, Int), z, x, n)
   return z
end

function length(x::fmpq_abs_series)
   return ccall((:fmpq_poly_length, :libflint), Int, (Ref{fmpq_abs_series},), x)
end

precision(x::fmpq_abs_series) = x.prec

zero(R::FmpqAbsSeriesRing) = R(0)

one(R::FmpqAbsSeriesRing) = R(1)

function gen(R::FmpqAbsSeriesRing)
   z = fmpq_abs_series([fmpq(0), fmpq(1)], 2, max_precision(R))
   z.parent = R
   return z
end

function deepcopy_internal(a::fmpq_abs_series, dict::ObjectIdDict)
   z = fmpq_abs_series(a)
   z.prec = a.prec
   z.parent = parent(a)
   return z
end

function isgen(a::fmpq_abs_series)
   return precision(a) == 0 || ccall((:fmpq_poly_is_x, :libflint), Bool,
                            (Ref{fmpq_abs_series},), a)
end

iszero(a::fmpq_abs_series) = length(a) == 0

isunit(a::fmpq_abs_series) = valuation(a) == 0 && isunit(coeff(a, 0))

function isone(a::fmpq_abs_series)
   return precision(a) == 0 || ccall((:fmpq_poly_is_one, :libflint), Bool,
                                (Ref{fmpq_abs_series},), a)
end

# todo: write an fmpq_poly_valuation
function valuation(a::fmpq_abs_series)
   for i = 1:length(a)
      if !iszero(coeff(a, i - 1))
         return i - 1
      end
   end
   return precision(a)
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FmpqAbsSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fmpq_abs_series}) = show_minus_one(fmpq)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fmpq_abs_series)
   z = parent(x)()
   ccall((:fmpq_poly_neg, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}),
               z, x)
   z.prec = x.prec
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpq_abs_series, b::fmpq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function -(a::fmpq_abs_series, b::fmpq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   z = parent(a)()
   z.prec = prec
   ccall((:fmpq_poly_sub_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function *(a::fmpq_abs_series, b::fmpq_abs_series)
   check_parent(a, b)
   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(a)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   z = parent(a)()
   z.prec = prec

   if lena == 0 || lenb == 0
      return z
   end

   lenz = min(lena + lenb - 1, prec)

   ccall((:fmpq_poly_mullow, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, b, lenz)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpq_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fmpq_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, y, x)
   return z
end

function *(x::fmpz, y::fmpq_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fmpq_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpz}),
               z, y, x)
   return z
end

function *(x::fmpq, y::fmpq_abs_series)
   z = parent(y)()
   z.prec = y.prec
   ccall((:fmpq_poly_scalar_mul_fmpq, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq}),
               z, y, x)
   return z
end

*(x::fmpq_abs_series, y::Int) = y*x

*(x::fmpq_abs_series, y::fmpz) = y*x

*(x::fmpq_abs_series, y::fmpq) = y*x

*(x::fmpq_abs_series, y::Integer) = x*fmpz(y)

*(x::Integer, y::fmpq_abs_series) = fmpz(x)*y

*(x::fmpq_abs_series, y::Rational) = x*fmpq(y)

*(x::Rational, y::fmpq_abs_series) = fmpq(x)*y

+(x::fmpq_abs_series, y::Rational) = x + fmpq(y)

+(x::Rational, y::fmpq_abs_series) = fmpq(x) + y

-(x::fmpq_abs_series, y::Rational) = x - fmpq(y)

-(x::Rational, y::fmpq_abs_series) = fmpq(x) - y

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpq_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   z.prec = x.prec + len
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, x, len)
   return z
end

function shift_right(x::fmpq_abs_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = length(x)
   z = parent(x)()
   if len >= xlen
      z.prec = max(0, x.prec - len)
   else
      z.prec = x.prec - len
      ccall((:fmpq_poly_shift_right, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, x, len)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fmpq_abs_series, prec::Int)
   prec < 0 && throw(DomainError())
   if x.prec <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   ccall((:fmpq_poly_set_trunc, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, x, prec)
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpq_abs_series, b::Int)
   b < 0 && throw(DomainError())
   # special case powers of x for constructing power series efficiently
   if precision(a) > 0 && isgen(a) && b > 0
      return shift_left(a, b - 1)
   elseif length(a) == 1
      z = parent(a)(coeff(a, 0)^b)
      set_prec!(z, precision(a))
      return z
   elseif b == 0
      z = one(parent(a))
      set_prec!(z, precision(a))
      return z
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit !=0
         z = z*z
         if (UInt(bit) & b) != 0
            z *= a
         end
         bit >>= 1
      end
      return z
   end
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fmpq_abs_series, y::fmpq_abs_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)

   n = max(length(x), length(y))
   n = min(n, prec)

   return Bool(ccall((:fmpq_poly_equal_trunc, :libflint), Cint,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               x, y, n))
end

function isequal(x::fmpq_abs_series, y::fmpq_abs_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || length(x) != length(y)
      return false
   end
   return Bool(ccall((:fmpq_poly_equal, :libflint), Cint,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               x, y, length(x)))
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(x::fmpq_abs_series, y::Rational{T}) where T <: Union{Int, BigInt} = x == fmpq(y)

==(x::fmpq_abs_series, y::Integer) = x == fmpz(y)

==(x::Rational{T}, y::fmpq_abs_series) where T <: Union{Int, BigInt} = y == x

==(x::Integer, y::fmpq_abs_series) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpq_abs_series, y::fmpq_abs_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   v2 = valuation(y)
   v1 = valuation(x)
   if v2 != 0
      if v1 >= v2
         x = shift_right(x, v2)
         y = shift_right(y, v2)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec, y.prec - v2 + v1)
   z = parent(x)()
   z.prec = prec
   ccall((:fmpq_poly_div_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, x, y, prec)
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpq_abs_series, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fmpq_poly_scalar_div_si, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, x, y)
   return z
end

function divexact(x::fmpq_abs_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fmpq_poly_scalar_div_fmpz, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpz}),
               z, x, y)
   return z
end

function divexact(x::fmpq_abs_series, y::fmpq)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   ccall((:fmpq_poly_scalar_div_fmpq, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq}),
               z, x, y)
   return z
end

divexact(x::fmpq_abs_series, y::Integer) = divexact(x, fmpz(y))

divexact(x::fmpq_abs_series, y::Rational{T}) where T <: Union{Int, BigInt} = divexact(x, fmpq(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpq_abs_series)
  iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ccall((:fmpq_poly_inv_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

function Base.exp(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in exp")
   if length(a) == 0 || a.prec == 1
      return parent(a)([fmpq(1)], 1, a.prec)
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_exp_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   log(a::fmpq_abs_series)
> Return log$(a)$. Requires the constant term to be one.
"""
function log(a::fmpq_abs_series)
   !isone(coeff(a, 0)) && error("Constant term not one in log")
   if length(a) == 1 || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_log_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   tan(a::fmpq_abs_series)
> Return tan$(a)$. Requires a zero constant term.
"""
function tan(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in tan")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_tan_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   tanh(a::fmpq_abs_series)
> Return tanh$(a)$. Requires a zero constant term.
"""
function tanh(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in tanh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_tanh_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   sin(a::fmpq_abs_series)
> Return sin$(a)$. Requires a zero constant term.
"""
function sin(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in sin")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_sin_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   sinh(a::fmpq_abs_series)
> Return sinh$(a)$. Requires a zero constant term.
"""
function sinh(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in sinh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_sinh_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   cos(a::fmpq_abs_series)
> Return cos$(a)$. Requires a zero constant term.
"""
function cos(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in cos")
   if length(a) == 0 || a.prec == 1
      return one(parent(a))
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_cos_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   cosh(a::fmpq_abs_series)
> Return cosh$(a)$. Requires a zero constant term.
"""
function cosh(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in cosh")
   if length(a) == 0 || a.prec == 1
      return one(parent(a))
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_cosh_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   asin(a::fmpq_abs_series)
> Return asin$(a)$. Requires a zero constant term.
"""
function asin(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in asin")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_asin_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   asinh(a::fmpq_abs_series)
> Return asinh$(a)$. Requires a zero constant term.
"""
function asinh(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in asinh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_asinh_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   atan(a::fmpq_abs_series)
> Return atan$(a)$. Requires a zero constant term.
"""
function atan(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in atan")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_atan_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   atanh(a::fmpq_abs_series)
> Return atanh$(a)$. Requires a zero constant term.
"""
function atanh(a::fmpq_abs_series)
   !iszero(coeff(a, 0)) && error("Constant term not zero in atanh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_atanh_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

doc"""
   sqrt(a::fmpq_abs_series)
> Return the power series square root of $a$. Requires a constant term equal to
> one.
"""
function sqrt(a::fmpq_abs_series)
   !isone(coeff(a, 0)) && error("Constant term not one in sqrt")
   z = parent(a)()
   z.prec = a.prec
   ccall((:fmpq_poly_sqrt_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, a.prec)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function setcoeff!(z::fmpq_abs_series, n::Int, x::fmpq)
   ccall((:fmpq_poly_set_coeff_fmpq, :libflint), Void,
                (Ref{fmpq_abs_series}, Int, Ref{fmpq}),
               z, n, x)
   return z
end

function mul!(z::fmpq_abs_series, a::fmpq_abs_series, b::fmpq_abs_series)
   lena = length(a)
   lenb = length(b)

   aval = valuation(a)
   bval = valuation(b)

   prec = min(a.prec + bval, b.prec + aval)
   prec = min(prec, max_precision(parent(a)))

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = min(lena + lenb - 1, prec)
   if lenz < 0
      lenz = 0
   end

   z.prec = prec
   ccall((:fmpq_poly_mullow, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               z, a, b, lenz)
   return z
end

function addeq!(a::fmpq_abs_series, b::fmpq_abs_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenz = max(lena, lenb)
   a.prec = prec
   ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Ref{fmpq_abs_series}, Int),
               a, a, b, lenz)
   return a
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fmpq_abs_series}, ::Type{T}) where {T <: Integer} = fmpq_abs_series

promote_rule(::Type{fmpq_abs_series}, ::Type{Rational{T}}) where T <: Union{Int, BigInt} = fmpq_abs_series

promote_rule(::Type{fmpq_abs_series}, ::Type{fmpz}) = fmpq_abs_series

promote_rule(::Type{fmpq_abs_series}, ::Type{fmpq}) = fmpq_abs_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FmpqAbsSeriesRing)()
   z = fmpq_abs_series()
   z.prec = a.prec_max
   z.parent = a
   return z
end

function (a::FmpqAbsSeriesRing)(b::Integer)
   if b == 0
      z = fmpq_abs_series()
      z.prec = a.prec_max
   else
      z = fmpq_abs_series([fmpq(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FmpqAbsSeriesRing)(b::fmpz)
   if iszero(b)
      z = fmpq_abs_series()
      z.prec = a.prec_max
   else
      z = fmpq_abs_series([fmpq(b)], 1, a.prec_max)
   end
   z.parent = a
   return z
end

function (a::FmpqAbsSeriesRing)(b::fmpq)
   if iszero(b)
      z = fmpq_abs_series()
      z.prec = a.prec_max
   else
      z = fmpq_abs_series([b], 1, a.prec_max)
   end
   z.parent = a
   return z
end

(a::FmpqAbsSeriesRing)(b::Rational{T}) where T <: Union{Int, BigInt} = a(fmpq(b))

function (a::FmpqAbsSeriesRing)(b::fmpq_abs_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FmpqAbsSeriesRing)(b::Array{fmpq, 1}, len::Int, prec::Int)
   z = fmpq_abs_series(b, len, prec)
   z.parent = a
   return z
end

###############################################################################
#
#   PowerSeriesRing constructor
#
###############################################################################

function PowerSeriesRing(R::FlintRationalField, prec::Int, s::AbstractString; model=:capped_relative, cached = true)
   S = Symbol(s)

   if model == :capped_relative
      parent_obj = FmpqRelSeriesRing(prec, S, cached)
   elseif model == :capped_absolute
      parent_obj = FmpqAbsSeriesRing(prec, S, cached)
   end

   return parent_obj, gen(parent_obj)
end


text:
using Distributions
using MonteCarloMeasurements

export sourceImportanceSampler
function sourceImportanceSampler(p,q;ℓ=:ℓ)
    procp(p, st::Follows) = :($ℓ += logpdf($(st.rhs), $(st.x)))
    procp(p, st::Let)     = convert(Expr, st)
    procp(p, st::Return)  = nothing
    procp(p, st::LineNumber) = convert(Expr, st)

    procq(q, st::Follows) = @q begin
        $(st.x) = rand($(st.rhs))
        $ℓ -= logpdf($(st.rhs), $(st.x))
    end
    procq(q, st::Let)     = convert(Expr, st)
    procq(q, st::Return)  = nothing
    procq(q, st::LineNumber) = convert(Expr, st)

    pbody = buildSource(p, procp) |> striplines
    qbody = buildSource(q, procq) |> striplines

    kwargs = freeVariables(q) ∪ arguments(p)
    kwargsExpr = Expr(:tuple,kwargs...)

    stochExpr = begin
        vals = map(stochastic(q)) do x Expr(:(=), x,x) end
        Expr(:tuple, vals...)
    end

    @gensym logimportance
    result = @q function $logimportance(pars)
        @unpack $kwargsExpr = pars
        $ℓ = 0.0
        $qbody
        $pbody
        return ($ℓ, $stochExpr)
    end

    flatten(result)
end

export sourceParticleImportance
function sourceParticleImportance(p,q;ℓ=:ℓ)
    p = canonical(p)
    q = canonical(q)
    @gensym N
    # This determines how to initialize a Particle for a given expression
    vars(expr) = (bound(p) ∪ bound(q) ∪ stochastic(p) ∪ stochastic(q)) ∩ variables(expr)

    procp(p, st::Follows) = :($ℓ += logpdf($(st.rhs), $(st.x)))
    procp(p, st::Let)     = convert(Expr, st)
    procp(p, st::Return)  = nothing
    procp(p, st::LineNumber) = convert(Expr, st)

    function procq(q, st::Follows)
        if isempty(vars(st.rhs))
            @q begin
                $(st.x) = Particles($N, $(st.rhs))
                $ℓ -= logpdf($(st.rhs), $(st.x))
            end
        else
            @q begin
                $(st.x) = rand($(st.rhs))
                $ℓ -= logpdf($(st.rhs), $(st.x))
            end
        end
    end
    procq(q, st::Let)     = convert(Expr, st)
    procq(q, st::Return)  = nothing
    procq(q, st::LineNumber) = convert(Expr, st)

    pbody = buildSource(p, procp) |> striplines
    qbody = buildSource(q, procq) |> striplines

    kwargs = freeVariables(q) ∪ arguments(p)
    kwargsExpr = Expr(:tuple,kwargs...)

    stochExpr = begin
        vals = map(stochastic(q)) do x Expr(:(=), x,x) end
        Expr(:tuple, vals...)
    end

    @gensym particlemportance
    result = @q function $particlemportance($N, pars)
        @unpack $kwargsExpr = pars
        $ℓ = 0.0 * Particles($N, Uniform())
        $qbody
        $pbody
        return ($ℓ, $stochExpr)
    end

    flatten(result)
end

#     @gensym rand

#     flatten(@q (
#         function $rand(args...;kwargs...)
#             @unpack $argsExpr = kwargs
#             # kwargs = Dict(kwargs)
#             $body
#             $stochExpr
#         end
#     ))

# end

export makeImportanceSampler
function makeImportanceSampler(p,q)
    fpre = @eval $(sourceImportanceSampler(p,q))

    function f(r;kwargs...)
        qsample = r(;kwargs...)
        ℓ = Base.invokelatest(fpre, merge(kwargs, pairs(qsample)))
        return (qsample, ℓ)
    end

    return f
end

# p = @model μ begin
#     x ~ Normal(μ, 1)
# end

# q = @model μ begin
#     x ~ Cauchy(μ)
# end

# julia> sourceImportanceLogWeights(p,q)
# :(function ##logdensity#669(pars)
#       @unpack (x, μ) = pars
#       ℓ = 0.0
#       ℓ += logpdf(Normal(μ, 1), x)
#       ℓ -= logpdf(Cauchy(μ), x)
#       return ℓ
#   end)

# impsamp =  makeImportanceSampler(p,q)

# impsamp(;μ=3.0)


text:
using SQLite
using DataFrames


text:
############################################################
# joAbstractDAparallelLinearOperator #################################
############################################################

# includes
include("joAbstractDAparallelLinearOperator/constructors.jl")
include("joAbstractDAparallelLinearOperator/base_functions.jl")
include("joAbstractDAparallelLinearOperator/extra_functions.jl")