text:
export batch_normalize,
    batch_normalize!,
    rotmat,
    rand_hermitian,
    rand_unitary,
    sprand_hermitian,
    sprand_unitary,
    general_controlled_gates,
    general_c1_gates,
    linop2dense,
    # kron
    hilbertkron,
    batched_kron!,
    batched_kron,
    # norms
    trnorm,
    nucnorm,
    # fidelity
    pure_state_fidelity,
    density_fidelity,
    purification_fidelity,
    # matrix tools
    autostatic,
    rot_mat

using LuxurySparse, LinearAlgebra, BitBasis, SparseArrays
import LinearAlgebra: svdvals

"""
nucnorm(m)

Computes the nuclear norm of a matrix `m`.
"""
function nucnorm(m::AbstractMatrix)
    norm(svdvals(m), 1)
end

"""
trnorm(m)

Computes the trace norm of a matrix `m`.
"""
trnorm(m::AbstractMatrix) = nucnorm(m)

"""
    batch_normalize!(matrix)

normalize a batch of vector.
"""
function batch_normalize!(s::AbstractMatrix, p::Real = 2)
    B = size(s, 2)
    for i = 1:B
        normalize!(view(s, :, i), p)
    end
    s
end

"""
    batch_normalize

normalize a batch of vector.
"""
function batch_normalize(s::AbstractMatrix, p::Real = 2)
    ts = copy(s)
    batch_normalize!(ts, p)
end

"""
    hilbertkron(num_bit::Int, gates::Vector{AbstractMatrix}, locs::Vector{Int}) -> AbstractMatrix

Return general kronecher product form of gates in Hilbert space of `num_bit` qubits.

* `gates` are a list of matrices.
* `start_locs` should have the same length as `gates`, specifing the gates starting positions.
"""
function hilbertkron(num_bit::Int, ops::Vector{<:AbstractMatrix}, start_locs::Vector{Int})
    sizes = [log2dim1(op) for op in ops]
    start_locs = num_bit .- start_locs .- sizes .+ 2

    order = sortperm(start_locs)
    sorted_ops = ops[order]
    sorted_start_locs = start_locs[order]
    num_ids = vcat(
        sorted_start_locs[1] - 1,
        diff(push!(sorted_start_locs, num_bit + 1)) .- sizes[order],
    )

    _wrap_identity(sorted_ops, num_ids)
end

# kron, and wrap matrices with identities.
function _wrap_identity(
    data_list::Vector{T},
    num_bit_list::Vector{Int},
) where {T<:AbstractMatrix}
    length(num_bit_list) == length(data_list) + 1 || throw(ArgumentError())
    ⊗ = kron
    reduce(
        zip(data_list, num_bit_list[2:end]);
        init = IMatrix(1 << num_bit_list[1]),
    ) do x, y
        x ⊗ y[1] ⊗ IMatrix(1 << y[2])
    end
end

batched_kron(a, b, c, xs...) =
    Base.afoldl(batched_kron, (batched_kron)((batched_kron)(a, b), c), xs...)

function batched_kron(A::AbstractArray{T,3}, B::AbstractArray{S,3}) where {T,S}
    @assert size(A, 3) == size(B, 3) "batch size mismatch"
    C = Array{Base.promote_op(*, T, S),3}(
        undef,
        size(A, 1) * size(B, 1),
        size(A, 2) * size(B, 2),
        size(A, 3),
    )
    return batched_kron!(C, A, B)
end

function batched_kron!(
    C::Array{T,3},
    A::AbstractArray{T1,3},
    B::AbstractArray{T2,3},
) where {T,T1,T2}
    @assert !Base.has_offset_axes(A, B)
    m, n = size(A)
    p, q = size(B)
    @inbounds for k = 1:size(C, 3)
        for s = 1:n, r = 1:m, w = 1:q, v = 1:p
            C[p*(r-1)+v, q*(s-1)+w, k] = A[r, s, k] * B[v, w, k]
        end
    end
    return C
end

@static if !@isdefined(kron!)
    export kron!
    # NOTE: JuliaLang/julia/pull/31069 includes this function
    function kron!(
        C::AbstractMatrix{T},
        A::AbstractMatrix{T1},
        B::AbstractMatrix{T2},
    ) where {T,T1,T2}
        @assert !Base.has_offset_axes(A, B)
        m = 1
        @inbounds for j = 1:size(A, 2), l = 1:size(B, 2), i = 1:size(A, 1)
            aij = A[i, j]
            for k = 1:size(B, 1)
                C[m] = aij * B[k, l]
                m += 1
            end
        end
        return C
    end
end

"""
    general_controlled_gates(num_bit::Int, projectors::Vector{Tp}, cbits::Vector{Int}, gates::Vector{AbstractMatrix}, locs::Vector{Int}) -> AbstractMatrix

Return general multi-controlled gates in hilbert space of `num_bit` qubits,

* `projectors` are often chosen as `P0` and `P1` for inverse-Control and Control at specific position.
* `cbits` should have the same length as `projectors`, specifing the controling positions.
* `gates` are a list of controlled single qubit gates.
* `locs` should have the same length as `gates`, specifing the gates positions.
"""
function general_controlled_gates(
    n::Int,
    projectors::Vector{<:AbstractMatrix},
    cbits::Vector{Int},
    gates::Vector{<:AbstractMatrix},
    locs::Vector{Int},
)
    IMatrix(1 << n) - hilbertkron(n, projectors, cbits) +
    hilbertkron(n, vcat(projectors, gates), vcat(cbits, locs))
end

"""
    general_c1_gates(num_bit::Int, projector::AbstractMatrix, cbit::Int, gates::Vector{AbstractMatrix}, locs::Vector{Int}) -> AbstractMatrix

general (low performance) construction method for control gate on different lines.
"""
general_c1_gates(
    num_bit::Int,
    projector::Tp,
    cbit::Int,
    gates::Vector{Tg},
    locs::Vector{Int},
) where {Tg<:AbstractMatrix,Tp<:AbstractMatrix} =
    hilbertkron(num_bit, [IMatrix(2) - projector], [cbit]) +
    hilbertkron(num_bit, vcat([projector], gates), vcat([cbit], locs))

"""
    rotmat(M::AbstractMatrix, θ::Real)

Returns rotated `M`: ``exp(-\\frac{imθ}{2} M)``.
"""
rotmat(M::AbstractMatrix, θ::Real) = exp(-im * θ / 2 * M)

"""
    linop2dense([T=ComplexF64], linear_map!::Function, n::Int) -> Matrix

Returns the dense matrix representation given linear map function.
"""
linop2dense(linear_map!::Function, n::Int) = linop2dense(ComplexF64, linear_map!, n)
linop2dense(::Type{T}, linear_map!::Function, n::Int) where {T} =
    linear_map!(Matrix{T}(I, 1 << n, 1 << n))

################### Fidelity ###################

"""
    density_fidelity(ρ1, ρ2)

General fidelity (including mixed states) between two density matrix for qubits.

# Definition

```math
F(ρ, σ)^2 = tr(ρσ) + 2 \\sqrt{det(ρ)det(σ)}
```
"""
function density_fidelity(ρ1::AbstractMatrix, ρ2::AbstractMatrix)
    return sqrt(tr(ρ1 * ρ2) + 2 * sqrt(det(ρ1) * det(ρ2)))
end

"""
    pure_state_fidelity(v1::Vector, v2::Vector)

fidelity for pure states.
"""
pure_state_fidelity(v1::Vector, v2::Vector) = abs(v1' * v2)

"""
    purification_fidelity(m1::Matrix, m2::Matrix)

Fidelity for mixed states via purification.

Reference:
    http://iopscience.iop.org/article/10.1088/1367-2630/aa6a4b/meta
"""
function purification_fidelity(m1::Matrix, m2::Matrix)
    O = m1' * m2
    return tr(sqrt(O * O'))
end

"""
    rand_unitary([T=ComplexF64], N::Int) -> Matrix

Create a random unitary matrix.
"""
rand_unitary(N::Int) = rand_unitary(ComplexF64, N)
rand_unitary(::Type{T}, N::Int) where {T} = qr(randn(T, N, N)).Q |> Matrix

"""
    sprand_unitary([T=ComplexF64], N::Int, density) -> SparseMatrixCSC

Create a random sparse unitary matrix.
"""
sprand_unitary(N::Int, density::Real) = sprand_unitary(ComplexF64, N, density)
sprand_unitary(::Type{T}, N::Int, density::Real) where {T} =
    SparseMatrixCSC(qr(sprandn(T, N, N, density)).Q)

"""
    rand_hermitian([T=ComplexF64], N::Int) -> Matrix

Create a random hermitian matrix.
"""
rand_hermitian(N::Int) = rand_hermitian(ComplexF64, N)

function rand_hermitian(::Type{T}, N::Int) where {T}
    A = randn(T, N, N)
    A + A'
end

"""
    sprand_hermitian([T=ComplexF64], N, density)

Create a sparse random hermitian matrix.
"""
sprand_hermitian(N::Int, density) = sprand_hermitian(ComplexF64, N, density)

function sprand_hermitian(::Type{T}, N::Int, density::Real) where {T}
    A = sprandn(T, N, N, density)
    return A + A'
end

@static if VERSION < v"1.1.0"
    function SparseArrays.sprandn(::Type{T}, n::Int, m::Int, density::Real) where {T<:Real}
        return T.(sprandn(n, m, density))
    end

    function SparseArrays.sprandn(
        ::Type{Complex{T}},
        n::Int,
        m::Int,
        density::Real,
    ) where {T<:Real}
        return T.(sprandn(n, m, density)) + im * T.(sprandn(n, m, density))
    end
end

"""
    autostatic(A[; threshold=8])

Staticize dynamic array `A` by a `threshold`.
"""
autostatic(A::AbstractVecOrMat; threshold::Int = 8) =
    length(A) > (1 << threshold) ? A : staticize(A)

# General definition
function rot_mat(::Type{T}, gen::AbstractMatrix, theta::Real) where {N,T}
    I = IMatrix{size(gen, 1),T}()
    m = I * cos(theta / 2) - im * sin(theta / 2) * gen
    if eltype(m) != T
        m2 = similar(m, T)
        copyto!(m2, m)
        return m2
    else
        return m
    end
end

BitBasis.unsafe_reorder(A::IMatrix, orders::NTuple{N,<:Integer}) where {N} = A

function BitBasis.unsafe_reorder(A::PermMatrix, orders::NTuple{N,<:Integer}) where {N}
    od = Vector{Int}(undef, 1 << length(orders))
    for (i, b) in enumerate(ReorderedBasis(orders))
        @inbounds od[i] = 1 + b
    end

    perm = similar(A.perm)
    vals = similar(A.vals)

    @simd for i = 1:length(perm)
        @inbounds perm[od[i]] = od[A.perm[i]]
        @inbounds vals[od[i]] = A.vals[i]
    end

    return PermMatrix(perm, vals)
end


text:
module Playground
using Base: Cartesian
typealias Feature UInt8

# Enumerates possible features on the grid:
# - FLUID: Site on which to perform standard lattice-boltzmann
# - SOLID: Wall or boundary of some kind
# - INLET: input into the simulation domain
# - OUTLET: output out of the simulation domain
# - NOTHING: convenience value when constructing grid
# - NULL: means streaming to /dev/null
const NOTHING = convert(Feature, 1)
const FLUID   = convert(Feature, 2)
const SOLID   = convert(Feature, 3)
const INLET   = convert(Feature, 4)
const OUTLET  = convert(Feature, 5)
const NULL    = convert(Feature, 6)
const NUMBER_OF_FEATURE_TYPES = convert(Feature, 6)

for N = 2:3
    @eval begin
        function initialize(characterize::Function, grid::Array{Feature, $N})
            Cartesian.@nloops $N i grid begin
                site::Feature = characterize((Cartesian.@ntuple $N i)...)
                if site != NOTHING
                    (Cartesian.@nref $N grid i) = site
                end
            end
            grid
        end
    end
end

# Goes over the playground grid, using the function to set it up
# The function is passed the indices to the playground grid for each site in
# the grid.
# It should return an instance of Feature
# If the return is equal to NOTHING, then that site is left untouched.
# Otherwise, the site is set to that value.
# This mechanism makes it possible to call initialize multiple times with
# different functions.
initialize(characterize::Function, dimensions) = initialize(
    characterize,
    zeros(Feature, dimensions)
)

end

ext: jl


text:
############################################################
# judiWavefield ##############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: June 2017

export judiWavefield, fft, ifft

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

mutable struct judiWavefield{T} <: judiMultiSourceVector{T}
    nsrc::Integer
    dt::Vector{T}
    data::Vector{<:Union{Array{T, N}, PyArray}} where N
end

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

## outer constructors

"""
judiWavefield
        nsrc::Integer
        dt::AbstractFloat
        data

Abstract vector for seismic wavefields.

Constructors
============

Construct wavefield vector from an info structure, a cell array of wavefields and the computational \\
time step dt:

    judiWavefield(nsrc, dt, data)

"""
function judiWavefield(nsrc::Integer, dt::AbstractFloat, data::Array{T, N}) where {T<:Number, N}
	# length of vector
	dataCell = [convert(Array{Float32, N}, data) for j=1:nsrc]
	return judiWavefield{Float32}(nsrc, [Float32(dt) for i=1:nsrc], dataCell)
end

function judiWavefield(dt::AbstractFloat, data::Vector{Array{T, N}}) where {T, N}
	# length of vector
	nsrc = length(data)
	T != Float32 && (data = tof32.(data))
	return judiWavefield{Float32}(nsrc, [Float32(dt) for i=1:nsrc], data)
end

judiWavefield(dt::AbstractFloat, nsrc::Integer, data::Array{T, N}) where {T<:Number, N} = judiWavefield(nsrc, dt, data)
judiWavefield(dt::AbstractFloat, data::Vector{Any}) = judiWavefield(dt, tof32.(data))
judiWavefield(dt::AbstractFloat, data::Array{T, N}) where {T<:Number, N} = judiWavefield(1, dt, data)

conj(w::judiWavefield{T}) where {T<:Complex} = judiWavefield{R}(w.nsrc, w.dt, conj(w.data))

############################################################
## overloaded multi_source functions
time_sampling(jv::judiWavefield) = jv.dt

####################################################################
# JOLI conversion
jo_convert(::Type{T}, jw::judiWavefield{T}, ::Bool) where {T<:Number} = jw
jo_convert(::Type{T}, jw::judiWavefield{vT}, B::Bool) where {T<:Number, vT} = judiWavefield{T}(jw.nsrc, jv.dt, jo_convert.(T, jw.data, B))
zero(::Type{T}, v::judiWavefield{vT}; nsrc::Integer=v.nsrc) where {T, vT} = judiWavefield{T}(nsrc, v.dt, T(0) .* v.data[1:nsrc])
zero(::Type{T}, v::judiWavefield{vT}; nsrc::Integer=v.nsrc) where {T<:AbstractFloat, vT<:Complex} = judiWavefield{T}(nsrc, v.dt, T(0) .* real(v.data[1:nsrc]))
(w::judiWavefield)(x::Vector{<:Array}) = judiWavefield(w.dt, x)

function copy!(jv::judiWavefield, jv2::judiWavefield)
    v.data .= jv2.data
    jv.dt = jv2.dt
    jv
end

copyto!(jv::judiWavefield, jv2::judiWavefield) = copy!(jv, jv2)
make_input(w::judiWavefield) = w.data[1]
check_compat(ms::Vararg{judiWavefield, N}) where N = all(y -> y.dt == first(ms).dt, ms)

getindex(a::judiWavefield{T}, srcnum::RangeOrVec) where T = judiWavefield{T}(length(srcnum), a.dt[srcnum], a.data[srcnum])
####################################################################

function push!(a::judiWavefield{T}, b::judiWavefield{T}) where T
    append!(a.data, b.data)
    append!(a.dt, b.dt)
    a.nsrc += b.nsrc
end

# DFT operator for wavefields, acts along time dimension
function fft(x_in::judiWavefield{T}) where T
    x = similar(x_in, Complex{Float32})
    for i=1:x_in.nsrc
        x.data[i] = fft(x_in.data[i], 1)/sqrt(size(x_in.data[i], 1))
    end
    return x
end

function ifft(x_in::judiWavefield{T}) where T
    x = similar(x_in, Float32)
    for i=1:x_in.nsrc
        x.data[i] = real(ifft(x_in.data[i], 1)) * sqrt(size(x_in.data[i], 1))
    end
    return x
end


text:
using RecordedArrays
using Random
using EvolutionaryModelingTools
using EvolutionaryModelingTools.Scalar

# parameters
const T = 10.0 # max time
const r = 0.5   # growth rate
const d = 0.1   # death rate
const K = 100   # carrying capacity
const n = 10    # initial population size

@cfunc growth_c(r, n) = r * n
@ufunc growth_u!(ind, n) = n[ind] += 1

@cfunc death_c(d, n) = d * n
@ufunc death_u!(ind, n) = n[ind] -= 1

@cfunc comp_c(r, K, n) = r * n * n' / K
@ufunc comp_u!(ind::CartesianIndex{0}, n) = n[ind] -= 1

# the ind must be CartesianIndex{0}, this is just for test
@ufunc check_extinction(ind::CartesianIndex{0}, n) = n[ind] == 0 ? :extinct : :finnish

growth = Reaction(growth_c, growth_u!)
death = Reaction(death_c, death_u!)
competition = Reaction(comp_c, comp_u!)

_maybe_scalar(n::Number) = scalar(n)
_maybe_scalar(n) = n

function run_gillespie(rng=Random.GLOBAL_RNG, T=T, r=r, d=d, K=K, n=n)
    ps = (; r, d, K, n=_maybe_scalar(n))
    rs = (growth, death, competition)
    return gillespie(check_extinction, rng, T, ps, rs)
end

function run_gillespie_record(rng=Random.GLOBAL_RNG, T=T, r=r, d=d, K=K, n=n)
    clock = ContinuousClock(T)
    ps = (; r, d, K, n= recorded(DynamicEntry, clock, n))
    rs = (growth, death, competition)
    return gillespie(check_extinction, rng, clock, ps, rs)
end

function run_manually(rng=Random.GLOBAL_RNG, T=T, r=r, d=d, K=K, n=n)
    t = 0
    while t <= T
        # calculate rate
        growth_rate = r * n              # intrinsic growth
        death_rate = d * n               # intrinsic death
        competition_rate = r * n * n / K # resource competition
        summed = growth_rate + death_rate + competition_rate  # sum a_i
        # break if summed is zero
        summed == 0 && break
        # update current time
        t += -log(rand(rng)) / summed
        # sample a reaction and adjust population size
        rn = rand(rng) * summed
        if rn < growth_rate
            n += 1
        else
            n -= 1
        end
    end
    return n
end

function run_manually_record(rng=Random.GLOBAL_RNG, T=T, r=r, d=d, K=K, n=n)
    clock = ContinuousClock(T)
    n = recorded(DynamicEntry, clock, n)
    for _ in clock
        # calculate rate
        growth_rate = r * n              # intrinsic growth
        death_rate = d * n               # intrinsic death
        competition_rate = r * n * n / K # resource competition
        summed = growth_rate + death_rate + competition_rate  # sum a_i
        # break if summed is zero
        summed == 0 && break
        # update current time
        increase!(clock, -log(rand(rng)) / summed)
        # sample a reaction and adjust population size
        rn = rand(rng) * summed
        if rn < growth_rate
            n[] += 1
        else
            n[] -= 1
        end
    end
    return n
end


text:
"""
```julia
pa=SeisInvExpt(attrib_mod,attrib_inv, attrib; rfields, parameterization)
```
Predefined gallery of `SeisInvExpt`. Choose `attrib::Symbol`
* `=:pizza` an experiment fast enough to be run on a laptop
* `=:downhole` sources and receivers on a drill-string
"""
function SeisInvExpt(attrib_mod::Union{FdtdAcou, FdtdAcouBorn},attrib_inv::Union{LS,LS_prior,Migr,Migr_FD}, attrib::Symbol;
		     rfields=[:vx,:vz,:p], α=0.0, parameterization=[:χvp, :χrho, :null])

	if(attrib==:pizza)
		# starting model
		model = Medium(:pizza);
		model0 = similar(model)

		ageom=AGeom(model.mgrid, SSrcs(5), Srcs(1), Recs(20))
		update!(ageom, SSrcs(),[0,0],900.0,[0, 2π])
		update!(ageom, Recs(),[0,0],900.0,[0, 2π])

		wav, tgrid=ricker(model, 3, 1.0)
		srcwav = SrcWav(tgrid, ageom, [:p])
		update!(srcwav, [:p], wav)

		igrid=broadcast(x->range(x[1],stop=x[end],step=50.),model.mgrid)
		igrid_interp_scheme=:B2
	elseif(attrib==:downhole)
		gri = [range(-120.,stop=30.,step=1.0), range(-30.,stop=30.,step=1.0)]
		model = Medium(gri)
		update!(model, [:vp,:rho], [[2500., 3500.], [2500., 3500.]])
		fill!(model)
		model0=deepcopy(model)
		print(model)
		update!(model, [:vp,:rho], randn_perc=0.1)
		update!(model0, [:vp,:rho], randn_perc=0.1)
		update!(model, [:vp], ellip_rad=[2000., 5.], ellip_loc=[0.,0.],ellip_pert=0.1, ellip_α=α)

		ageom=AGeom(model.mgrid, SSrcs(2), Srcs(1), Recs(20))

		update!(ageom, SSrcs(),[-75,0],[-50,0])
		update!(ageom, Recs(),[-100,0],[-50,0])

		wav, tgrid=ricker(model, 8, 0.8)
		srcwav = SrcWav(tgrid, ageom, [:p])
		update!(srcwav, [:p], wav)

		igrid=[range(-40.,stop=30.,length=100), range(-30.,length=1,step=60.)]
		igrid_interp_scheme=:B2
	else
		error("invalid attrib")
	end

	return  SeisInvExpt(attrib_mod, attrib_inv, srcwav=srcwav, ageom=ageom, tgrid=tgrid,
		       rfields=rfields,
		       mediumm=model0, mediumm0=model0, mediumm_obs=model,
		       igrid_interp_scheme=igrid_interp_scheme,
		       igrid=igrid, parameterization=parameterization, verbose=false);
end


text:
using Random
using Test
using Gaugefields
using Wilsonloop

function gradientflow_test_4D(NX,NY,NZ,NT,NC)
    Dim = 4
    Nwing = 1

    Random.seed!(123)

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end
    factor = 1/(comb*U[1].NV*U[1].NC)

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("0 plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("0 polyakov loop = $(real(poly)) $(imag(poly))")

    #Plaquette term
    loops_p = Wilsonline{Dim}[]
    for μ=1:Dim
        for ν=μ:Dim
            if ν == μ
                continue
            end
            loop1 = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)],Dim = Dim)
            push!(loops_p,loop1)
            #loop1 = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)],Dim = Dim)
            #push!(loops_p,loop1)
        end
    end

    #Rectangular term
    loops = Wilsonline{Dim}[]
    for μ=1:Dim
        for ν=μ:Dim
            if ν == μ
                continue
            end
            loop1 = Wilsonline([(μ,1),(ν,2),(μ,-1),(ν,-2)],Dim = Dim)
            push!(loops,loop1)
            loop1 = Wilsonline([(μ,2),(ν,1),(μ,-2),(ν,-1)],Dim = Dim)

            push!(loops,loop1)
        end
    end

    listloops = [loops_p,loops]
    listvalues = [1+im,0.1]
    g = Gradientflow_general(U,listloops,listvalues,eps = 0.01)

    for itrj=1:100
        flow!(U,g)
        if itrj % 10 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end
    return plaq_t

end

function gradientflow_test_2D(NX,NT,NC)
    Dim = 2
    Nwing = 1

    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("0 plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("0 polyakov loop = $(real(poly)) $(imag(poly))")

    #g = Gradientflow(U,eps = 0.01)
    #listnames = ["plaquette"]
    #listvalues = [1]
    loops_p = Wilsonline{Dim}[]
    for μ=1:Dim
        for ν=μ:Dim
            if ν == μ
                continue
            end

            loop1 = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)],Dim = Dim)

            push!(loops_p,loop1)
            #loop1 = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)],Dim = Dim)
            #push!(loops_p,loop1)
        end
    end

    loops = Wilsonline{Dim}[]
    for μ=1:Dim
        for ν=μ:Dim
            if ν == μ
                continue
            end
            loop1 = Wilsonline([(μ,1),(ν,2),(μ,-1),(ν,-2)],Dim = Dim)
            push!(loops,loop1)
            loop1 = Wilsonline([(μ,2),(ν,1),(μ,-2),(ν,-1)],Dim = Dim)

            push!(loops,loop1)
        end
    end

    listloops = [loops_p,loops]
    listvalues = [1+im,0.1] #complex coefficient
    g = Gradientflow_general(U,listloops,listvalues,eps = 0.01)
    #g = Gradientflow(U,eps = 0.01)

    for itrj=1:100
        flow!(U,g)
        if itrj % 10 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

    return plaq_t

end

#const eps = 0.1

println("2D system")
@testset "2D" begin
    NX = 4
    #NY = 4
    #NZ = 4
    NT = 4
    Nwing = 0

    @testset "NC=1" begin
        β = 2.3
        NC = 1
        println("NC = $NC")
        #val =0.6414596466929057
        val =   0.9993254431181984
        @time plaq_t = gradientflow_test_2D(NX,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end
    #error("d")

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.6414596466929057
        val =  0.9768786716327604
        @time plaq_t = gradientflow_test_2D(NX,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        #val = 0.5779454661484242
        val = 0.9656356864814539
        @time plaq_t = gradientflow_test_2D(NX,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.19127260002797497
        val = 0.8138836242603148
        @time plaq_t = gradientflow_test_2D(NX,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end
end

println("4D system")
@testset "4D" begin
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.6414596466929057
        #val = 0.5920897445000382
        val  =0.9440125563836135
        @time plaq_t = gradientflow_test_4D(NX,NY,NZ,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        #val = 0.5779454661484242
        #val  =0.9440125563836135
        #val = 0.5385142466966718
        val = 0.8786515255315753
        @time plaq_t = gradientflow_test_4D(NX,NY,NZ,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.19127260002797497
        #val = 0.1904815857904191
        val = 0.7301232810349298
        @time plaq_t =gradientflow_test_4D(NX,NY,NZ,NT,NC)
        #@test abs(plaq_t-val)/abs(val) < eps
    end

end


text:
@testset "Language" begin
    a = c_const("a")
    b = c_const("b")

    @test isa(a, Constant)
    @test isa(b, Constant)

    f = c_functor("f", 1)
    p = c_pred("p", 2)

    @test isa(f, Functor)
    @test isa(p, Predicate)

    s = f(a)
    @test isa(s, Structure)

    lit = p(a,b)
    @test isa(lit, Literal)

    vX = c_var("X")
    vY = c_var("Y")

    q = c_pred("q", 2)

    h = q(vX, vY)
    @test isa(h, Literal)

    clause = p(vX) <= q(vX, vY)
    @test isa(clause, Clause)

    clause2 = p(vX) <= q(vX, vY) & p(vY)
    @test isa(clause2, Clause)

    n = c_pred("n", 1)
    lit_n1 = n(1)
    @test isa(lit_n1, Literal)
end


text:
module Keys

export Key

import SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

@enum Key::UInt32 begin
  UNKNOWN = SDL2.SDLK_UNKNOWN;
  BACKSPACE = SDL2.SDLK_BACKSPACE;
  TAB = SDL2.SDLK_TAB;
  CLEAR = SDL2.SDLK_CLEAR;
  RETURN = SDL2.SDLK_RETURN;
  PAUSE = SDL2.SDLK_PAUSE;
  ESCAPE = SDL2.SDLK_ESCAPE;
  SPACE = SDL2.SDLK_SPACE;
  EXCLAIM = SDL2.SDLK_EXCLAIM;
  QUOTEDBL = SDL2.SDLK_QUOTEDBL;
  HASH = SDL2.SDLK_HASH;
  DOLLAR = SDL2.SDLK_DOLLAR;
  AMPERSAND = SDL2.SDLK_AMPERSAND;
  QUOTE = SDL2.SDLK_QUOTE;
  LEFTPAREN = SDL2.SDLK_LEFTPAREN;
  RIGHTPAREN = SDL2.SDLK_RIGHTPAREN;
  ASTERISK = SDL2.SDLK_ASTERISK;
  PLUS = SDL2.SDLK_PLUS;
  COMMA = SDL2.SDLK_COMMA;
  MINUS = SDL2.SDLK_MINUS;
  PERIOD = SDL2.SDLK_PERIOD;
  SLASH = SDL2.SDLK_SLASH;
  K_0 = SDL2.SDLK_0;
  K_1 = SDL2.SDLK_1;
  K_2 = SDL2.SDLK_2;
  K_3 = SDL2.SDLK_3;
  K_4 = SDL2.SDLK_4;
  K_5 = SDL2.SDLK_5;
  K_6 = SDL2.SDLK_6;
  K_7 = SDL2.SDLK_7;
  K_8 = SDL2.SDLK_8;
  K_9 = SDL2.SDLK_9;
  COLON = SDL2.SDLK_COLON;
  SEMICOLON = SDL2.SDLK_SEMICOLON;
  LESS = SDL2.SDLK_LESS;
  EQUALS = SDL2.SDLK_EQUALS;
  GREATER = SDL2.SDLK_GREATER;
  QUESTION = SDL2.SDLK_QUESTION;
  AT = SDL2.SDLK_AT;
  LEFTBRACKET = SDL2.SDLK_LEFTBRACKET;
  BACKSLASH = SDL2.SDLK_BACKSLASH;
  RIGHTBRACKET = SDL2.SDLK_RIGHTBRACKET;
  CARET = SDL2.SDLK_CARET;
  UNDERSCORE = SDL2.SDLK_UNDERSCORE;
  BACKQUOTE = SDL2.SDLK_BACKQUOTE;
  A = SDL2.SDLK_a;
  B = SDL2.SDLK_b;
  C = SDL2.SDLK_c;
  D = SDL2.SDLK_d;
  E = SDL2.SDLK_e;
  F = SDL2.SDLK_f;
  G = SDL2.SDLK_g;
  H = SDL2.SDLK_h;
  I = SDL2.SDLK_i;
  J = SDL2.SDLK_j;
  K = SDL2.SDLK_k;
  L = SDL2.SDLK_l;
  M = SDL2.SDLK_m;
  N = SDL2.SDLK_n;
  O = SDL2.SDLK_o;
  P = SDL2.SDLK_p;
  Q = SDL2.SDLK_q;
  R = SDL2.SDLK_r;
  S = SDL2.SDLK_s;
  T = SDL2.SDLK_t;
  U = SDL2.SDLK_u;
  V = SDL2.SDLK_v;
  W = SDL2.SDLK_w;
  X = SDL2.SDLK_x;
  Y = SDL2.SDLK_y;
  Z = SDL2.SDLK_z;
  DELETE = SDL2.SDLK_DELETE;
  KP0 = SDL2.SDLK_KP_0;
  KP1 = SDL2.SDLK_KP_1;
  KP2 = SDL2.SDLK_KP_2;
  KP3 = SDL2.SDLK_KP_3;
  KP4 = SDL2.SDLK_KP_4;
  KP5 = SDL2.SDLK_KP_5;
  KP6 = SDL2.SDLK_KP_6;
  KP7 = SDL2.SDLK_KP_7;
  KP8 = SDL2.SDLK_KP_8;
  KP9 = SDL2.SDLK_KP_9;
  KP_PERIOD = SDL2.SDLK_KP_PERIOD;
  KP_DIVIDE = SDL2.SDLK_KP_DIVIDE;
  KP_MULTIPLY = SDL2.SDLK_KP_MULTIPLY;
  KP_MINUS = SDL2.SDLK_KP_MINUS;
  KP_PLUS = SDL2.SDLK_KP_PLUS;
  KP_ENTER = SDL2.SDLK_KP_ENTER;
  KP_EQUALS = SDL2.SDLK_KP_EQUALS;
  UP = SDL2.SDLK_UP;
  DOWN = SDL2.SDLK_DOWN;
  RIGHT = SDL2.SDLK_RIGHT;
  LEFT = SDL2.SDLK_LEFT;
  INSERT = SDL2.SDLK_INSERT;
  HOME = SDL2.SDLK_HOME;
  END = SDL2.SDLK_END;
  PAGEUP = SDL2.SDLK_PAGEUP;
  PAGEDOWN = SDL2.SDLK_PAGEDOWN;
  F1 = SDL2.SDLK_F1;
  F2 = SDL2.SDLK_F2;
  F3 = SDL2.SDLK_F3;
  F4 = SDL2.SDLK_F4;
  F5 = SDL2.SDLK_F5;
  F6 = SDL2.SDLK_F6;
  F7 = SDL2.SDLK_F7;
  F8 = SDL2.SDLK_F8;
  F9 = SDL2.SDLK_F9;
  F10 = SDL2.SDLK_F10;
  F11 = SDL2.SDLK_F11;
  F12 = SDL2.SDLK_F12;
  F13 = SDL2.SDLK_F13;
  F14 = SDL2.SDLK_F14;
  F15 = SDL2.SDLK_F15;
  NUMLOCK = SDL2.SDLK_NUMLOCKCLEAR;
  CAPSLOCK = SDL2.SDLK_CAPSLOCK;
  SCROLLOCK = SDL2.SDLK_SCROLLLOCK;
  RSHIFT = SDL2.SDLK_RSHIFT;
  LSHIFT = SDL2.SDLK_LSHIFT;
  RCTRL = SDL2.SDLK_RCTRL;
  LCTRL = SDL2.SDLK_LCTRL;
  RALT = SDL2.SDLK_RALT;
  LALT = SDL2.SDLK_LALT;
  # RMETA = SDL2.SDLK_RGUI;
  # LMETA = SDL2.SDLK_LGUI;
  # LSUPER = SDL2.SDLK_LGUI;
  # RSUPER = SDL2.SDLK_RGUI;
  MODE = SDL2.SDLK_MODE;
  HELP = SDL2.SDLK_HELP;
  SYSREQ = SDL2.SDLK_SYSREQ;
  MENU = SDL2.SDLK_MENU;
  POWER = SDL2.SDLK_POWER;
  EURO = SDL2.SDLK_CURRENCYUNIT;
end

end

using .Keys


text:
module FFTLog

using FFTW
using Base: @kwdef
using SpecialFunctions: gamma
import Base: *

export prepare_FFTLog!, evaluate_FFTLog, evaluate_FFTLog!
export prepare_Hankel!, evaluate_Hankel, evaluate_Hankel!
export mul!
export get_y

abstract type AbstractPlan end

set_num_threads(N::Int) = FFTW.set_num_threads(N)

"""
    _c_window(N::AbstractArray, NCut::Int)

This function returns the smoothing window function as defined
in Eq. (C.1) of [McEwen et al. (2016)](https://arxiv.org/abs/1603.04826).
"""
function _c_window(N::AbstractArray, NCut::Int)
    NRight = last(N) - NCut
    NR = filter(x->x>=NRight, N)
    ThetaRight = (last(N).-NR) ./ (last(N) - NRight - 1)
    W = ones(length(N))
    W[findall(x->x>=NRight, N)] = ThetaRight .- 1 ./ (2*π) .* sin.(2 .* π .*
	ThetaRight)
    return W
end

"""
    FFTLogPlan()

This struct contains all the elements necessary to evaluate the integral with one Bessel function.
"""
@kwdef mutable struct FFTLogPlan{T,C} <: AbstractPlan
    x::Vector{T}
    y::Matrix{T} = zeros(10,10)
    fy::Matrix{T} = zeros(10,10)
    hm::Matrix{C} = zeros(ComplexF64,10,10)
    hm_corr::Matrix{C} = zeros(ComplexF64,10,10)
    d_ln_x::T = log(x[2]/x[1])
    fy_corr::Matrix{T} = zeros(10,10)
    original_length::Int = length(x)
    gl::Matrix{C} = zeros(ComplexF64,100,100)
    ν::T = 1.01
    n_extrap_low::Int = 0
    n_extrap_high::Int = 0
    c_window_width::T = 0.25
    n_pad::Int = 0
    n::Int = 0
    N::Int = original_length+n_extrap_low+n_extrap_high+2*n_pad
    m::Vector{T} = zeros(N)
    cm::Vector{C} = zeros(ComplexF64,N)
    ηm::Vector{T} = zeros(N)
    plan_rfft::FFTW.rFFTWPlan = plan_rfft(randn(1024))
    plan_irfft = plan_irfft(randn(Complex{Float64}, 2,
    Int((original_length+n_extrap_low+n_extrap_high+2*n_pad)/2) +1),
    original_length+n_extrap_low+n_extrap_high+2*n_pad, 2)
end

@kwdef mutable struct HankelPlan{T,C} <: AbstractPlan
    x::Vector{T}
    y::Matrix{T} = zeros(10,10)
    fy::Matrix{T} = zeros(10,10)
    hm::Matrix{C} = zeros(ComplexF64,10,10)
    hm_corr::Matrix{C} = zeros(ComplexF64,10,10)
    d_ln_x::T = log(x[2]/x[1])
    fy_corr::Matrix{T} = zeros(10,10)
    original_length::Int = length(x)
    gl::Matrix{C} = zeros(ComplexF64,100,100)
    ν::T = 1.01
    n_extrap_low::Int = 0
    n_extrap_high::Int = 0
    c_window_width::T = 0.25
    n_pad::Int = 0
    n::Int = 0
    N::Int = original_length+n_extrap_low+n_extrap_high+2*n_pad
    m::Vector{T} = zeros(N)
    cm::Vector{C} = zeros(ComplexF64,N)
    ηm::Vector{T} = zeros(N)
    plan_rfft::FFTW.rFFTWPlan = plan_rfft(randn(1024))
    plan_irfft = plan_irfft(randn(Complex{Float64}, 2,
    Int((original_length+n_extrap_low+n_extrap_high+2*n_pad)/2) +1),
    original_length+n_extrap_low+n_extrap_high+2*n_pad, 2)
end

function _eval_cm!(plan::AbstractPlan, fx)
    plan.cm = plan.plan_rfft * (fx .* plan.x .^ (-plan.ν))
    plan.cm .*= _c_window(plan.m, floor(Int, plan.c_window_width*plan.N/2))
end

function _eval_ηm!(plan::AbstractPlan)
    #TODO: #9 since we know the length of the initial array, we could use this info here to
    #remove the necessity of DLnX
    plan.ηm = 2 .* π ./ (plan.N .* plan.d_ln_x) .* plan.m
end

function _eval_gl(ell, z::Vector, n::Int)
    gl = ((-1)^n) .* 2 .^ (z .-n) .* gamma.((ell .+ z .- n)/2) ./
    gamma.((3 .+ ell .+ n .- z)/2)
    if n != 0
        for i in 1:n
            gl .*= (z .- i)
        end
    end
    return gl
end

function _logextrap(x::Vector, n_extrap_low::Int, n_extrap_high::Int)
    d_ln_x_low = log(x[2]/x[1])
    d_ln_x_high= log(reverse(x)[1]/reverse(x)[2])
    if n_extrap_low != 0
        X = vcat(x[1] .* exp.(d_ln_x_low .* Array(-n_extrap_low:-1)), x)
    end
    if n_extrap_high != 0
        X = vcat(X, last(X) .* exp.(d_ln_x_high .* Array(1:n_extrap_high)))
    end
    return X
end

function _zeropad(x::Vector, n_pad::Int)
    return vcat(zeros(n_pad), x, zeros(n_pad))
end

function _eval_y!(plan::AbstractPlan, ell::Vector)
    reverse!(plan.x)

    plan.y = zeros(length(ell), length(plan.x))
    plan.fy_corr = zeros(size(plan.y))

    @inbounds for myl in 1:length(ell)
        plan.y[myl,:] = (ell[myl] + 1) ./ plan.x
        plan.fy_corr[myl,:] =
        plan.y[myl,:] .^ (-plan.ν) .* sqrt(π) ./4
    end

    plan.fy = zeros(size(plan.y))

    reverse!(plan.x)
end

function _eval_gl_hm!(plan::AbstractPlan, ell::Vector)
    z = plan.ν .+ im .* plan.ηm
    plan.gl = zeros(length(ell), length(z))

    plan.hm_corr = zeros(ComplexF64, length(ell), Int(plan.N/2+1))

    @inbounds for myl in 1:length(ell)
        plan.hm_corr[myl,:] =
        (plan.x[1] .* plan.y[myl,1] ) .^ (-im .*plan.ηm)
        plan.gl[myl,:] = _eval_gl(ell[myl], z, plan.n)
    end
    plan.hm = zeros(ComplexF64, size(plan.gl))

end

function prepare_FFTLog!(plan::AbstractPlan, ell::Vector)
    plan.x = _logextrap(plan.x, plan.n_extrap_low + plan.n_pad,
	  plan.n_extrap_high + plan.n_pad)

    _eval_y!(plan, ell)

    plan.m = Array(0:length(plan.x)/2)
    _eval_ηm!(plan)

    _eval_gl_hm!(plan, ell)

    plan.plan_rfft = plan_rfft(plan.x)
    plan.plan_irfft = plan_irfft(randn(Complex{Float64}, length(ell),
    Int((length(plan.x)/2) +1)),
    plan.original_length+plan.n_extrap_low+plan.n_extrap_high+2*plan.n_pad, 2)
end

function prepare_Hankel!(hankplan::HankelPlan, ell::Vector)
    prepare_FFTLog!(hankplan, ell .-0.5)
end

function get_y(plan::AbstractPlan)
    return plan.y[:,plan.n_extrap_low+plan.n_pad+1:plan.n_extrap_low+
	plan.n_pad+plan.original_length]
end

function evaluate_FFTLog(plan::AbstractPlan, fx)
    fy = similar(get_y(plan))
    evaluate_FFTLog!(fy, plan, fx)
    return fy
end

function evaluate_FFTLog!(fy, plan::AbstractPlan, fx)
    fx = _logextrap(fx, plan.n_extrap_low,
	plan.n_extrap_high)
    fx = _zeropad(fx, plan.n_pad)
    _eval_cm!(plan, fx)

    @inbounds for myl in 1:length(plan.y[:,1])
        plan.hm[myl,:] = plan.cm .* @view plan.gl[myl,:]
        plan.hm[myl,:] .*= @view plan.hm_corr[myl, :]
    end

    plan.fy[:,:] .= plan.plan_irfft * conj!(plan.hm)
    plan.fy[:,:] .*= @view plan.fy_corr[:,:]

    fy[:,:] .= @view plan.fy[:,plan.n_extrap_low+plan.n_pad+
	1:plan.n_extrap_low+plan.n_pad+plan.original_length]
end

function evaluate_Hankel(hankplan::HankelPlan, fx)
    fy = evaluate_FFTLog(hankplan, fx .*(
        hankplan.x[hankplan.n_extrap_low+hankplan.n_pad+1:hankplan.n_extrap_low+
	hankplan.n_pad+hankplan.original_length]).^(5/2) )
    fy .*= sqrt.(2*get_y(hankplan)/π)
    return fy
end

function evaluate_Hankel!(fy, hankplan::HankelPlan, fx)
    evaluate_FFTLog!(fy, hankplan, fx .*(
        hankplan.x[hankplan.n_extrap_low+hankplan.n_pad+1:hankplan.n_extrap_low+
	hankplan.n_pad+hankplan.original_length]).^(5/2) )
    fy .*= sqrt.(2*get_y(hankplan)/π)
    return fy
end

function mul!(Y, Q::FFTLogPlan, A)
    evaluate_FFTLog!(Y, Q, A)
end

function mul!(Y, Q::HankelPlan, A)
    Y[:,:] .= evaluate_Hankel!(Y, Q, A)
end

function *(Q::FFTLogPlan, A)
    evaluate_FFTLog(Q, A)
end

function *(Q::HankelPlan, A)
    evaluate_Hankel(Q, A)
end

end # module


text:
#=
Simulation for ESSE

Ignacio Quintero Mächler

t(-_-t)

January 12 2017
=#

"""
    simulate_sse(λ       ::Array{Float64,1},
                 μ       ::Array{Float64,1},
                 l       ::Array{Float64,1},
                 g       ::Array{Float64,1},
                 q       ::Array{Float64,1},
                 β       ::Array{Float64,1},
                 x       ::Array{Float64,1},
                 y       ::Array{Float64,N},
                 cov_mod ::NTuple{M,String};
                 δt      ::Float64 = 1e-4,
                 nspp_min::Int64   = 1,
                 nspp_max::Int64   = 200_000,
                 retry_ext::Bool   = true))

Simulate tree according to EGeoHiSSE.
"""
function simulate_sse(λ       ::Array{Float64,1},
                      μ       ::Array{Float64,1},
                      l       ::Array{Float64,1},
                      g       ::Array{Float64,1},
                      q       ::Array{Float64,1},
                      β       ::Array{Float64,1},
                      x       ::Array{Float64,1},
                      y       ::Array{Float64,N},
                      cov_mod ::NTuple{M,String};
                      δt      ::Float64 = 1e-4,
                      nspp_min::Int64   = 1,
                      nspp_max::Int64   = 200_000,
                      retry_ext::Bool   = true) where {M,N}

  # make simulation
  ed, el, st, n, S, k =
    simulate_edges(λ, μ, l, g, q, β, x, y, δt, cov_mod, nspp_max)

  @info "Tree with $n species successfully simulated"

  if retry_ext
    in(0.0, el) && @warn "a lineage speciated at time 0.0..."
    while ed == 0 || n < nspp_min || n > nspp_max || in(0.0, el)
      ed, el, st, n, S, k =
        simulate_edges(λ, μ, l, g, q, β, x, y, δt, cov_mod, nspp_max)
      @info ("Tree with $n species successfully simulated")
    end
  else
    if ed == 0
      return 0, 0, 0
    end
  end

  # organize in postorder
  ed     = numberedges(ed, n)
  ed, el = postorderedges(ed, el, n)

  ## round branch lengths
  # find out order of simulation
  i = 1
  while !isone(δt*10^i)
    i += 1
  end

  el = map(x -> round(x; digits = i+1), el)

  # organize states
  tip_val = Dict(i => st[i] for i = 1:n)

  tip_val = states_to_values(tip_val, S, k)

  return tip_val, ed, el
end

"""
    simulate_edges(λ       ::Array{Float64,1},
                   μ       ::Array{Float64,1},
                   l       ::Array{Float64,1},
                   g       ::Array{Float64,1},
                   q       ::Array{Float64,1},
                   β       ::Array{Float64,1},
                   x       ::Array{Float64,1},
                   y       ::Array{Float64, N},
                   δt      ::Float64,
                   cov_mod ::NTuple{M,String},
                   nssp_max::Int64)

Simulate edges tree according to EGeoHiSSE.
"""
function simulate_edges(λ       ::Array{Float64,1},
                        μ       ::Array{Float64,1},
                        l       ::Array{Float64,1},
                        g       ::Array{Float64,1},
                        q       ::Array{Float64,1},
                        β       ::Array{Float64,1},
                        x       ::Array{Float64,1},
                        y       ::Array{Float64, N},
                        δt      ::Float64,
                        cov_mod ::NTuple{M,String},
                        nssp_max::Int64) where {M,N}

  h = Int64((sqrt(length(q)*4 + 1) + 1)/2)
  k = div(length(l), h)
  ny = size(y, 2)

  # if multidimensional
  md = ny > 1

  # check number of parameters and data
  model = id_mod(cov_mod, k, h, ny, λ, μ, l, g, q, β)

  # total simulation time
  simt = x[end]

  # make approximate time function
  af! = make_af(x, y, Val(size(y,2)))

  # preallocate af! result vector
  r = Array{Float64,1}(undef, size(y,2))

  # get function estimates at `simt`
  af!(simt, r)

  # create states
  S = create_states(k, h)

  ns = length(S)

  # areas & hidden states
  as = 1:k
  hs = 0:(h-1)

  # vector of initial state probabilities
  isp = Array{Float64,1}(undef, length(S))

  # assign initial state probabilities
  init_states_pr!(isp, λ, l, g, q, μ, β, S, r, k, h, ny, model, md, as, hs)

  # preallocate vector of individual area probabilities
  spr = Array{Float64,1}(undef, ns)

  # sample initial state
  si = prop_sample(spr, isp, ns)

  # n = 2 starting lineages (assume no stem branch)
  n = 2

  # edges alive
  ea = [1, 2]

  # edge array
  ed = zeros(Int64, nssp_max*2, 2)
  ed[ea,:] = [1 2;
              1 3]

  # edge lengths
  el = zeros(nssp_max*2)

  # make probability vectors
  λpr, μpr, gpr, qpr, Sλpr, Sμpr, Sgpr, Sqpr =
   event_probs(λ, μ, l, g, q, β, r, S, k, h, ny, model, md, δt, ns, as, hs)

  # make λ, μ, g & q probability functions
  updλpr! = make_updλpr!(λ, β, λpr, Sλpr, δt, k, model, md, ny)
  updμpr! = make_updμpr!(μ, l, β, μpr, Sμpr, δt, k, h, model, md, ny)
  updgpr! = make_updgpr!(g, β, gpr, Sgpr, δt, k, h, model, md, as, S, ny)
  updqpr! = make_updqpr!(Sqpr)

  # preallocate states probabilities for specific lengths of each event
  svλ = Array{Float64,1}[]
  svμ = Array{Float64,1}[]
  svg = Array{Float64,1}[]
  svq = Array{Float64,1}[]
  for i in Base.OneTo(ns)
    push!(svλ, zeros(length(λpr[i])))
    push!(svμ, zeros(length(μpr[i])))
    push!(svg, zeros(length(gpr[i])))
    push!(svq, zeros(length(qpr[i])))
  end

  # make state change vectors
  gtos, μtos, qtos, λtos = makecorresschg(gpr, μpr, qpr, λpr, S, as, hs, k, ns)

  # model first speciation event for daughter inheritance
  st = λtos[si][prop_sample(svλ[si], λpr[si], length(svλ[si]))]

  # random assignment to daughters
  if rand() < 0.5
    reverse!(st)
  end

  # start simulation
  while true

    # keep track of time
    simt -= δt

    # estimate `z(simt)`
    af!(simt, r)

    # one time step
    for i in Base.OneTo(n)

      # update edge length
      el[ea[i]] += δt

      sti = st[i]

      #=
        gain event?
      =#
      if rand() < updgpr!(sti, S[sti], r)
        # sample from the transition probabilities
        @inbounds st[i] = gtos[sti][prop_sample(svg[sti], gpr[sti], length(svg[sti]))]

        # no more events at this time
        continue
      end

      #=
        μ or loss event?
      =#
      if rand() < updμpr!(sti, S[sti], r)

        # if extinction
        if isone(length(S[sti].g))

          # update time in other extant lineages
          el[ea[(i+1):n]] .+= δt

          # node to remove
          nod = ed[ea[i],1]

          # top or bottom edge of node
          top = nod == ed[ea[i]+1, 1]

          # remove edge (only preserving persisting species)
          ned = findfirst(isequal(0), ed)[1]

          if ned == 3
            printstyled("What would you do if an endangered animal is eating an endangered plant? Sometimes nature is too cruel... \n",
              color=:light_red)
            printstyled("tree went extinct... rerun simulation \n",
              color=:light_red)
            return 0, 0, 0, 0, 0, 0
          end

          @views ed[ea[i]:(ned-1),:] = ed[(ea[i]+1):ned,:]

          # remove edge length
          @views el[ea[i]:(ned-1)] = el[(ea[i]+1):ned]

          # update alive lineages
          deleteat!(ea,i)
          ea[i:end] .-= 1

          # remove node and extra edge
          @views pr = findfirst(isequal(nod), ed[:,1])
          @views da = findfirst(isequal(nod), ed[:,2])

          if isnothing(da)
            da = 0
          end

          if da != 0
            ed[da,2] = ed[pr,2]
            el[da]  += el[pr]
          end

          # remove intervening edge
          @views ed[pr:(ned-2),:] = ed[(pr+1):(ned-1),:]

          # remove edge lengths of intervening node
          @views el[pr:(ned-2)] = el[(pr+1):(ned-1)]

          ## update alive lineages
          # is the remaining node terminal?
          if da == 0 || in(ed[da,2], ed[:,1])
            ea[i:end] .-= 1
          else
            ea[findfirst(isequal(pr),ea)] = da
            sort!(ea)
            if top
              ea[(i+1):end] .-= 1
            else
              ea[i:end]     .-= 1
            end
          end

          # update living species states
          deleteat!(st,i)

          # update n species
          n = lastindex(ea)

          # break loop
          break

        # if local extinction (state change)
        else
          @inbounds st[i] = μtos[sti][prop_sample(svμ[sti], μpr[sti], length(svμ[sti]))]

          # no more events at this time
          continue
        end

      end

      #=
          q event?
      =#
      if rand() < updqpr!(sti)
        @inbounds st[i] = qtos[sti][prop_sample(svq[sti], qpr[sti], length(svq[sti]))]

        # no more events at this time
        continue
      end

      #=
          λ event?
      =#
      if rand() < updλpr!(sti, S[sti], r)

        # update time in other extant lineages
        el[ea[(i+1):n]] .+= δt

        ### add new edges
        # start node
        ed[ea[end] + 1,1] = ed[ea[end] + 2,1] = ed[ea[i],2]

        # end nodes
        ed[ea[end] + 1,2] = maximum(ed)+1
        ed[ea[end] + 2,2] = maximum(ed)+1

        # update living edges
        push!(ea, ea[end]+1, ea[end]+2)
        deleteat!(ea, i)

        # update living states according to speciation event
        nst1, nst2 =
          λtos[sti][prop_sample(svλ[sti], λpr[sti], length(svλ[sti]))]

        # random assignment to daughters
        if rand() < 0.5
          push!(st,nst1,nst2)
        else
          push!(st,nst2,nst1)
        end

        deleteat!(st,i)

        # update number of species alive
        n = lastindex(ea)

        # no more events for any of the remaining lineages
        break
      end

    end

    simt < 0.0 && break

    if n == nssp_max
      @warn "more than $nssp_max species"
      return 0, 0, 0, n, 0, 0
    end
  end

  # remove 0s
  ed = ed[1:(2n-2),:]
  el = el[1:(2n-2)]

  return ed, el, st, n, S, k
end

"""
    id_mod(cov_mod::NTuple{N,String},
           k      ::Int64,
           h      ::Int64,
           nzt    ::Int64,
           λ      ::Array{Float64,1},
           μ      ::Array{Float64,1},
           l      ::Array{Float64,1},
           g      ::Array{Float64,1},
           q      ::Array{Float64,1},
           β      ::Array{Float64,1})

Check if number of parameters and `z(t)` functions
is consistent with specified model.
"""
function id_mod(cov_mod::NTuple{N,String},
                k      ::Int64,
                h      ::Int64,
                ny    ::Int64,
                λ      ::Array{Float64,1},
                μ      ::Array{Float64,1},
                l      ::Array{Float64,1},
                g      ::Array{Float64,1},
                q      ::Array{Float64,1},
                β      ::Array{Float64,1}) where {N}

  all_ok = true

  # check consistent length of parameter vectors proposed
   if (lastindex(λ) != (k+1)*h   ||
       lastindex(μ) != k*h       ||
       lastindex(l) != k*h       ||
       lastindex(g) != k*(k-1)*h ||
       lastindex(q) != h*(h-1))
     printstyled("$k areas with $h hidden states should have the following parameter lengths:
       λ = $((k+1)*h) parameters
       μ = $(k*h) parameters
       l = $(k*h) parameters
       g = $(k*(k-1)*h) parameter matrix
       q = $(h*(h-1)) parameters")
    error("Parameter vector lengths not consistent with the number of states")
  end

   model = [false, false, false]

  for m in cov_mod
    # if speciation
    if occursin(r"^[s|S][A-za-z]*", m)
      model[1] = true
    end
    # if extinction
    if occursin(r"^[e|E][A-za-z]*", m)
      model[2] = true
    end
    # if dispersal
    if occursin(r"^[t|T|r|R|q|Q][A-za-z]*", m)
      model[3] = true
    end
  end

  # beta length should be
  yppar = ny == 1 ? 1 : div(ny,
         model[1]*k +
         model[2]*k +
         model[3]*k*(k-1))

  betaE = model[1]*k*h*yppar + model[2]*k*h*yppar + model[3]*k*(k-1)*h*yppar

  if lastindex(β) == betaE
    printstyled("Number of parameters are consistent with specified model \n",
      color=:green)
  else
    all_ok = false
    printstyled("this covariate model with $k areas, $ny covariates and $h hidden states should have a β vector length of $betaE \n")
  end

  mexp = "$(model[1] ? "speciation," : "")$(model[2] ? "extinction," : "")$(model[3] ? "transition," : "")"
  mexp = replace(mexp, "," => ", ")
  mexp = mexp[1:(end-2)]

  if all_ok
    printstyled("Simulating $mexp covariate SSE model with $k states, $h hidden states and $ny covariates \n",
      color=:green)
  else
    error("Parameter and z(t) function number not consistent with specified model")
  end

  return tuple(model...)
end

"""
    init_states_pr!(isp  ::Array{Float64,1},
                    λ    ::Array{Float64,1},
                    l    ::Array{Float64,1},
                    g    ::Array{Float64,1},
                    q    ::Array{Float64,1},
                    μ    ::Array{Float64,1},
                    β    ::Array{Float64,1},
                    S    ::Array{Sgh,1},
                    k    ::Int64,
                    h    ::Int64,
                    ny   ::Int64,
                    model::NTuple{3,Bool},
                    md   ::Bool,
                    as   ::UnitRange{Int64},
                    hs   ::UnitRange{Int64})

 Estimate starting transition probabilities for initial states
"""
function init_states_pr!(isp  ::Array{Float64,1},
                         λ    ::Array{Float64,1},
                         l    ::Array{Float64,1},
                         g    ::Array{Float64,1},
                         q    ::Array{Float64,1},
                         μ    ::Array{Float64,1},
                         β    ::Array{Float64,1},
                         S    ::Array{Sgh,1},
                         r    ::Array{Float64,1},
                         k    ::Int64,
                         h    ::Int64,
                         ny   ::Int64,
                         model::NTuple{3,Bool},
                         md   ::Bool,
                         as   ::UnitRange{Int64},
                         hs   ::UnitRange{Int64})

  # expected number of covariates
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting values for models
  m2s = model[1]*k*h
  m3s = m2s + model[2]*k*h
  y2s = model[1]*yppar*k
  y3s = y2s + model[2]*yppar*k

  for si in Base.OneTo(length(S))

    s = S[si]

    # allopatric speciation rates
    sr = 0.0
    if length(s.g) != k
      sr += λ[(s.h+1)*(k+1)]
    end

    # loss rates
    lr = 0.0
    for i = setdiff(as, s.g)
      lr += l[s.h*k+i]
    end

    # gain rates
    if model[3]
      gr = 0.0
      if length(s.g) > 1
        for ta = s.g, fa = setdiff(s.g, ta)
          gr += expf(g[s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                     β[m3s + s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                     md ? r[y3s + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)] : r[1])
        end
      end
    else
      gr = 0.0
      if length(s.g) > 1
        for ta = s.g, fa = setdiff(s.g, ta)
          gr += g[s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)]
        end
      end
    end

    # hidden states
    hr = 0.0
    for fa = setdiff(hs, s.h)
        hr += q[fa*(h-1) + (s.h > fa ? s.h : s.h+1)]
    end

    # add all and hidden rates
    isp[si] = sr + lr + gr + hr
  end

  return nothing
end

"""
    event_probs(λ    ::Array{Float64,1},
                μ    ::Array{Float64,1},
                l    ::Array{Float64,1},
                g    ::Array{Float64,1},
                q    ::Array{Float64,1},
                β    ::Array{Float64,1},
                r    ::Array{Float64,1},
                S    ::Array{Sgh,1},
                k    ::Int64,
                h    ::Int64,
                ny   ::Int64,
                model::NTuple{3,Bool},
                md   ::Bool,
                δt   ::Float64,
                ns   ::Int64,
                as   ::UnitRange{Int64},
                hs   ::UnitRange{Int64})

Create event probabilities for each state given the input parameters.
"""
function event_probs(λ    ::Array{Float64,1},
                     μ    ::Array{Float64,1},
                     l    ::Array{Float64,1},
                     g    ::Array{Float64,1},
                     q    ::Array{Float64,1},
                     β    ::Array{Float64,1},
                     r    ::Array{Float64,1},
                     S    ::Array{Sgh,1},
                     k    ::Int64,
                     h    ::Int64,
                     ny   ::Int64,
                     model::NTuple{3,Bool},
                     md   ::Bool,
                     δt   ::Float64,
                     ns   ::Int64,
                     as   ::UnitRange{Int64},
                     hs   ::UnitRange{Int64})

  # expected number of covariates
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting values for models
  m2s = model[1]*k
  m3s = m2s + model[2]*k
  y2s = model[1]*yppar*k
  y3s = y2s + model[2]*yppar*k

  @inbounds begin
    ### make fixed vectors of approximate probabilities for each state
    ## λ
    λpr = Array{Float64,1}[]
    for s in S
      push!(λpr, zeros(length(s.g) + div(length(vicsubsets(s.g)),2)))
    end
    # fill it up
    for si in Base.OneTo(ns)
      s = S[si]
      na = 0
      # within-area speciation
      if model[1]
        for a in s.g
          na += 1
          λpr[si][na] = expf(λ[(k+1)*s.h + a], β[s.h*ncov + a], md ? r[a] : r[1])*δt
        end
      else
        for a in s.g
          na += 1
          λpr[si][na] = λ[(k+1)*s.h + a]*δt
        end
      end
      # between-area speciation
      if na > 1
        λpr[si][na+1:end] .= λ[k+1 + (k+1)*s.h]*δt
      end
    end

    ## μ and loss
    μpr = Array{Float64,1}[]
    for s in S
      push!(μpr, zeros(length(s.g)))
    end
    # fill it up
    for si in Base.OneTo(ns)
      s = S[si]
      if isone(length(s.g))
        if model[2]
          for a in s.g
            μpr[si][1] = expf(μ[k*s.h + a], β[m2s + s.h*ncov + a], md ? r[y2s + a] : r[1])*δt
          end
        else
          for a in s.g
            μpr[si][1] = μ[k*s.h + a]*δt
          end
        end
      else
        na = 0
        for a in s.g
          na += 1
          μpr[si][na] = l[k*s.h + a]*δt
        end
      end
    end

    ## gain
    gpr = Array{Float64,1}[]
    for s in S
      push!(gpr, zeros(length(setdiff(as,s.g))))
    end
    # fill it up
    for si in Base.OneTo(ns)
      s = S[si]
      if model[3]
        for (i,ta) = enumerate(setdiff(as,s.g)), fa = s.g
          gpr[si][i] +=
            expf(q[s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                 β[m3s + s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                 md ? r[y3s + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)] : r[1])*δt
        end
      else
        for (i,ta) = enumerate(setdiff(as,s.g)), fa = s.g
          gpr[si][i] +=
            g[s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)]*δt
        end
      end
    end

    ## hidden states
    qpr = Array{Float64,1}[]
    for s in S
      push!(qpr, zeros(h-1))
    end
    # fill it up
    for si in Base.OneTo(ns)
      s = S[si]
       for (i, th) = enumerate(setdiff(hs, s.h))
        qpr[si][i] = q[s.h*(h-1) + (s.h > th ? s.h : s.h+1)]*δt
      end
    end

    # make λ, μ and g probability sums
    Sλpr = zeros(ns)
    Sμpr = zeros(ns)
    Sgpr = zeros(ns)
    Sqpr = zeros(ns)
    for i in Base.OneTo(ns)
      Sλpr[i] = sum(λpr[i])
      Sμpr[i] = sum(μpr[i])
      Sgpr[i] = sum(gpr[i])
      Sqpr[i] = sum(qpr[i])
    end
  end

  return λpr, μpr, gpr, qpr, Sλpr, Sμpr, Sgpr, Sqpr
end

"""
    makecorresschg(gpr::Array{Array{Float64,1},1},
                   μpr::Array{Array{Float64,1},1},
                   qpr::Array{Array{Float64,1},1},
                   λpr::Array{Array{Float64,1},1},
                   S  ::Array{Sgh, 1},
                   as::UnitRange{Int64},
                   hs::UnitRange{Int64},
                   k ::Int64,
                   ns::Int64)

Make vector of vectors of corresponding states changes
for **gains** and **looses** and **hidden states**
"""
function makecorresschg(gpr::Array{Array{Float64,1},1},
                        μpr::Array{Array{Float64,1},1},
                        qpr::Array{Array{Float64,1},1},
                        λpr::Array{Array{Float64,1},1},
                        S  ::Array{Sgh, 1},
                        as::UnitRange{Int64},
                        hs::UnitRange{Int64},
                        k ::Int64,
                        ns::Int64)
  # for *gains*
  gtos = Array{Int64,1}[]
  for i in Base.OneTo(ns)
    push!(gtos, zeros(Int64,length(gpr[i])))
  end
  for si in Base.OneTo(ns)
    s = S[si]
    for (i,ta) = enumerate(setdiff(as,s.g))
      gtos[si][i] = findfirst(r -> isSghequal(r,Sgh(union(s.g, ta),s.h)), S)
    end
  end

  # for *losses*
  μtos = Array{Int64,1}[]
  for i in Base.OneTo(ns)
    push!(μtos, zeros(Int64,length(μpr[i])))
  end
  for si in Base.OneTo(ns)
    s = S[si]
    isone(length(s.g)) && continue
    na = 0
    for a in s.g
      na += 1
      μtos[si][na] = findfirst(r -> isSghequal(r,Sgh(setdiff(s.g, a),s.h)), S)
    end
  end

  # for *hidden states*
  qtos = Array{Int64,1}[]
  for i in Base.OneTo(ns)
    push!(qtos, zeros(Int64,length(qpr[i])))
  end
  for si in Base.OneTo(ns)
    s = S[si]
     for (i, th) = enumerate(setdiff(hs, s.h))
      qtos[si][i] = th*(2^k-1) + si - s.h*(2^k-1)
    end
  end

  # for *hidden states*
  qtos = Array{Int64,1}[]
  for i in Base.OneTo(ns)
    push!(qtos, zeros(Int64,length(qpr[i])))
  end
  for si in Base.OneTo(ns)
    s = S[si]
     for (i, th) = enumerate(setdiff(hs, s.h))
      qtos[si][i] = th*(2^k-1) + si - s.h*(2^k-1)
    end
  end

  # for *speciation*
  λtos = Array{Array{Int64,1},1}[]
  for i in Base.OneTo(ns)
    push!(λtos, fill([0,0],length(λpr[i])))
  end
  for si in Base.OneTo(ns)
    s = S[si]
    na = 0
    if isone(length(s.g))
      for a in s.g
        na += 1
        λtos[si][na] = [si,si]
      end
    else
      # within-area speciation
      for a in s.g
        na += 1
        λtos[si][na] = [a + (2^k-1)*s.h,si]
      end
      # between-area speciation
      vs = vicsubsets(s.g)[1:div(end,2)]
      for v in vs
        na += 1
        λtos[si][na] =
          [findfirst(r -> isSghequal(r,Sgh(v[1],s.h)), S),
           findfirst(r -> isSghequal(r,Sgh(v[2],s.h)), S)]
      end
    end
  end

  return gtos, μtos, qtos, λtos
end

"""
    make_updλpr!(λ    ::Array{Float64,1},
                 β    ::Array{Float64,1},
                 λpr  ::Array{Array{Float64,1},1},
                 Sλpr ::Array{Float64,1},
                 δt   ::Float64,
                 k    ::Int64,
                 model::NTuple{3,Bool},
                 md   ::Bool,
                 ny   ::Int64)

Make λ instantaneous probability function for each state.
"""
function make_updλpr!(λ    ::Array{Float64,1},
                      β    ::Array{Float64,1},
                      λpr  ::Array{Array{Float64,1},1},
                      Sλpr ::Array{Float64,1},
                      δt   ::Float64,
                      k    ::Int64,
                      model::NTuple{3,Bool},
                      md   ::Bool,
                      ny   ::Int64)

  # expected number of covariates
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting values for models
  m2s = model[1]*k
  m3s = m2s + model[2]*k
  y2s = model[1]*yppar*k
  y3s = y2s + model[2]*yppar*k

  # if dependent on `z(t)` and multidimensional
  function f1(si ::Int64,
              s  ::Sgh,
              r::Array{Float64,1})

    @inbounds begin
      na = 0
      # within-area speciation
      for a in s.g
        na += 1
        λpr[si][na] = expf(λ[(k+1)*s.h + a], β[s.h*ncov + a], md ? r[a] : r[1])*δt
      end
      # between-area speciation
      if na > 1
        λpr[si][k+1:end] .= λ[k+1 + (k+1)*s.h]*δt
      end

      return Sλpr[si] = sum(λpr[si])
    end
  end

  # if dependent on `z(t)` and unidimensional
  function f2(si ::Int64,
              s  ::Sgh,
              r::Array{Float64,1})
    return Sλpr[si]
  end

  if model[1]
    return f1
  else
    return f2
  end
end

"""
    make_updμpr!(μ    ::Array{Float64,1},
                 l    ::Array{Float64,1},
                 β    ::Array{Float64,1},
                 μpr  ::Array{Array{Float64,1},1},
                 Sμpr ::Array{Float64,1},
                 δt   ::Float64,
                 k    ::Int64,
                 h    ::Int64,
                 model::NTuple{3,Bool},
                 md   ::Bool,
                 ny   ::Int64)

Make `μ` instantaneous probability function
"""
function make_updμpr!(μ    ::Array{Float64,1},
                      l    ::Array{Float64,1},
                      β    ::Array{Float64,1},
                      μpr  ::Array{Array{Float64,1},1},
                      Sμpr ::Array{Float64,1},
                      δt   ::Float64,
                      k    ::Int64,
                      h    ::Int64,
                      model::NTuple{3,Bool},
                      md   ::Bool,
                      ny   ::Int64)

  # expected number of covariates
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting values for models
  m2s = model[1]*k
  m3s = m2s + model[2]*k
  y2s = model[1]*yppar*k
  y3s = y2s + model[2]*yppar*k

  function f1(si::Int64,
              s ::Sgh,
              r ::Array{Float64,1})

    @inbounds begin
      if isone(length(s.g))
        for a in s.g
          μpr[si][1] =
            expf(μ[k*s.h + a], β[m2s + s.h*ncov + a], md ? r[y2s + a] : r[1])*δt
        end
      else
        for a in s.g
          μpr[si][a] =
            expf(μ[k*s.h + a], β[m2s + s.h*ncov + a], md ? r[y2s + a] : r[1])*δt
        end
      end

      return Sμpr[si] = sum(μpr[si])
    end
  end

  # if dependent on `z(t)` and unidimensional
  function f2(si ::Int64,
              s  ::Sgh,
              r::Array{Float64,1})
    return Sμpr[si]
  end

  if model[2]
    return f1
  else
    return f2
  end
end

"""
    make_updgpr!(g    ::Array{Float64,1},
                 β    ::Array{Float64,1},
                 gpr  ::Array{Array{Float64,1},1},
                 Sgpr ::Array{Float64,1},
                 δt   ::Float64,
                 k    ::Int64,
                 h    ::Int64,
                 model::NTuple{3,Bool},
                 md   ::Bool,
                 as   ::UnitRange{Int64},
                 S    ::Array{Sgh,1},
                 ny   ::Int64)

Make gain instantaneous probability function
"""
function make_updgpr!(g    ::Array{Float64,1},
                      β    ::Array{Float64,1},
                      gpr  ::Array{Array{Float64,1},1},
                      Sgpr ::Array{Float64,1},
                      δt   ::Float64,
                      k    ::Int64,
                      h    ::Int64,
                      model::NTuple{3,Bool},
                      md   ::Bool,
                      as   ::UnitRange{Int64},
                      S    ::Array{Sgh,1},
                      ny   ::Int64)

  # expected number of covariates
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting values for models
  m3s = model[1]*k + model[2]*k
  y3s = model[1]*yppar*k + model[2]*yppar*k

  # do setdiff before
  sdf = Array{Int64,1}[]
  for s = S
    push!(sdf, setdiff(as,s.g))
  end

  function f1(si ::Int64,
              s  ::Sgh,
              r::Array{Float64,1})
    @inbounds begin
      for (i,ta) = enumerate(sdf[si])
        gpr[si][i] = 0.0
        for fa = s.g
          gpr[si][i] +=
            expf(g[s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                 β[m3s + s.h*k*(k-1) + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)],
                 md ? r[y3s + (fa-1)*(k-1) + (ta > fa ? ta - 1 : ta)] : r[1])*δt
        end
      end
      return Sgpr[si] = sum(gpr[si])
    end
  end

  function f2(si ::Int64,
              s  ::Sgh,
              r::Array{Float64,1})
    return Sgpr[si]
  end

  if model[3]
    return f1
  else
    return f2
  end
end

"""
    make_updqpr!(Sqpr::Array{Float64,1})

Make hidden states instantaneous probability function
"""
function make_updqpr!(Sqpr::Array{Float64,1})

  function f(si ::Int64)
    return Sqpr[si]
  end

  return f
end

"""
    numberedges(ed::Array{Int64,2}, nspp::Int64)

Change numbering scheme so that tips are `1:numerofspecies` followed
by node numbering. MRCA is `numerofspecies+1`.
"""
function numberedges(ed::Array{Int64,2}, nspp::Int64)
  nel = 1
  nnd = nspp + 1

  edc         = zeros(Int64,size(ed))
  edc[1:2,1] .= nnd
  nnd        += 1

  for i in axes(ed,1)
    #if internal node
    if in(ed[i,2],ed[:,1])
      e1                = findfirst(isequal(ed[i,2]), ed[:,1])
      edc[e1:(e1+1),1] .= nnd
      edc[i,2]          = nnd
      nnd              += 1
    else
      edc[i,2] = nel
      nel     += 1
    end
  end

  return edc
end

"""
    expf(α::Float64, β::Float64, x::Float64)

Exponential regression function for rates with base rate `α`,
coefficient `β` and covariate `x`.
"""
expf(α::Float64, β::Float64, x::Float64) = α * exp(β * x)

"""
    prop_sample(s::Array{Float64,1}, prv::Array{Float64,1}, ns::Int64)

Sample one state with the respective (relative) probabilities
given `prv`.
"""
function prop_sample(s::Array{Float64,1}, prv::Array{Float64,1}, ns::Int64)

  cumsum!(s, prv)
  for i in Base.OneTo(ns)
    s[i] /= s[ns]
  end

  return searchsortedfirst(s, rand())
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# partial constraint development
using Base: Test
using IncrementalInference, Distributions
# going to introduce two new constraint types
import IncrementalInference: getSample

type DevelopPartial <: IncrementalInference.FunctorSingleton
  x::Distribution
  partial::Tuple
end
getSample(dpl::DevelopPartial, N::Int=1) = (rand(dpl.x, N)', )

type DevelopDim2 <: IncrementalInference.FunctorSingleton
  x::Distribution
end
getSample(dpl::DevelopDim2, N::Int=1) = (rand(dpl.x, N), )

N=50
fg = emptyFactorGraph()

v1 = addNode!(fg,:x1,randn(2,N),N=N)

pr = DevelopDim2(MvNormal([0.0;0.0],0.01*eye(2)))
f1  = addFactor!(fg,[:x1],pr)

dp = DevelopPartial(Normal(2.0, 1.0),(1,))
f2  = addFactor!(fg,[:x1],dp)

println("test evaluation of full constraint prior")
pts = evalFactor2(fg, f1, v1.index, N=N)
@test size(pts,1) == 2
@test size(pts,2) == N
@test norm(Base.mean(pts,2)[1]-[0.0]) < 0.3

println("test evaluation of partial constraint prior")
X1pts = getVal(v1)
pts = evalFactor2(fg, f2, v1.index, N=N)

@test size(pts, 1) == 2
@test size(pts,2) == N
@test norm(Base.mean(pts,2)[1]-[2.0]) < 0.75
# ensure the correct response from
@test norm(X1pts[1,:] - pts[1,:]) > 2.0
@test norm(X1pts[2,:] - pts[2,:]) < 1e-10
memcheck = getVal(v1)
@test norm(X1pts - memcheck) < 1e-10

tree = wipeBuildNewTree!(fg)

inferOverTreeR!(fg,tree, N=N)
pts = getVal(fg, :x1)
@test norm(Base.mean(pts,2)[1]-[0.0]) < 0.25
@test norm(Base.mean(pts,2)[2]-[0.0]) < 0.25

# plotKDE(getVertKDE(fg, :x1),levels=3)

type DevelopPartialPairwise <: IncrementalInference.FunctorPairwise
  x::Distribution
  partial::Tuple
  DevelopPartialPairwise(x::Distribution) = new(x, (2,))
end
getSample(dpl::DevelopPartialPairwise, N::Int=1) = (rand(dpl.x, N)', )

function (dp::DevelopPartialPairwise)(res::Vector{Float64},
                              idx::Int,
                              meas::Tuple, #{RowVector{Float64,Array{Float64,1}}}, #Tuple{Array{Float64,2}},
                              x1::Array{Float64},
                              x2::Array{Float64}  )
  #
  res[1] = meas[1][1,idx] - (x2[2,idx]-x1[2,idx])
  nothing
end

println("test evaluation of multiple simultaneous partial constraints")

v2 = addNode!(fg,:x2,100*randn(2,N),N=N)

valx2 = getVal(fg, :x2)

dpp = DevelopPartialPairwise(Normal(10.0, 1.0))
f3  = addFactor!(fg,[:x1;:x2],dpp)

dp2 = DevelopPartial( Normal(-20.0, 1.0), (1,) )
f4  = addFactor!(fg,[:x2],dp2)

pts = evalFactor2(fg, f3, v2.index, N=N)
# plotKDE(kde!(pts))
@test size(pts,1) == 2
@test norm(Base.mean(pts,2)[2]-[10.0]) < 3.0
@test norm(valx2[1,:] - pts[1,:]) < 1e-5

pts = evalFactor2(fg, f4, v2.index, N=N)
@test size(pts,1) == 2
@test norm(Base.mean(pts,2)[1]-[-20.0]) < 0.75
@test (Base.std(pts,2)[1]-1.0) < 0.4

# keep previous values to ensure funciton evaluation is modifying correct data fields

println("test findRelatedFromPotential...")

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl", 129)

X2pts = getVal(v2)
p3 = findRelatedFromPotential(fg, f3, v2.index, N)
@test Ndim(p3) == 2
pts = KernelDensityEstimate.getPoints(p3)
@test size(pts,2) == N

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2pts[1,:] - pts[1,:]) < 1e-10
@test norm(X2pts[2,:] - pts[2,:]) > 10*N*0.5
memcheck = getVal(v2)
@test norm(X2pts - memcheck) < 1e-10

X2pts = getVal(v2)
p4 = findRelatedFromPotential(fg, f4, v2.index, N)
@test Ndim(p4) == 2
pts = KernelDensityEstimate.getPoints(p3)
@test size(pts,2) == N

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2pts[1,:] - pts[1,:]) < 1e-10
@test norm(X2pts[2,:] - pts[2,:]) > 10*N*0.5
memcheck = getVal(v2)
@test norm(X2pts - memcheck) < 1e-10

println("test belief prediction...")
# partial prior
X2pts = getVal(v2)
val = predictbelief(fg, v2, [f4], N=N)
@test norm(X2pts[2,:] - val[2,:]) < 1e-10
@test 0.0 < norm(X2pts[1,:] - val[1,:])
@test norm(Base.mean(val[1,:])+20.0) < 0.75

# partial pairwise
X2pts = getVal(v2)
val = predictbelief(fg, v2, [f3], N=N)
@test norm(X2pts[1,:] - val[1,:]) < 1e-10
@test 0.0 < norm(X2pts[2,:] - val[2,:])
@test abs(Base.mean(val[2,:] - getVal(v1)[2,:])-10.0) < 0.75

# combination of partials
val = predictbelief(fg, v2, [f3;f4], N=N)
# plotKDE(kde!(val),levels=3)
@test norm(Base.mean(val,2)[1]-[-20.0]) < 2.0
@test norm(Base.mean(val,2)[2]-[10.0]) < 2.0
@test (Base.std(val,2)[1]-1.0) < 3.0
@test (Base.std(val,2)[2]-1.0) < 3.0

tree = wipeBuildNewTree!(fg )#, drawpdf=true)
# run(`evince bt.pdf`)

inferOverTreeR!(fg,tree, N=N)

pts = getVal(fg, :x1)
@test norm(Base.mean(pts,2)[1]-[0.0]) < 0.5
@test norm(Base.mean(pts,2)[2]-[0.0]) < 0.5

pts = getVal(fg, :x2)
@test norm(Base.mean(pts,2)[1]-[-20.0]) < 2.0
@test norm(Base.mean(pts,2)[2]-[10.0]) < 2.0
@test (Base.std(pts,2)[1]-1.0) < 3.0
@test (Base.std(pts,2)[2]-1.0) < 3.0

# plotKDE(getVertKDE(fg, :x2),levels=3)

# spyCliqMat(tree, :x2)

#


text:
using Test

using RealNeuralNetworks.Neurons
using RealNeuralNetworks.NBLASTs
using RealNeuralNetworks.Utils.VectorClouds
using RealNeuralNetworks.Utils.RangeIndexingArrays
using RealNeuralNetworks.Neurons.Segments

using CSV

ASSET_DIR = joinpath(@__DIR__, "../asset/")

NEURON_ID1 = 77625
NEURON_ID2 = 77641

const k = 20

@testset "test NBLAST module..." begin
    neuron1 = Neurons.load(joinpath(ASSET_DIR, "$(NEURON_ID1).swc.bin"));
    neuron2 = Neurons.load(joinpath(ASSET_DIR, "$(NEURON_ID2).swc.bin"));

    vectorCloud1 = NBLASTs.VectorCloud(neuron1)
    vectorCloud2 = NBLASTs.VectorCloud(neuron2)

    # transform to micron
    vectorCloud1[1:3, :] ./= 1000
    vectorCloud2[1:3, :] ./= 1000

    # read the precomputed score matrix, which is the joint distribution of the scores
    println("\nread the precomputed score table...")
    SMAT_FCWB_PATH = joinpath(ASSET_DIR, "smat_fcwb.csv")
    df = CSV.read(SMAT_FCWB_PATH)
    ria = RangeIndexingArray{Float32}(df)

    println("\ncompute nblast score...")
    @time score = NBLASTs.nblast(vectorCloud2, vectorCloud1; ria=ria)
    @time score = NBLASTs.nblast(vectorCloud2, vectorCloud1; ria=ria)
    println("query $(NEURON_ID1) against target $(NEURON_ID2): ", score)
    @test isapprox(score, Float32(50887.82f0))

    println("\ncompute nblast score...")
    @time score = NBLASTs.nblast(vectorCloud1, vectorCloud2; ria=ria)
    println("query $(NEURON_ID1) against target $(NEURON_ID2): ", score)
    @test isapprox(score, Float32(53697.95f0))

    #@time score = NBLASTs.nblast(vectorCloud1, vectorCloud1; ria=ria)
    #println("query $(NEURON_ID1) against itself: ", score)
    #@test isapprox(score, Float32(86507.24f0))

    neuronList = [neuron1, neuron2]
    # the result from R NBLAST is :
    # ID    77625	    77641
    # 77625	86501.20 	53696.72
    # 77641	50891.03 	101011.08
    RNBLAST_RESULT = Float32[86501.20 	53696.72; 50891.03 	101011.08]
    println("\ncompute similarity matrix...")
    @time rawSimilarityMatrix, normalizedSimilarityMatrix, meanSimilarityMatrix = NBLASTs.nblast_allbyall(neuronList; ria=ria, k=k,
                                                        downscaleFactor=1000)
    @show rawSimilarityMatrix
    @test isapprox.(rawSimilarityMatrix,  RNBLAST_RESULT) |> all

    #println("\ncompute normalised similarity matrix...")
    #RNBLAST_RESULT = Float32[1.0 0.5315924; 0.5883275 1.0]
    #@show normalizedSimilarityMatrix
    #@test isapprox.(normalizedSimilarityMatrix,  RNBLAST_RESULT) |> all
    #
    #RNBLAST_RESULT = Float32[1.0 0.5599599; 0.5599599 1.0]
    #println("\ncompute mean similarity matrix...")
    #@show meanSimilarityMatrix
    #@test isapprox.(meanSimilarityMatrix,  RNBLAST_RESULT) |> all
    #
    #vectorCloudList = [vectorCloud1, vectorCloud2]
    #RNBLAST_RESULT = Float32[86501.20 	53696.72; 50891.03 	101011.08]
    #@time rawSimilarityMatrix = NBLASTs.nblast_allbyall(vectorCloudList; ria=ria)
    #@show rawSimilarityMatrix
    #@test isapprox.(rawSimilarityMatrix,  RNBLAST_RESULT) |> all

    #
    #println("\ncompute nblast score using default zfish score table ...")
    ## this test will not work since neuron reading from swc/swc.bin
    ## do not support node class yet
    #vectorCloud1 = NBLASTs.VectorCloud(neuron1; class=Segments.DENDRITE_CLASS)
    #vectorCloud2 = NBLASTs.VectorCloud(neuron2; class=Segments.DENDRITE_CLASS)
    #@time score = NBLASTs.nblast(vectorCloud2, vectorCloud1;)
    #println("query $(NEURON_ID1) against target $(NEURON_ID2): ", score)
    ## @test isapprox(score, -63722.48f0)
    #
    #println("\ntest similarity measurement for dendrites...")
    #@time rawSimilarityMatrix, normalizedSimilarityMatrix, meanSimilarityMatrix =
    #            NBLASTs.nblast_allbyall(neuronList; ria=ria, k=k, downscaleFactor=1000)
    #@show rawSimilarityMatrix
    ##println("similarity matrix of dendrites: ", similarityMatrix)

    #println("compute similarity score with semantic small-to-big...")
    #semanticSmall2bigSimilarityMatrix =
    #            NBLASTs.nblast_allbyall_small2big(neuronList; k=k, ria=ria, semantic=true);

end # end of test module


text:
@testset "Hybrid System tests" begin
    test_sys = PSB.build_system(PSB.PSITestSystems, "c_sys14"; add_forecasts = false)

    h_sys = HybridSystem(
        name = "Test H",
        available = true,
        status = true,
        bus = get_component(Bus, test_sys, "Bus 1"),
        active_power = 1.0,
        reactive_power = 1.0,
        thermal_unit = ThermalStandard(nothing),
        electric_load = PowerLoad(nothing),
        storage = GenericBattery(nothing),
        renewable_unit = RenewableDispatch(nothing),
        base_power = 100.0,
        operation_cost = TwoPartCost(nothing),
    )
    add_component!(test_sys, h_sys)

    initial_time = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)
    other_time = initial_time + resolution
    name = "test"
    horizon = 24
    data = Dict(initial_time => ones(horizon), other_time => 5.0 * ones(horizon))
    forecast = Deterministic(name, data, resolution)
    add_time_series!(test_sys, h_sys, forecast)
    ts = get_time_series(Deterministic, h_sys, "test")
    @test isa(ts, Deterministic)
    @test_throws ErrorException to_json(test_sys, "./test_sys.json", force = true)
end


text:
export Neighbourhood, LocalNeighbourhood, GlobalNeighbourhood, HierachicalNeighbourhood

abstract type Neighbourhood{I<:Integer} <: AbstractVector{I} end
# Base.eltype(::Type{<:Neighbourhood{I}}) = I
# Base.eltype(n::Neighbourhood{I}) = I
Base.IndexStyle(::Type{<:Neighbourhood}) = IndexLinear()
Base.size(n::Neighbourhood) = (n.particle_number,)
Base.length(n::Neighbourhood) = n.particle_number
Base.getindex(n::Neighbourhood, i::Integer) = getneighbour(n, i)
Base.setindex!(n::Neighbourhood, v, i::Integer) = error("Not supported!")
rearrange!(n::Neighbourhood, results_best::AbstractVector{<:Number}, compare::Function) = nothing

"""
    LocalNeighbourhood

The LocalNeighbourhood organizes the particle in a ring topology. For width=1 for example,
each particle has three neighbours, one to the left, one to the right and itself.
The width says to how many particles on the left and on the right it is connected.
"""
struct LocalNeighbourhood{I} <: Neighbourhood{I}
    particle_number::I
    neighbours::AbstractVector{<:AbstractVector{<:I}}
end
function LocalNeighbourhood(particle_number::I, width::I=one(I)) where I <: Integer
    @assert (particle_number ≥ 2*width + 1) "For $particle_number particles the width can be $(Int(round((particle_number-1)/2))) at most."
    all_neigs = Vector{Vector{I}}(undef, particle_number)
    for current in one(particle_number):particle_number
        neighbours = collect((current-width):(current+width))
        neighbours[neighbours.<one(I)] .+= particle_number
        neighbours[neighbours.>particle_number] .-= particle_number
        all_neigs[current] = neighbours
    end
    LocalNeighbourhood(particle_number, all_neigs)
end
getneighbour(l::LocalNeighbourhood, i::Integer) = l.neighbours[i]

"""
    GlobalNeighbourhood

In the GlobalNeighbourhood all particle are connected to all other particles.
"""
struct GlobalNeighbourhood{I} <: Neighbourhood{I}
    particle_number::I
    neighbours::AbstractVector{<:I}
end
function GlobalNeighbourhood(particle_number::Integer)
    GlobalNeighbourhood(particle_number, collect(one(particle_number):particle_number))
end
getneighbour(g::GlobalNeighbourhood, i::Integer) = g.neighbours

"""
    HierachicalNeighbourhood

Inspired by this paper: https://ieeexplore.ieee.org/abstract/document/1299745/
It arranges all particles in a tree. The parent of a node is the only neighbour. The root has itself as neighbours.
After every iteration, the tree will be updated. Particles can climbe up and down in the tree.

# Initialization
HierachicalNeighbourhood(particle_number::Integer, branching_degree::Integer)
HierachicalNeighbourhoodByTreeHeight(tree_height::Integer, branching_degree::Integer)
        - The particle number results from the full tree.
"""
struct HierachicalNeighbourhood{I} <: Neighbourhood{I}
    single_particle::AbstractVector{<:AbstractVector{<:I}}
    particle_number::I
    branching_degree::I
    childs::AbstractVector{<:AbstractVector{<:I}}
    parent::AbstractVector{<:I}
    particle_to_graph::AbstractVector{<:I}
    graph_to_particle::AbstractVector{<:I}
end
function HierachicalNeighbourhood(particle_number::I, branching_degree::I) where I<:Integer
    # @assert particle_number ≤ max_particle_number "With d=$branching_degree and h=$tree_height you can use $max_particle_number particles at most."
    childs = getchilds(particle_number, branching_degree)
    parent = getparents(childs)
    particle_to_graph = collect(one(I):particle_number)
    graph_to_particle = collect(one(I):particle_number)
    HierachicalNeighbourhood([[i] for i in one(I):particle_number], particle_number, branching_degree, childs, parent, particle_to_graph, graph_to_particle)
end
function HierachicalNeighbourhoodByTreeHeight(tree_height::I, branching_degree::I) where I<:Integer
    particle_number = convert(Int, (branching_degree^tree_height - one(I))/(branching_degree - one(I)))
    HierachicalNeighbourhood(particle_number, branching_degree)
end
getneighbour(g::HierachicalNeighbourhood, i::Integer) = g.single_particle[g.graph_to_particle[g.parent[g.particle_to_graph[i]]]]

"""
    getchilds(particle_number::T, branching_degree::T) where T <: Integer

Get a list of all childs for each particle number.
"""
function getchilds(particle_number::I, branching_degree::I) where I <: Integer
    childs = [I[] for _ in one(I):particle_number]
    last_childs = [one(I)]
    cummmulative_width = one(I)
    current_height = one(I)
    while true
        width = branching_degree^current_height
        for current_branch in one(I):width
            child_index = cummmulative_width + current_branch
            child_index > particle_number && return childs
            push!(childs[last_childs[(current_branch - one(I)) % length(last_childs) + one(I)]], child_index)
        end
        last_childs = vcat(childs[last_childs]...)
        cummmulative_width += width
        current_height += one(I)
    end
end

"""
    getparents(childs::AbstractVector{<:AbstractVector{<:Integer}})

Get a list of parents from a list of cilds.
"""
function getparents(childs::AbstractVector{<:AbstractVector{<:I}}) where I<:Integer
    parents = ones(eltype(childs[one(I)]), length(childs))
    for parent_index in (one(I)+one(I)):length(childs)
        isempty(childs[parent_index]) && continue
        parents[childs[parent_index]] .= parent_index
    end
    return parents
end

"""
    function bestchild(graph_to_particle::AbstractVector{I}, childs_graph::AbstractVector{I}, personal_best::AbstractVector{<:Number}, compare::Function) where I<:Integer

Resturns the "best" child of a node in the neighbourhood tree. What the best is, is determined by compare function. If more than one particle have the "best" personal best value, the first in the row is taken.
"""
@inline function bestchild(graph_to_particle::AbstractVector{I}, childs_graph::AbstractVector{I}, personal_best::AbstractVector{<:Number}, compare::Function) where I<:Integer
    @inbounds child_min_particle = graph_to_particle[childs_graph[one(I)]]
    @inbounds minval = personal_best[child_min_particle]
    @inbounds for child in childs_graph
        childs_particle = graph_to_particle[child]
        nextval = personal_best[childs_particle]
        if compare(nextval, minval)
            minval = nextval
            child_min_particle = childs_particle
        end
    end
    return child_min_particle
end

"""
    function swapifnecessary!(n::HierachicalNeighbourhood{I}, personal_best::AbstractVector{<:Number}, child_min_particle::I, parent_particle::I, child_min_graph::I, parent_graph::I, compare::Function) where I<:Integer

Swaps to particles in the neighbourhood tree if the personal_best of the particle in the subtree is "better" than the node. "better" means if the `compare` function results in true.
"""
@inline function swapifnecessary!(n::HierachicalNeighbourhood{I}, personal_best::AbstractVector{<:Number}, child_min_particle::I, parent_particle::I, child_min_graph::I, parent_graph::I, compare::Function) where I<:Integer
    @inbounds if compare(personal_best[child_min_particle], personal_best[parent_particle])
        swap!(n.particle_to_graph, child_min_particle, parent_particle)
        swap!(n.graph_to_particle, child_min_graph, parent_graph)
    end
    return
end

"""
    rearrange!(n::HierachicalNeighbourhood, personal_best::AbstractVector{<:Number}, compare::Function)

If a parent node has a worse (higher or lesser) evaluation than the best child, they will be swapped.
The tree will be checked from the top to the bottom. A particle can at most climbe one level up,
but can travel down from the top to the bottom in one function call.
"""
function rearrange!(n::HierachicalNeighbourhood, personal_best::AbstractVector{<:Number}, compare::Function)
    n_childs::Int = n.particle_number
    @inbounds for parent_graph in 1:n_childs
        childs_graph = n.childs[parent_graph]
        isempty(childs_graph) && continue
        child_min_particle = bestchild(n.graph_to_particle, childs_graph, personal_best, compare)
        child_min_graph = n.particle_to_graph[child_min_particle]
        parent_particle = n.graph_to_particle[parent_graph]
        swapifnecessary!(n, personal_best, child_min_particle, parent_particle, child_min_graph, parent_graph, compare)
    end
    return
end


text:
export hiddenLayer

type hiddenLayer <:Layer
  bottom ::Array{Layer}
  top    ::Array{Layer}
  topRange ::Vector

  rᵢ     ::AFArray{Float32} #input rate
  W      ::AFArray{Float32}
  ∂W     ::AFArray{Float32}
  b      ::AFArray{Float32}
  ∂b     ::AFArray{Float32}
  I      ::AFArray{Float32}
  A      ::Function # activation function
  ∂A     ::Function # derivative of activation function
  rₒ     ::AFArray{Float32} # output rate

  ∂rᵢ_L  ::AFArray{Float32}
  ∂w_L   ::AFArray{Float32}
  ∂b_L   ::AFArray{Float32}
  ∂i_L   ::AFArray{Float32}
  ∂rₒ_L  ::AFArray{Float32}

  lock   ::Bool
end
hiddenLayer(inDim::Int,outDim::Int)=
  hiddenLayer(Layer[],
              Layer[],
              Vector{AFArray}(),

              rand(AFArray{Float32},inDim,batchsize),
              0.01*randn(AFArray{Float32},outDim,inDim),
              zeros(AFArray{Float32},outDim,inDim),
              zeros(AFArray{Float32},outDim,1),
              zeros(AFArray{Float32},outDim,1),
              randn(AFArray{Float32},outDim,batchsize),
              recLinear,#sigmoid
              ∂recLinear,#∂sigmoid
              rand(AFArray{Float32},outDim,batchsize),

              randn(AFArray{Float32},inDim,batchsize),
              randn(AFArray{Float32},outDim,inDim,batchsize),
              randn(AFArray{Float32},outDim,batchsize),
              randn(AFArray{Float32},outDim,batchsize),
              randn(AFArray{Float32},outDim,batchsize),

              false)


text:
using Symbolics, Test, DiffRules, SpecialFunctions

@sym x y z m ν

@testset "             Derivatives" begin
    f(x) = x^3
    g(x) = x^2
    @test D(f+g)(x) |> simplify == 3 * x ^ 2 + 2x
    @test (D^2)(f+g)(x) == 3 * (2x) + 2
    @test (D^3)(f+g)(x) == 6
    @test D(x^2 + y^2 + z^2, x) == 2x
    @test D(D(2x*y, x), y) == 2
    # test diff rules
    @test D(x -> atan(x,y))(x) ==  y/(x^2 + y^2)
    @test D(y -> atan(x,y))(y) == -x/(x^2 + y^2)
    @test D(x -> hypot(x,y))(x) == x/hypot(x,y)
    @test D(y -> hypot(x,y))(y) == y/hypot(x,y)

    @test D(x -> besselj(ν,x))(x) == (besselj(ν - 1, x) - besselj(ν + 1, x))/2
    @test D(x -> besseli(ν,x))(x) == (besseli(ν - 1, x) + besseli(ν + 1, x))/2
    @test D(x -> bessely(ν,x))(x) == (bessely(ν - 1, x) - bessely(ν + 1, x))/2
    @test D(x -> besselk(ν,x))(x) == (besselk(ν - 1, x) + besselk(ν + 1, x))/2

    @test D(x -> hankelh1(ν,x))(x) == (hankelh1(ν - 1, x) - hankelh1(ν + 1, x))/2
    @test D(x -> hankelh2(ν,x))(x) == (hankelh2(ν - 1, x) - hankelh2(ν + 1, x))/2

    @test D(x -> polygamma(m,x))(x) == polygamma(m + 1, x)

    @test D(x -> beta(x,y))(x) == beta(x,y)*(digamma(x)-digamma(x+y))
    @test D(y -> beta(x,y))(y) == beta(x,y)*(digamma(y)-digamma(x+y))

    @test D(x -> lbeta(x,y))(x) == digamma(x)-digamma(x+y)
    @test D(y -> lbeta(x,y))(y) == digamma(y)-digamma(x+y)

    for (M, f, arity) in DiffRules.diffrules()
        if arity == 1 && (M == :Base || M == :SpecialFunctions) && f ∉ [:inv, :+, :-, :abs, ] # [:bessely0, :besselj0, :bessely1, :besselj1]
            deriv = DiffRules.diffrule(M, f, :x)
            @eval begin
                @test D($M.$f)(x) == $deriv
            end
        end
    end
end


text:
@testset "comms_t" begin
    utime = 123456
    format = "format_name_here"
    format_version_major = 1
    format_version_minor = 32
    data = rand(UInt8, 1000)
    msg = DrakeVisualizer.Comms.CommsT(
        utime,
        format,
        format_version_major,
        format_version_minor,
        data)
    encoded = DrakeVisualizer.Comms.encode(msg)
    decoded = DrakeVisualizer.Comms.decode(encoded, DrakeVisualizer.Comms.CommsT)

    for field in fieldnames(typeof(msg))
        @test getfield(msg, field) == getfield(decoded, field)
    end
end


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

include("../perfutil.jl")

include("level1.jl")
include("level2.jl")
include("level3.jl")

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
function commander(c::Client, m::Message, ::Val{:julia_con})
    command = extract_command("jc", m.content)
    args = split(command)
    @debug "parse result" command args

    if length(args) == 0 ||
        args[1] ∉ ["2021", "now", "today", "tomorrow", "day"] # next is not available yet
        help_commander(c, m, :julia_con)
        return
    end

    arg = (occursin(r"[0-9]+", args[1]) ? "jc" : "") * args[1]
    jc_execute(c, m, Symbol(arg), args[2:end])

end

jc_execute(c, m, arg::Symbol, args) = jc_execute(c, m, Val(arg), args)

function jc_execute(c, m, ::Val{:jc2021}, args)
    @info "2021 was required" m.content m.author.username m.author.discriminator
    jc_juliacon2021(c, m)
end

function jc_execute(c, m, ::Val{:now}, args)
    @info "now was required" m.content m.author.username m.author.discriminator args

    tz_arg = isempty(args) || !TimeZones.istimezone(args[1]) ? "UTC" : args[1]
    current = ZonedDateTime(now(), TimeZone(tz_arg))
    # current = ZonedDateTime(Dates.DateTime("2021-07-30T21:30:00.000"), TimeZone(tz_arg))
    @info "now command" tz_arg current

    if isempty(JuliaCon.get_running_talks(now = current))
        not_now = "There is no JuliaCon program now."
        @info not_now
        today = " Try `jc today`"
        schedule = "(Full schedule: https://live.juliacon.org/agenda)"
        reply(c, m, not_now * today * "\n" * schedule)
    else
        reply(c, m, JuliaCon.now(now = current, output=:text))
    end
    return nothing
end

function jc_execute(c, m, ::Val{:today}, args)
    tz_arg = isempty(args) || !TimeZones.istimezone(args[1]) ? "UTC" : args[1]
    dt = now(TimeZone(tz_arg))
    jc_today(c, m, dt)
    return nothing
end

function jc_execute(c, m, ::Val{:tomorrow}, args)
    tz_arg = isempty(args) || !TimeZones.istimezone(args[1]) ? "UTC" : args[1]
    dt = now(TimeZone(tz_arg)) + Day(1)
    jc_today(c, m, dt)
    return nothing
end

function jc_execute(c, m, ::Val{:next}, args)
    @info "Calling :next"
end

function jc_execute(c, m, ::Val{:day}, args)
    zonedDT = ZonedDateTime(Date(args[1]), tz"UTC")
    jc_today(c, m, zonedDT)
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:julia_con})
    reply(c, m, """
        How to look for information about JuliaCon 2021 with the `jc` command:
        ```
        jc help
        jc 2021
        jc now
        jc today [timezone]
        jc tomorrow [timezone]
        jc day <xxxx-mm-dd>
        ```
        `help` prints this message
        `2021` greets the users and provide a link to the official website of year 2021's edition
        `now` details the talks and events occuring right now on the different tracks
        `today`, `tomorrow`, `day` lists the talks and event of the day. [timezone] is an optional argument which need to be valid following the TimeZones.jl package

        Examples
        ```
        jc today Asia/Tokyo
        jc tomorrow America/Juneau
        jc day 2021-07-27
        ```
        """)
end

function jc_juliacon2021(c::Client, m::Message)
    reply(c, m, "Welcome to JuliaCon 2021! Find more information at https://juliacon.org/2021/.")
    return nothing
end

function split_pretty_table(str)
    lines = split(str, "\n")
    acc = "```\n"
    strings = Vector{String}()
    for line in lines
        if sum(length, [acc, line]) ≤ 1996
            acc *= "\n" * line
        else
            push!(strings, acc * "\n```")
            acc = "```\n" * line
        end
    end
    push!(strings, acc * "\n```")
    return strings
end

function jc_today(c::Client, m::Message, zonedDT::ZonedDateTime)
    responses = jc_today(zonedDT)
    foreach(s -> reply(c, m, s), responses)
end

function jc_today(zonedDT::ZonedDateTime)
    @info "jc_today" zonedDT
    aux = JuliaCon.today(now = zonedDT, output = :text)
    if aux === nothing
        not_today = "There is no JuliaCon program at $(zonedDT)"
        schedule = "(Full schedule: https://live.juliacon.org/agenda)"
        return [not_today * "\n" * schedule]
    else
        strings = Vector{String}()
        tables, legend = aux[1:end-1], aux[end]
        for t in tables, str in split_pretty_table(t)
            push!(strings, str)
        end
        push!(strings, legend)
        return strings
    end
end