text:
#=
read_rdb_hdr.jl
based on
https://gitlab.com/fMRI/toppe/blob/master/+toppe/+utils/read_rdb_hdr.m

Matlab version is:
Copyright (c) 2012 by General Electric Company. All rights reserved.

2019-05-22 Julia version by Jeff Fessler
=#

export read_rdb_hdr

# using MIRTio: header_*, rdb_hdr_26_002_def

"""
    ht = read_rdb_hdr(fid::IOStream)

Read GE raw (RDB) header for MRI scans

Returns NamedTuple `ht` with header values accessible by `ht.key`
"""
function read_rdb_hdr(fid::IOStream)
    seek(fid, 0)
    rdbm_rev = read(fid, Float32) # raw header (RDBM) revision number

    if rdbm_rev == Float32(26.002)
        hd = rdb_hdr_26_002_def()
        ht = header_read(fid, hd)
        ht = header_string(ht) # convert Array{Cuchar} to String
    else
        throw("unknown RDBM rev $rdbm_rev")
    end

    (ht.rdbm_rev != rdbm_rev) && throw("rev mismatch: $(ht.rdbm_rev) != $rdbm_rev")

    return ht
end

"""
    ht = read_rdb_hdr(file::String)
Read from `file`
"""
function read_rdb_hdr(file::String)
    return open(read_rdb_hdr, file)
end


text:
function triangle_centrality1(A)
    T = mul(A, A', Semirings.PLUS_TIMES[Float64], mask=A)
    y = reduce(Monoids.PLUS_MONOID[Float64], T, dims=2)
    k = reduce(Monoids.PLUS_MONOID[Float64], y)
    return (3 * mul(A, y) - 2 * mul(T, y) .+ y) ./ k
end

function triangle_centrality2(A)
    T = mul(A, A', Semirings.PLUS_PAIR, mask=A, desc=Descriptors.S)
    y = reduce(Monoids.PLUS_MONOID[Float64], T, dims=2)
    k = reduce(Monoids.PLUS_MONOID[Float64], y)
    return map(BinaryOps.TIMES, eadd(eadd(3.0 * mul(A, y), 2.0 * mul(T, y, Semirings.PLUS_SECOND), BinaryOps.MINUS), y), 1/k)
end

ext: jl


text:
"""
    Catwalk.jl

An optimizing Just In Time compiler written in Julia.
"""
module Catwalk

export @jit, JIT

include("typelist.jl")
include("frequencies.jl")
include("compile.jl")
include("profile.jl")
include("costmodel.jl")
include("optimize.jl")
include("explore.jl")

mutable struct CallBoost
    fnsym::Symbol
    profilestrategy::ProfileStrategy
    optimizer::Optimizer
    currentprofiler::Profiler
    round::Int
    CallBoost(fnsym; profilestrategy=SparseProfile(), optimizer=TopNOptimizer()) = new(fnsym, profilestrategy, optimizer, NoProfiler(), 1)
end

function step!(boost::CallBoost)
    boost.round += 1
end

function ctx(boost::CallBoost)
    newprofiler = getprofiler(boost.profilestrategy, boost.round)
    context = CallCtx{typeof(newprofiler), fixtypes(boost.optimizer, boost.currentprofiler)}(newprofiler)
    empty!(newprofiler)
    boost.currentprofiler = newprofiler
    return context
end

struct CallCtx{TProfiler, TFixtypes}
    profiler::TProfiler
    CallCtx() = new{NoProfiler, EmptyTypeList}()
    CallCtx{TProfiler, TFixtypes}(profiler) where {TProfiler, TFixtypes} = new{TProfiler, TFixtypes}(profiler)
end

fixtypes(::Type{CallCtx{TProfiler, TFixtypes}}) where {TProfiler, TFixtypes} = TFixtypes
profiler(ctx::CallCtx) = ctx.profiler

struct JIT
    id::Int
    callboosts::Vector{CallBoost}
    explorer::Explorer
end
function JIT(callboosts...; explorerfactory=BasicExplorer)
    id = rand(Int)
    jit = JIT(id, [], explorerfactory(id))
    for boost in callboosts
        add_boost!(jit, boost)
    end
    return jit
end

optimizerid(jit::JIT) = jit.id

function add_boost!(jit::JIT, boost)
    push!(jit.callboosts, boost)
    register_callsite!(jit.explorer, boost.fnsym)
end

function step!(jit::JIT)
    update_callboosts(jit)
    step!.(jit.callboosts)
    step!(jit.explorer)
end

function update_callboosts(jit::JIT)
    currentsyms = Set(map(b -> b.fnsym, jit.callboosts))
    for newsite in setdiff(get_freshcallsites!(jit.id), currentsyms)
        add_boost!(jit, CallBoost(newsite))
    end
end

struct OptimizerCtx{TCallCtxs, TExplorer}
    callctxs::TCallCtxs
    explorer::TExplorer
    OptimizerCtx() = new{NamedTuple{}, NoExplorer}(NamedTuple(), NoExplorer())
    OptimizerCtx(optimizerid, callctxs, explorer) = new{typeof(callctxs), typeof(explorer)}(callctxs, explorer)
end

function ctx(jit::JIT)
    callctxs = (;map(boost -> (boost.fnsym, ctx(boost)), jit.callboosts)...)
    return OptimizerCtx(optimizerid(jit), callctxs, jit.explorer)
end

function callctx(::Type{OptimizerCtx{TCallCtxs, TExplorer}}, key) where {TCallCtxs, TExplorer}
    return callctx(TCallCtxs, key)
end

function callctx(ctxs::Type{NamedTuple{TNames, TVals}}, key) where {TNames, TVals}
    name = nameof(key)
    idx = findfirst(n -> n == name, TNames)
    if isnothing(idx)
        error("$name (type: $(typeof(name))) not found in $ctxs")
    end
    return fieldtypes(TVals)[idx]
end

function explorer(::Type{
        OptimizerCtx{
            TCallCtxs,
            TExplorer
        }
    }) where {TCallCtxs, TExplorer}
    return TExplorer
end

end # module


text:
#=
Copyright (c) 2015, Roberto Garrappa
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
    * Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in
      the documentation and/or other materials provided with the distribution
    * Neither the name of the Department of Mathematics - University of Bari - Italy nor the names
      of its contributors may be used to endorse or promote products derived
      from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
=#

"""
@article{Garrappa2015NumericalEO,
  title={Numerical Evaluation of Two and Three Parameter Mittag-Leffler Functions},
  author={Roberto Garrappa},
  journal={SIAM J. Numer. Anal.},
  year={2015},
  volume={53},
  pages={1350-1369}
}
"""

myeps(x) = x |> one |> float |> eps
myeps(x::Complex) =  x |> real |> myeps

"""
    mittleffderiv(α, β, z)

Compute the derivative of the Mittag-Leffler function at `z` for parameters `α, β`.
"""
function mittleffderiv(α, β, z)
    q = 1//2
    #case 1, small z
    if abs(z) <= q
        ω = α + β - 3//2
        D = α^2 - 4*α*β + 6*α + 1
        #k1
        if α>1
            k₁ = ((2 - α - β)/(α - 1)) + 1
        elseif α > 0 && α <= 1 && D <= 0
            k₁ = ((3 - α - β)/α) + 1
        else
            k₁ = maximum([((3-α-β)/α) + 1, ((1-2*ω*α + sqrt(D))/(2*(α^2)))+1])
        end
        k₀ = maximum([k₁, log(myeps(z)*(1-abs(z)))/log(abs(z))])
        k₀ = ceil(Int, k₀)
        out = zero(z)
        for k in 0:k₀
            out = out + ((k + 1)*z^k)/(gamma(α + β + α*k))
        end
    #case 2, larger z
    else
        out = (mittleff(α, β-1, z) - (β - 1)*mittleff(α, β, z))/(α*z)
    end
    return out
end

"""
    mittleffderiv(α, z)

Compute mittleffderiv(α, 1, z)
"""
mittleffderiv(α, z) = mittleffderiv(α, 1, z)

"""
    mittleff(α, β, γ, z)

Compute three-parametric mittleff(α, β, γ, z).

```tex
@article{2015,
   title={Numerical Evaluation of Two and Three Parameter Mittag-Leffler Functions},
   volume={53},
   ISSN={1095-7170},
   url={http://dx.doi.org/10.1137/140971191},
   DOI={10.1137/140971191},
   number={3},
   journal={SIAM Journal on Numerical Analysis},
   publisher={Society for Industrial & Applied Mathematics (SIAM)},
   author={Garrappa, Roberto},
   year={2015},
   month={Jan},
   pages={1350–1369}
}
```
"""
function mittleff(α, β, γ, z)
    log_epsilon = log(10^(-15))
    E = zeros(Float64)
    abs(z) < 1e-15 ? E=1/gamma(β) : E=LTInversion(1, z, α, β, γ, log_epsilon)
    return E
end

"""
    mittleff(α, β, z)

Compute two-parametric Mittag Leffler function
"""
mittleff(α, β, z) = mittleff(α, β, 1, z)

"""
    mittleff(α, z)

Compute single-parametric Mittag Leffler function
"""
mittleff(α, z) = mittleff(α, 1, 1, z)

mittleff(α, vec::AbstractArray) = map(x -> mittleff(α, 1, x), vec)
mittleff(α, β, vec::AbstractArray) = map(x -> mittleff(α, β, x), vec)
mittleff(α, β, γ, vec::AbstractArray) = map(x -> mittleff(α, β, γ, x), vec)

function LTInversion(t, lambda, alpha, beta, gama, log_epsilon)
    theta = angle(lambda)
    kmin = ceil(Int, -alpha/2 - theta/2/pi)
    kmax = floor(Int, alpha/2 - theta/2/pi)
    k_vett = collect(kmin:kmax)
    s_star = abs(lambda)^(1/alpha) * exp.(im*(theta.+2*k_vett*pi)/alpha)

    phi_s_star = (real(s_star).+ abs.(s_star))/2
    index_s_star = sortperm(phi_s_star)
    phi_s_star = sort(phi_s_star)
    s_star = s_star[index_s_star]
    index_save = phi_s_star .> 1e-15

    s_star = s_star[index_save]
    phi_s_star = phi_s_star[index_save]
    s_star = [0; s_star]
    phi_s_star = [0; phi_s_star]
    J1 = length(s_star)
    J = J1 - 1
    p = [max(0,-2*(alpha*gama-beta+1)) ; ones(J)*gama]
    q = [ones(J)*gama ; +Inf]
    phi_s_star = [phi_s_star; +Inf]
    admissible_regions = findall((phi_s_star[1:end-1] .< (log_epsilon - log(eps()))/t) .& (phi_s_star[1:end-1] .< phi_s_star[2:end]))
    JJ1 = admissible_regions[end]
    mu_vett = ones(ComplexF64, JJ1)*Inf
    N_vett = ones(ComplexF64, JJ1)*Inf
    h_vett = ones(ComplexF64, JJ1)*Inf
    find_region = 0

    while find_region==0
        for j1 = admissible_regions
            muj=0; hj=0; Nj=0
            if j1 < J1
                j1=Int64(j1)
                (muj, hj, Nj) = OptimalParam_RB(t, phi_s_star[j1], phi_s_star[j1+1], p[j1], q[j1], log_epsilon) ;
            else
                (muj, hj, Nj) = OptimalParam_RU(t,phi_s_star[j1],p[j1],log_epsilon) ;
            end
            mu_vett[j1] = muj; h_vett[j1] = hj; N_vett[j1] = Nj
        end
        if minimum(real.(N_vett[:])) > 200
            log_epsilon = log_epsilon + log(10)
        else
            find_region = 1
        end
    end

    (N, iN) = findmin(real.(N_vett)); mu = mu_vett[iN]; h = h_vett[iN];

    k = collect(-N:N)
    u = h.*k
    z = mu*(1*im*u.+1).^2
    zd = -2*mu*u .+ 2*mu*1*im
    zexp = exp.(z*t)
    F = z.^(alpha*gama-beta)./(z.^alpha .- lambda).^gama.*zd
    S = zexp.*F
    Integral = h*sum(S)/2/pi/im

    ss_star = s_star[(iN[1]+1):end]
    Residues = sum(1/alpha*(ss_star).^(1-beta).*exp.(t*ss_star))

    E = Integral + Residues
    if isreal(lambda)
        E = real(E)
    end
    return E
end

function OptimalParam_RU(t, phi_s_star_j, pj, log_epsilon)
    sq_phi_s_star_j = sqrt(phi_s_star_j)
    phi_s_star_j > 0 ? phibar_star_j = phi_s_star_j*1.01 : phibar_star_j = 0.01
    sq_phibar_star_j = sqrt(phibar_star_j)

    f_min=1
    f_max=10
    f_tar=5

    stop=0
    sq_muj=0
    A=0
    Nj=0

    while stop==0
        phi_t = phibar_star_j*t
        log_eps_phi_t = log_epsilon/phi_t
        Nj = Complex(ceil(real(phi_t/pi*(1-3*log_eps_phi_t/2+sqrt(Complex(1-2*log_eps_phi_t))))), ceil(imag(phi_t/pi*(1-3*log_eps_phi_t/2+sqrt(Complex(1-2*log_eps_phi_t))))))
        A = pi*Nj/phi_t
        sq_muj = sq_phibar_star_j*abs(4-A)/abs(7-sqrt(1+12*A))
        fbar=((sq_phibar_star_j-sq_phi_s_star_j)/sq_muj)^(-pj)
        stop=(pj<1e-14) || (f_min<fbar && fbar<f_max)
        if stop == 0
            sq_phibar_star_j = f_tar^(-1/pj)*sq_muj+sq_phi_s_star_j
            phibar_star_j = sq_phibar_star_j^2
        end
    end
    muj = sq_muj^2
    hj=(-3*A-2+2*sqrt(1+12*A))/(4-A)/Nj

    log_eps=log(eps())
    threshold = (log_epsilon-log_eps)/t
    if muj > threshold
        abs(pj)<1e-14 ? Q=0 : Q=f_tar^(-1/pj)*sqrt(muj)
        phibar_star_j = (Q+sqrt(phi_s_star_j))^2
        if phibar_star_j<threshold
            w=sqrt(log_eps/(log_eps-log_epsilon))
            u=sqrt(-phibar_star_j*t/log_eps)
            muj = threshold
            Nj=ceil(w*log_epsilon/2/pi/(u*w-1))
            hj=sqrt(log_eps/(log_eps-log_epsilon))/Nj
        else
            Nj=+Inf
            hj=0
        end
    end
    return muj, hj, Nj
end

function OptimalParam_RB(t, phi_s_star_j, phi_s_star_j1, pj, qj, log_epsilon)
    log_eps = -36.043653389117154
    fac = 1.01
    conservative_error_analysis = 0

    f_max = exp(log_epsilon - log_eps)

    sq_phi_star_j = sqrt(phi_s_star_j)
    threshold = 2*sqrt((log_epsilon - log_eps)/t)
    sq_phi_star_j1 = min(sqrt(phi_s_star_j1), threshold-sq_phi_star_j)

    if pj < 1.0e-14 && qj < 1.0e-14
        sq_phibar_star_j = sq_phi_star_j
        sq_phibar_star_j1 = sq_phi_star_j1
        adm_region = 1
    end

    if pj < 1.0e-14 && qj >= 1.0e-14
        sq_phibar_star_j = sq_phi_star_j
        sq_phi_star_j > 0 ? f_min = fac*(sq_phi_star_j/(sq_phi_star_j1-sq_phi_star_j))^qj : f_min = fac
        if f_min < f_max
            f_bar = f_min + f_min/f_max*(f_max-f_min)
            fq = f_bar^(-1/qj)
            sq_phibar_star_j1 = (2*sq_phi_star_j1-fq*sq_phi_star_j)/(2+fq)
            adm_region = 1
        else
            adm_region = 0
        end
    end

    if pj >= 1.0e-14 && qj < 1.0e-14
        sq_phibar_star_j1 = sq_phi_star_j1
        f_min = fac*(sq_phi_star_j1/(sq_phi_star_j1-sq_phi_star_j))^pj
        if f_min < f_max
            f_bar = f_min + f_min/f_max*(f_max-f_min)
            fp = f_bar^(-1/pj)
            sq_phibar_star_j = (2*sq_phi_star_j+fp*sq_phi_star_j1)/(2-fp)
            adm_region = 1
        else
            adm_region = 0
        end
    end

    if pj >= 1.0e-14 && qj >= 1.0e-14
        f_min = fac*(sq_phi_star_j+sq_phi_star_j1)/(sq_phi_star_j1-sq_phi_star_j)^max(pj,qj) ;
        if f_min < f_max
            f_min = max(f_min,1.5)
            f_bar = f_min + f_min/f_max*(f_max-f_min)
            fp = f_bar^(-1/pj)
            fq = f_bar^(-1/qj)
            conservative_error_analysis==0 ? w = -phi_s_star_j1*t/log_epsilon : w = -2*phi_s_star_j1*t/(log_epsilon-phi_s_star_j1*t)
            den = 2+w - (1+w)*fp + fq
            sq_phibar_star_j = ((2+w+fq)*sq_phi_star_j + fp*sq_phi_star_j1)/den
            sq_phibar_star_j1 = (-(1+w)*fq*sq_phi_star_j + (2+w-(1+w)*fp)*sq_phi_star_j1)/den
            adm_region = 1
        else
            adm_region = 0
        end
    end
    if adm_region==1
        log_epsilon = log_epsilon - log(f_bar)
        conservative_error_analysis==0 ? w = -sq_phibar_star_j1^2*t/log_epsilon : w = -2*sq_phibar_star_j1^2*t/(log_epsilon-sq_phibar_star_j1^2*t)
        muj = (((1+w)*sq_phibar_star_j + sq_phibar_star_j1)/(2+w))^2
        hj = -2*pi/log_epsilon*(sq_phibar_star_j1-sq_phibar_star_j)/((1+w)*sq_phibar_star_j + sq_phibar_star_j1)
        Nj = Complex(ceil(real(sqrt(Complex(1-log_epsilon/t/muj))/hj)), ceil(imag(sqrt(Complex(1-log_epsilon/t/muj))/hj)))
    else
        muj = 0; hj = 0; Nj = +Inf
    end
    return muj, hj, Nj
end

function mlds(z,al,be,k)
    max_gamma_arg = 171.624
    Jmax = floor(Int, (max_gamma_arg - be)/al)
    G = gamma.(al.*(collect(0:Jmax)).+be)
    jj = collect(k:Jmax)
    f = ones(size(jj))
    for l = 1 : k
        f = f.*(jj.-l.+1)
    end
    c = f./G[k+1:Jmax+1]
    E = zeros(size(z))
    Err_Round = zeros(size(z))
    Err_Round1 = copy(Err_Round)
    Err_Round2 = copy(Err_Round)
    for n = 1 : length(z)
        if abs(z[n]) < eps()
            E[n] = factorial(k)/G[k+1]
        else
            sum_arg = c.*z[n].^(jj.-k)
            abs_sum_arg = abs.(sum_arg)
            i_abs_sum_arg = abs_sum_arg .> eps()/2
            if any(i_abs_sum_arg) == 0
                i_abs_sum_arg[1] = 1
            end
            abs_sum_arg = abs_sum_arg[i_abs_sum_arg]
            sum_arg = sum_arg[i_abs_sum_arg]
            i_abs_sum_arg = sortperm(abs_sum_arg)
            abs_sum_arg = sort(abs_sum_arg)
            sum_arg = sum_arg[i_abs_sum_arg]
            if length(sum_arg) == 1
                E[n] = sum_arg
            else
                S = cumsum(sum_arg)
                S = S[2:end]
                E[n] = S[end]
            end
            J = length(sum_arg) - 1
            JJ = [J; collect(J:-1:1)]
            Err_Round1[n] = sum(JJ.*abs_sum_arg)*eps()
            length(sum_arg)==1 ? Err_Round2[n] = Err_Round1[n] : Err_Round2[n] = sum(abs.(S))*eps()
            Err_Round[n] = exp((log(Err_Round1[n])+log(Err_Round2[n]))/2)
        end
    end
    return E, Err_Round
end

function mldr(t, s, α, β, k)
    omega = zeros(k+2)
    omega[1] = 1
    jj = collect(0:k)
    p = (α.-jj)
    pr = 1
    for j = 1:k+1
        pr = pr*p[j]
        omega[j+1] = pr/factorial(j)
    end
    Hk = zeros(k+1)
    Hk[1] = 1
    for j = 1:k
        ll = collect(1:j)
        Hk[j+1] = -1 ./α.*sum(omega[ll.+2].*(k.*ll./j.+1).*Hk[j.-ll.+1])
    end
    ck = zeros(k+1)
    for j = 0:k
        temp = 0
        for l = 0 : k-j
            if l == 0
                p = 1
            else
                ll = 0:l-1
                p = prod(α.-β.-ll)
            end
            temp = temp .+ p*Hk[k-j-l+1]./factorial(l)
        end
        ck[j+1] = temp/factorial(j)
    end
    tep = Polynomial(ck)
    result = tep.(s)
    R = 1 ./α^(k+1)*exp.(t.*s).*(s//1)^(1-α.*k.-β).*result
    return R
end
#=
"""
```tex
@article{Garrappa2018ComputingTM,
  title={Computing the Matrix Mittag-Leffler Function with Applications to Fractional Calculus},
  author={Roberto Garrappa and Marina Popolizio},
  journal={Journal of Scientific Computing},
  year={2018},
  volume={77},
  pages={129-153}
}
```
"""
function mldlt(z,alpha,beta,k)
    log_epsilon = log(10^(-15)) ;
    E = zeros(size(z)) ;
    for j = 1 : length(z)
        if abs(z[j] < 1.0e-14)
            E[j] = factorial(k)/gamma(alpha*k+beta)
        else
            E[j] = LTIversion(1,z[j],alpha,beta,k,log_epsilon) ;
        end
        if isreal(z[j])
            E[j] = real(E[j])
        end
    end
    return E
end

function mld(z, alpha, beta, k)
    Ek = zeros(size(z))
    tau = 1.0e-14 ;
    max_gamma_arg = 171.624 ;
    Jmax = floor((max_gamma_arg - beta)/alpha) ;
    z_abs_max = (tau*gamma(alpha*Jmax+beta)/prod(Jmax.-collect(0:k-1))).^(1/(Jmax-k)) ;
    i_z_se = abs.(z).<z_abs_max ;
    z_se = z[i_z_se] ;
    E_se = zeros(size(z)) ; Err_Round = zeros(size(z)) ; e = ones(size(z)) ;
    (E_se[i_z_se], Err_Round[i_z_se]) = mlds(z_se,alpha,beta,k) ;
    i_z_se_accept = (Err_Round./(e.+abs.(E_se)) .< tau) .& ~(E_se==0)  ;
    Ek[i_z_se_accept] = E_se[i_z_se_accept]
    i_ze_lt = .!i_z_se_accept ;
    if k <= 3
        p = 0
    elseif k <= 7
        p = 1
    else
        p = 2
    end
    c = zeros(k-p+1, k-p+1) ;
    c[1, 1] = 1 ;
    for kk = 1 : k-p
        c[kk+1,1] = (1-alpha*(kk-1)-beta)*c[kk,1]
        for j = 1 : kk-1
            c[kk+1,j+1] = c[kk,j] + (j+1-alpha*(kk-1)-beta)*c[kk,j+1]
        end
        c[kk+1,kk+1] = 1 ;
    end
    for j = 0 : k-p
        if abs(c[k-p+1,j+1]) > 1.0e-14
            Ek[i_ze_lt] = Ek[i_ze_lt] .+ c[k-p+1,j+1]*mldlt(z[i_ze_lt], alpha, (k-p)*alpha+beta-j, p) ;
        end
    end
    Ek[i_ze_lt] = Ek[i_ze_lt]./alpha^(k-p) ;
    return Ek
end

function mittleff(A, alpha, beta)
    fun(z, k) = mld(z, alpha, beta, k)
    E=funm(A)
end
#FIXME: Didn't see any matrix function support in Julia😟.
# Yingbo has an implementation: https://github.com/YingboMa/Funm.jl
=#


text:
module UnicodePlotsSimple

using Crayons

export histogram
export distribution
export heatmap
export densitymap

eps(::Type{Any}) = Base.eps(Float32)
eps(::Type{Int64}) = Base.eps(Float32)
eps(::Type{Float32}) = Base.eps(Float32)
eps(::Type{Float64}) = Base.eps(Float64)
cpad(s, n::Integer, p=" ") = rpad(lpad(s,div(n+length(s),2),p),n,p)

# HEATMAPS styles
marker(format::Val{:HEATMAP}, cases, x,y,c) = "  "
marker(format::Val{:ONLYNUMBER}, cases, x,y,c) = lpad(string(max(cases[c][x, y],0) > 99 ? "…" : cases[c][x, y]),2)
marker(format::Val{:NUMBERANDCOLOR}, cases, x,y,c) = lpad(string(max(cases[c][x, y],0) > 99 ? "…" : cases[c][x, y]),2)

# DENSITY styles
marker(format::Val{:NUMBER}, cases, x,y,c) = lpad(string(max(cases[c][x, y],0) > 99 ? "…" : cases[c][x, y]),2)
marker(format::Val{:DISTRIBUTION}, cases, x,y,c) = lpad(string(max(cases[c][x, y],0) > 99 ? "…" : cases[c][x, y]),2)
marker(format::Val{:DENSITY}, cases, x,y,c) = ["░░", "▒▒", "▓▓", "██"][cases[c][x, y]]

# HEATMAPS styles
markercolor(format::Val{:HEATMAP}, c) = Crayon(;background=c)
markercolor(format::Val{:ONLYNUMBER}, c) = Crayon(;bold=true)
markercolor(format::Val{:NUMBERANDCOLOR}, c) = Crayon(;background=c,foreground=(33,33,33), bold=true)

# DENSITY styles
markercolor(format::Val{:NUMBER}, c) = Crayon(;background=c,foreground=(33,33,33), bold=true)
markercolor(format::Val{:DISTRIBUTION}, c) = Crayon(;background=c,foreground=(33,33,33), bold=true)
markercolor(format::Val{:DENSITY}, c) = Crayon(;foreground=c)

heatmap(xy::Matrix; format=Val(:HEATMAP) ,title=[], xticks=[], yticks=[], reversey=false) = heatmap([xy], title=title, xticks=xticks, yticks=yticks, reversey=reversey, format=format)
heatmap(xy::Vector; format=Val(:HEATMAP) ,title=[], xticks=[], yticks=[], reversey=false) = heatmap(xy,title,xticks,yticks,reversey, format)
function heatmap(xy::Vector, title, xticks, yticks, reversey, format)
	densitymap(xy, nothing, title, xticks, yticks, reversey, format)
end

densitymap(xy::Matrix, case::Matrix; title=[], xticks=[], yticks=[], reversey=false, format=Val(:DENSITY)) = densitymap([xy], [case], title=title, xticks=xticks, yticks=yticks, reversey=reversey, format=format)
densitymap(xy::Vector, cases; title=[], xticks=[], yticks=[], reversey=false, format=Val(:DENSITY)) =
densitymap(xy,cases,title,xticks,yticks,reversey,format)
function densitymap(xy::Vector{Matrix{T}}, cases, title, xticks, yticks, reversey, format) where T
	ϵ = eps(eltype(xy[1]))
	scale = 0:5:255
	charts_num = length(xy)
	width, height = size(xy[1])
	# freq_matrix = Array{Int,3}(undef, width, height, charts_num)
	unicode_matrix = Array{UInt8,3}(undef, width, height, charts_num)
	extremas = Tuple{T,T}[]
	for c in 1:charts_num
		minstrength, maxstrength = extrema(xy[c])
		push!(extremas, (minstrength, maxstrength))
		if minstrength == maxstrength
			unicode_matrix[:,:,c] .= UInt8(0)
		else
			for xi in 1:size(xy[c],1)
				for yi in 1:size(xy[c],2)
					strength = (xy[c][xi,yi]-minstrength) / ((maxstrength-minstrength) * (1+ϵ)) * length(scale)
					unicode_matrix[xi,yi,c] = convert(UInt, scale[floor(Int, strength)+1])
				end
			end
		end
	end

	length(title)!= 0 && (print(" ");for c in 1:charts_num
		print(cpad(title[c],2*length(xticks)))
	end;println())

	if isa(format, Val{:ONLYNUMBER}) || isa(format, Val{:NUMBERANDCOLOR})
		cases=[Int.(floor.(m.*1000)) for m in xy]  # we multiplicate with 1000 to avoid disappearing 0.01 and 0.2 and so on... (I know it is not perfect but for now it is enough for me...)
	end
	if isa(format, Val{:DISTRIBUTION}) || isa(format, Val{:ONLYNUMBER}) || isa(format, Val{:NUMBERANDCOLOR})
		maxcase = maximum.(cases)
		scaled_cases = [maxcase[c]>0 ? floor.(Int, case_m ./ (maxcase[c]*(1+ϵ)) * 98) .+ ifelse.(case_m .> 0,1,0) : zero(case_m) for (c,case_m) in enumerate(cases)]
	elseif isa(format, Val{:DENSITY})
		maxcase = maximum.(cases)
		scaled_cases = [maxcase[c]>0 ? floor.(Int, case_m ./ (maxcase[c]*(1+ϵ)) * 4) .+ 1 : zero(case_m) for (c,case_m) in enumerate(cases)]
	else
		scaled_cases = cases
	end

	for h in 1:height
		length(yticks) != 0 && print(yticks[reversey ? height+1-h : h])
		for c in 1:charts_num
			for w in 1:width
				# @show h,c,w, height
				m = marker(format, scaled_cases, w,reversey ? height+1-h : h,c, )
				rgb=(255, unicode_matrix[w,reversey ? height+1-h : h,c], 0)
				cray_col = markercolor(format, rgb)
				print(cray_col, m, Crayon(reset=true))
			end
			print(" ")
		end

		if h==1
			print(" " ^ (1), "min   …    max ");
		else
			length(extremas) >= h-1 && (print(" ");
			printstyled(lpad(round(extremas[h-1][1]; digits=6), 8); color=:green, bold=true);
			print(" … ");
			printstyled(lpad(round(extremas[h-1][2]; digits=8), 6); color=:green, bold=true))
		end
		println()
	end
	length(xticks) != 0 && (print("");
	for c in 1:charts_num
		for w in 1:width print(lpad(xticks[w],2)) end ; print(" ")
	end)
	println()
end

distribution(probs; height=3, values=[]) = distribution(probs, height, values)
function distribution(probs, height, values)
	width = displaysize(stdout)[2] # terminal_width
	ϵ = eps(eltype(probs))
	histbars_blocks = [' ', '▁', '▂', '▃', '▄', '▅', '▆', '▇', '█']
	blocknum = length(histbars_blocks)
	histbin_codes = fill('█',(blocknum-1)*height, height)
	for h in 1:height
		histbin_codes[(h-1)*(blocknum-1)+1:h*(blocknum-1),height-h+1] .= histbars_blocks[2:end]
	end
	histbin_codes[map(ci -> cld(ci[1],blocknum-1)+ci[2]<=height,CartesianIndices(histbin_codes))].= ' '

	maxheight = maximum(probs)
	# maxheight = isnan(maxheight) ? 1 : maxheight

	unicode_matrix = Matrix{Char}(undef, length(probs), height)
	for (i,bin) in enumerate(probs)
		strength = bin / (maxheight *(1+ ϵ)) * size(histbin_codes,1)
		unicode_matrix[i, :] .= histbin_codes[floor(Int, strength)+1,:]
	end
	v_mean = sum(probs) / length(probs)
	v_σ = length(probs) > 1 ? sqrt(sum((probs .- v_mean).^2) / (length(probs)-1)) : 0f0
	v_min = minimum(probs)
	v_max = maximum(probs)
	for wx in 1:width:length(probs)
		for h in 1:height
			print(join(unicode_matrix[wx:min(wx+width-1,end),h],""))
			# if wx>length(probs)-width && size(unicode_matrix,1)- wx + 2+9+2+8+3+8 > terminal_width
				if wx>length(probs)-width && h==height-2
					print("  mean ± σ:  ")
					printstyled(lpad(round(v_mean; digits=6), 8); color=:green, bold=true)
					print(" ± ")
					printstyled(rpad(round(v_σ; digits=6), 8); color=:green, bold=true)
				elseif wx>length(probs)-width && h==height-1
					print("  min … max: ")
					printstyled(lpad(round(v_min; digits=6), 8); color=:green, bold=true)
					print(" … ")
					printstyled(rpad(round(v_max; digits=6), 8); color=:green, bold=true)
				elseif wx>length(probs)-width && h==height && length(values)>0
					if length(values)<20
						print("  values:    ")
						printstyled(join(values,", "); color=:green, bold=true)
					else
						print("  values: Too many… ")
					end
				end
			# end
			println()
		end
	end
end

histogram(values; width=64, height=3, title="", printstat=true, reversescale=false, line1="", line2="", line3="", fixleftscale=-123444f0, fixrightscale=-123444f0) = histogram(values, width, height, title, printstat, reversescale, line1, line2, line3, fixleftscale, fixrightscale)
function histogram(values, width, height, title, printstat, reversescale, line1, line2, line3, fixleftscale, fixrightscale)
	ϵ = eps(eltype(values))
	length(values) == 0 && (println("$title  $line1  $line2 ---Empty histogram--- "); return)
	histbars_blocks = [' ', '▁', '▂', '▃', '▄', '▅', '▆', '▇', '█']
	blocknum = length(histbars_blocks)
	histbin_codes = fill('█',(blocknum-1)*height, height)
	for h in 1:height
		histbin_codes[(h-1)*(blocknum-1)+1:h*(blocknum-1),height-h+1] .= histbars_blocks[2:end]
	end
	histbin_codes[map(ci -> cld(ci[1],blocknum-1)+ci[2]<=height,CartesianIndices(histbin_codes))].= ' '
	# for h in 1:height
	# 	println(join(histbin_codes[:,h],"")," █")
	# end
	firstbin = fixrightscale == -123444f0 ? minimum(values) : fixrightscale
	endbin = fixleftscale == -123444f0 ? maximum(values) : fixleftscale

	binsize = firstbin == endbin ? 1. : (endbin - firstbin) / width * (1 + ϵ)
	histogram_bins = zeros(Int, width)
	for v in values
		!(0 < floor(Int,(v-firstbin) / binsize)+1 <= width) && continue
		histogram_bins[floor(Int,(v-firstbin) / binsize)+1] += 1
	end
	maxbin = maximum(histogram_bins)
	unicode_matrix = Matrix{Char}(undef, width, height)
	for (i,bin) in enumerate(histogram_bins)
		strength = floor(Int, bin / (maxbin * (1 + ϵ))* size(histbin_codes,1))
		unicode_matrix[i, :] .= histbin_codes[strength+1,:]
	end

	title !=="" && printstyled(" " ^ max(0,cld((width-length(title)),2)),title,"\n"; bold=true)
	for h in 1:height
		if reversescale == false
			print(join(unicode_matrix[:,h],""))
		else
			print(join(unicode_matrix[end:-1:1,h],""))
		end
		if h==1
      print(" ",line1)
		elseif h==2
      print(" ",line2)
		elseif h==3
      print(" ",line3)
		end
		println()
	end
	if reversescale
		firstbin, endbin = endbin, firstbin
	end
	printstyled(rpad(round(firstbin; digits=6),12), bold=true)
	# sign_width = (sign(firstbin) < 0 ? 1 : 0) + (sign(endbin) < 0 ? 1 : 0)
	# print(" "^(width-14-Int(ceil(log10(abs(firstbin)))+ceil(log10(abs(endbin))))- sign_width))
	print(" "^(width-24))
	printstyled(lpad(round(endbin; digits=6),12), bold=true)
	# println()
	v_mean = sum(values) / length(values)
	v_σ = sqrt(sum((values .- v_mean).^2) / length(values))
	if printstat
		print(" ")
		print("Mean ± σ: ")
		printstyled(lpad(round(v_mean; digits=5), 6); color=:green, bold=true)
		print(" ± ")
		printstyled(round(v_σ; digits=4); color=:green, bold=true)
		if fixrightscale == -123444f0 && fixleftscale == -123444f0
			println()
		else
			println(" Cutted: $(length(values)-sum(histogram_bins))")
		end
		# print("Range (min … max): ")
		# printstyled(lpad(round(firstbin; digits=6), 6); color=:green, bold=true)
		# print(" … ")
		# printstyled(lpad(round(endbin; digits=6), 6); color=:green, bold=true)
		println()
	end
	# v_mean, v_σ, firstbin, endbin
	return
end

function histogramdetailed(values, width=64, height=3)
	# TODO implement fined detailed histogram print...
	@assert false "Unimplemented fine grained sparseprint..."
	# 30×30 SparseMatrixCSC{Bool, Int64} with 445 stored entries:
	# ⡙⠅⡥⠹⢗⢊⢣⣽⠜⢬⣛⠇⡦⡮⣑
	# ⢦⣷⡏⣸⢐⠺⣲⢣⣵⣀⡞⢹⢏⣃⣁
	# ⡟⣻⢹⣅⡻⣌⠄⠳⣱⣴⠖⠰⣢⠣⠵
	# ⡸⡑⠜⢰⢠⣢⢤⢎⡷⠽⡩⣒⡢⡊⢋
	# ⡀⠀⠇⢟⢓⣊⣊⢱⠙⡯⡅⢇⠀⣁⠍
	# ⡢⡡⡦⠦⠇⣢⡉⣻⣣⣽⢥⠽⢯⣌⢵
	# ⠞⡸⣡⠲⠪⢟⠹⡀⢍⡢⡍⢼⢛⡅⢃
	# ⠉⠀⠙⠊⠛⠓⠚⠂⠈⠂⠁⠓⠙⠈⠁
end

end # module


text:
@muladd function calc_twopoint_random!(_dW, sqdt, dW)
  @.. _dW = ifelse(sign(dW) > false, sqdt, -sqdt)
  return nothing
end

@muladd function calc_twopoint_random(sqdt, dW)
  return ifelse(sign(dW) > false, sqdt, -sqdt)
end

function calc_threepoint_random!(_dW, sq3dt, quantile, dW_scaled)
  @. _dW = ifelse(abs(dW_scaled) > -quantile, ifelse(dW_scaled < quantile, -sq3dt, sq3dt), zero(sq3dt))
  return nothing
end

@muladd function calc_threepoint_random(sq3dt, quantile, dW_scaled)
  return ifelse(abs(dW_scaled) > -quantile, ifelse(dW_scaled < quantile, -sq3dt, sq3dt), zero(sq3dt))
end


text:
# %% imports

using Cambrian
using EvolutionaryStrategies

include("../Actors/cartpole_tryout.jl")

# %% build a class with all the matrices, the fitness associated ( mean reward from cartpole tryout )

mutable struct MatrixClust

    fitness::Array{Float64}
    matrices::Array{Array{Float64,2}}
    s_to_m::Array{Int64}
    segs::Array{Array{Float64}}

    function MatrixClust(elites::Array{AbstractESIndividual,1}, mean_score::Array{Float64})
        fitness_tmp = [elites[i].fitness[1] for i in 1:length(elites)]
        ar = [map(x->borne(x,-1.0, 1.0),elites[i].genes[1:8*4]) for i in 1:length(elites)]
        matrices_tmp = [reshape(ar[i], (8,4)) for i in 1:length(ar)]

        idx = findall(x -> x > 0., mean_score .- 195.)
        fitness = fitness_tmp[idx]
        matrices = matrices_tmp[idx]
        segs = Array{Float64}[]
        s_to_m = Int64[]
        for i in eachindex(matrices)
            m = matrices[i]
            for j in 1:size(m)[1]
                push!(segs, m[j,:])
                push!(s_to_m, i)
            end
        end
        new(fitness, matrices, s_to_m, segs)
    end

    function MatrixClust()
        fitness = Float64[]
        matrices = Array{Float64,2}[]
        s_to_m = Int64[]
        segs = Array{Float64}[]
        new(fitness, matrices, s_to_m, segs)
    end
end

function join_matclust(mc1::MatrixClust, mc2::MatrixClust)
    l = length(mc1.matrices)
    mc1.fitness = [mc1.fitness ; mc2.fitness]
    mc1.matrices = [mc1.matrices ; mc2.matrices]
    mc1.s_to_m = [mc1.s_to_m ; mc2.s_to_m .+ l]
    mc1.segs = [mc1.segs ; mc2.segs]
    return mc1
end

# %% main function

# Builds the class MatrixClust with all the best individual from cartpole sNES
# saves it as mc.jl
function build_matclust(n::Int64)
    mc = MatrixClust()
    for i in 1:n
        elites = load("Julia/Actors/cartpoleElites/sNESelites_cartpole_$(i).jld", "elites")
        means = test_random(elites)
        mc_tmp = MatrixClust(elites, means)
        mc = join_matclust(mc, mc_tmp)
    end
    save("Julia/MatrixCluster/mc.jld", "mc", mc)
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
md"""
# Crossover Operations

Crossover is the way two candidates are combined to create new candidates. In NaiveGAflux crossover always maps two candidates into two new candidates. Just as for mutation, NaiveGAflux does this while preserving (to whatever extent possible) the parameters and alignment between layers of the combined models.

Crossover operations might not seem to make much sense when using parameter inheritance (i.e the concept that children retain the parameters of their parents). Randomly combining layers from two very different models will most likely not result in a well performing model. There are however a few potentially redeeming effects:

* Early in the evolution process parameters are not yet well fitted and inheriting parameters is not worse than random initialization
* A mature population on the other hand will consist mostly of models which are close relatives and therefore have somewhat similar weights.

Whether these effects actually make crossover a genuinely useful operation when evolving neural networks is not yet proven though. For now it is perhaps best to view the crossover operations as being provided mostly for the sake of completeness.

The following basic crossover operations are currently supported:
1. Swap segments between two models using [`CrossoverSwap`](@ref).
2. Swap optimizers between two candidates using [`OptimizerCrossover`](@ref).
3. Swap learning rate between two candidates using [`LearningRateCrossover`](@ref).

Most of the mutation utilities also work with crossover operations. Here are a few examples:
"""

@testset "Crossover examples" begin #src
# Start by creating a model to play with.
Random.seed!(NaiveGAflux.rng_default, 0) #src

invertex = denseinputvertex("A.in", 3)
layer1 = fluxvertex("A.l1", Dense(nout(invertex), 4), invertex; layerfun=ActivationContribution)
layer2 = fluxvertex("A.l2", Dense(nout(layer1), 5), layer1; layerfun=ActivationContribution)
layer3 = fluxvertex("A.l3", Dense(nout(layer2), 3), layer2; layerfun=ActivationContribution)
layer4 = fluxvertex("A.l4", Dense(nout(layer3), 2), layer3; layerfun=ActivationContribution)
modelA = CompGraph(invertex, layer4)

# Create an exact copy to show how parameter alignment is preserved.
# Prefix names with B so we can show that something actually happened.
import Functors
modelB = Functors.fmap(x -> x isa String ? replace(x, r"^A.\.*" => "B.") : x, modelA)

indata = reshape(collect(Float32, 1:3*2), 3,2)
@test modelA(indata) == modelB(indata)

@test name.(vertices(modelA)) == ["A.in", "A.l1", "A.l2", "A.l3", "A.l4"]
@test name.(vertices(modelB)) == ["B.in", "B.l1", "B.l2", "B.l3", "B.l4"]

# [`CrossoverSwap`](@ref) takes ones vertex from each graph as input and swaps a random segment from each graph.
# By default it tries to make segments as similar as possible
swapsame = CrossoverSwap()

swapA = vertices(modelA)[4]
swapB = vertices(modelB)[4]
newA, newB = swapsame((swapA, swapB))

# It returns vertices of a new graph to be compatible with mutation utilities.
# Parent models are not modified.
@test newA ∉ vertices(modelA)
@test newB ∉ vertices(modelB)

# The function `regraph` is an internal utility which should not be needed in normal use cases, but here we use it to make comparison easier.
modelAnew = NaiveGAflux.regraph(newA)
modelBnew = NaiveGAflux.regraph(newB)

@test name.(vertices(modelAnew)) == ["A.in", "A.l1", "B.l2", "B.l3", "A.l4"]
@test name.(vertices(modelBnew)) == ["B.in", "B.l1", "A.l2", "A.l3", "B.l4"]

@test modelA(indata) == modelB(indata) == modelAnew(indata) == modelBnew(indata)

# Deviation parameter will randomly make segments unequal.
swapdeviation = CrossoverSwap(0.5)
modelAnew2, modelBnew2 = NaiveGAflux.regraph.(swapdeviation((swapA, swapB)))

@test name.(vertices(modelAnew2)) == ["A.in", "A.l1", "A.l2", "B.l1", "B.l2", "B.l3", "A.l4"]
@test name.(vertices(modelBnew2)) == ["B.in", "A.l3", "B.l4"]

# [`VertexCrossover`](@ref) applies the wrapped crossover operation to all vertices in a `CompGraph` and is the main API for doing crossover.
# It in addtion, it selects compatible pairs for us (i.e `swapA` and `swapB` above).
# It also takes an optional deviation parameter which is used when pairing.
crossoverall = VertexCrossover(swapdeviation, 0.5)

modelAnew3, modelBnew3 = crossoverall((modelA, modelB))

# I guess things got swapped back and forth so many times not much changed in the end.
@test name.(vertices(modelAnew3)) == ["A.in", "A.l2", "A.l4"]
@test name.(vertices(modelBnew3)) ==  ["B.in", "B.l3", "B.l1", "B.l2", "A.l1", "A.l3", "B.l4"]

# As advertised above, crossovers interop with most mutation utilities, just remember that input is a tuple for things which require a callback.
# Perform the swapping operation with a 30% probability for each valid vertex pair and log which vertices are being swapped.
crossoversome = VertexCrossover(
                    MutationProbability(
                        LogMutation(
                            ((v1,v2),) -> "Swap $(name(v1)) and $(name(v2))",
                            swapdeviation),
                    0.3)
                )

@test_logs( (:info, "Swap A.l1 and B.l1"),
            (:info, "Swap A.l2 and B.l2"),
            crossoversome((modelA, modelB)))
end #src


text:
module CBFToolbox

# Julia packages required for this module
using LinearAlgebra
using Zygote
using JuMP
using OSQP
using Plots
using PlotThemes
import DifferentialEquations: ODEProblem, solve

# Abstract types
abstract type System end
abstract type Policy end
abstract type CertificateFunction end

# Export Simulation type
export Simulation
export Trajectory

# Export systems
export System
export ControlAffineSystem
export step!
export initialize!

# Export Certificate functions
export LyapunovFunction
export ControlLyapunovFunction
export TimeVaryingCLF
export BarrierFunction
export ControlBarrierFunction
export HighOrderCBF

# Export control policies
export Policy
export FeedbackPolicy
export CLFQP
export TimeVaryingCLFQP
export CBFQP
export TimeVaryingCBFQP

# Export utility functions and types
export CircularObstacle
export custom_plots
export custom_theme
export circle_shape
export vec2mat

# System definitions
include("Systems/control_affine_system.jl")
include("Systems/feedback_policy.jl")
include("Systems/system_library.jl")

# Lyapunov functions
include("LyapunovFunctions/control_lyapunov_functions.jl")
include("LyapunovFunctions/time_varying_clfs.jl")
include("LyapunovFunctions/clf_quadratic_programs.jl")
include("LyapunovFunctions/time_varying_clfqp.jl")

# Barrier functions
include("BarrierFunctions/control_barrier_functions.jl")
include("BarrierFunctions/cbf_quadratic_programs.jl")
include("BarrierFunctions/high_order_cbfs.jl")
include("BarrierFunctions/hocbf_quadratic_programs.jl")
include("BarrierFunctions/time_varying_cbfqp.jl")

# Base simulation type
include("Simulations/simulation.jl")

# Utilities
include("Utils/data_utils.jl")
include("Utils/cbf_utils.jl")
include("Utils/plot_utils.jl")

end


text:
using Test

using Logging
using Scruff
using Scruff.Models
using Scruff.Utils
using Scruff.RTUtils
using Scruff.SFuncs
using Scruff.Operators

@testset "test loggers" begin
    x1 = Cat([1,2], [0.1, 0.9])()(:x1)
    x2 = Cat([1,2,3], [0.2, 0.3, 0.5])()(:x2)
    cpd2 = Dict((1,1) => [0.3, 0.7], (1,2) => [0.6, 0.4], (2,1) =>[0.4, 0.6],
                (2,2) => [0.7, 0.3], (3,1) => [0.5, 0.5], (3,2) => [0.8, 0.2])
    x3 = DiscreteCPT([1,2], cpd2)()(:x3)
    x4 = DiscreteCPT([1,2], Dict((1,) => [0.15, 0.85], (2,) => [0.25, 0.75]))()(:x4)
    x5 = DiscreteCPT([1,2], Dict((1,) => [0.35, 0.65], (2,) => [0.45, 0.55]))()(:x5)
    x6 = DiscreteCPT([1,2], Dict((1,) => [0.65, 0.35], (2,) => [0.75, 0.25]))()(:x6)

    fivecpdnet = InstantNetwork(Variable[x1,x2,x3,x4,x5], VariableGraph(x3=>[x2,x1], x4=>[x3], x5=>[x3]))

    @testset "runtime logger" begin
        run = Runtime(fivecpdnet)
        default_initializer(run)
        inst1 = current_instance(run, x1)
        inst3 = current_instance(run, x3)
        inst5 = current_instance(run, x5)

        @test_logs(
            (:info, r"distribute_messages!"),
            (:info, r"get_variable"),
            (:info, r"collect_messages"),
            (:info, r"get_variables"),
            (:info, r"set_value!"),
            (:info, r"current_instance"),
            match_mode=:any,
            trace_runtime(Scruff.Algorithms.three_pass_BP, run))
        @test_logs(
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            (:info, r"collect_messages"),
            trace_runtime(Scruff.Algorithms.three_pass_BP, run; fnamefilter=x->x == "collect_messages"))
    end

    #= This test takes too long
    @testset "log all" begin
        run = Runtime(fivecpdnet)
        default_initializer(run)
        inst1 = current_instance(run, x1)
        inst3 = current_instance(run, x3)
        inst5 = current_instance(run, x5)

        @test_logs(
            (:info, r"distribute_messages!"),
            (:info, r"has_state"),
            (:info, r"get_network"),
            (:info, r"get_state"),
            (:info, r"collect_messages"),
            (:info, r"current_time"),
            (:info, r"get_variables"),
            (:info, r"set_value!"),
            (:info, r"outgoing_pis"),
            (:info, r"get_parents"),
            (:info, r"get_children"),
            (:info, r"compute_bel"),
            (:info, r"current_instance"),
            match_mode=:any,
            trace_algorithm(Scruff.Algorithms.three_pass_BP, run))
    end
    =#

    @testset "timing logger" begin
        run = Runtime(fivecpdnet)
        default_initializer(run)
        inst1 = current_instance(run, x1)
        inst3 = current_instance(run, x3)
        inst5 = current_instance(run, x5)

        @test_logs(
            (:info, "three_pass_BP"),
            (:info, "outgoing_pis"),
            (:info, "get_sfunc"),
            match_mode=:any,
            time_algorithm(Scruff.Algorithms.three_pass_BP, run))
    end
end


text:
"""
# TwoDots

The goal of TwoDots is to earn as many points as possible by connecting dots of the same color. If the connected dots form a rectangle, all of the dots of that color count towards your points. Click on the dots you want to connect and click submit.

# Example

```julia
using TwoDots
start()
```
"""
module TwoDots
    using Gtk, Cairo, Gtk.ShortNames
    export Dot, Game, GUI, start
    export click_dot!, click_submit!, can_select

    include("structs.jl")
    include("functions.jl")
    include("gui.jl")
end


text:
using Microbiome
using Distances
using LinearAlgebra
using DataFrames
using Clustering
using Random
using Test

@testset "Abundances" begin
    # Constructors
    M = rand(100, 10)
    df = hcat(DataFrame(x=collect(1:100)), DataFrame(M))

    a1 = abundancetable(M)
    a2 = abundancetable(df)
    @test typeof(a1) <:AbstractComMatrix
    @test typeof(a2) <:AbstractComMatrix

    abund = abundancetable(
        M, ["sample_$x" for x in 1:10],
        ["feature_$x" for x in 1:100])

    @test a1.occurrences == a2.occurrences == abund.occurrences

    # Normalization functions
    relab_fract = relativeabundance(abund)
    @test typeof(relab_fract) <: AbstractComMatrix

    relab_perc = relativeabundance(abund, kind=:percent)
    @test typeof(relab_perc) <: AbstractComMatrix

    rnorm = rownormalize(abund)
    cnorm = colnormalize(abund)

    @test size(abund) == (nfeatures(abund), nsamples(abund))
    @test size(relab_fract) == (nfeatures(relab_fract), nsamples(relab_fract))
    @test size(relab_perc) == (nfeatures(relab_perc), nsamples(relab_perc))
    @test size(rnorm) == (nfeatures(rnorm), nsamples(rnorm))
    @test size(cnorm) == (nfeatures(cnorm), nsamples(cnorm))

    for j in 1:10
        @test sum(getsample(relab_fract, j)) ≈ 1
        @test sum(getsample(relab_perc, j)) ≈ 100
    end

    for i in 1:nfeatures(rnorm)
        @test maximum(getfeature(rnorm,i)) ≈ 1
    end

    for j in 1:nsamples(cnorm)
        @test maximum(getsample(cnorm, j)) ≈ 1
    end

    # Filtering Functions
    filt = filterabund(relab_fract, 5)
    @test typeof(filt) <: AbstractComMatrix
    @test typeof(filterabund(df, 5)) <: AbstractComMatrix

    @test size(filt) == (6, 10)
    for i in 1:10
        @test sum(getsample(filt, i)) ≈ 1
    end

    @test featurenames(filt)[end] == "other"
end

@testset "Distances" begin
    # Constructors
    Random.seed!(1)
    M = rand(100, 10)
    df = hcat(DataFrame(x=collect(1:100)), DataFrame(M))
    abund = abundancetable(
        M, ["sample_$x" for x in 1:10],
        ["feature_$x" for x in 1:100])

    dm1 = getdm(M, BrayCurtis())
    dm2 = getdm(df, BrayCurtis())
    dm = getdm(abund, BrayCurtis())

    @test dm.dm == dm1.dm == dm2.dm

    rowdm1 = getrowdm(M, BrayCurtis())
    rowdm2 = getrowdm(df, BrayCurtis())
    rowdm = getrowdm(abund, BrayCurtis())

    @test rowdm.dm == rowdm1.dm == rowdm2.dm

    @test size(dm) == (10, 10)
    @test size(rowdm) == (100, 100)
    for i in 1:10; @test dm[i,i] == 0; end
    for i in 1:100; @test rowdm[i,i] == 0; end

    # PCoA
    p = pcoa(dm, correct_neg=true)
    @test sum(p.variance_explained) ≈ 1
    for i in 1:size(p, 2)
        @test eigenvalue(p, i) > 0
        @test typeof(eigenvalue(p, i)) <: Real
    end

    @test sum([variance(p, i) for i in 1:size(p,2)]) ≈ 1
    @test sort(variance(p, 1:size(p,2)), rev=true) == variance(p, 1:size(p,2))

    @test length(principalcoord(p, 1)) == size(dm, 1)
    @test principalcoord(p, 1:size(p,2)) == p.eigenvectors

    # Diversity indicies
    R = 100
    s1 = rand(R) # high diversity
    s2 = [i % 10 == 0 ? s1[i] : 0 for i in 1:R] # low diversity
    s3 = ones(R) # uniform
    s4 = [1., zeros(R-1)...] # no diversity

    @test shannon(s1) > shannon(s2)
    @test shannon(s3) ≈ log(R)
    @test shannon(s4) ≈ 0.

    @test ginisimpson(s1) > ginisimpson(s2)
    @test ginisimpson(s3) ≈ 1. - 1/R
    @test ginisimpson(s4) ≈ 0.
end

@testset "Leaf Ordering" begin
    Random.seed!(42)
    m = rand(100, 10)

    dm = pairwise(BrayCurtis(), m)
    h = hclust(dm, :single);

    ordered = optimalorder(h, dm)

    @test ordered.order == [7, 3, 1, 9, 2, 6, 10, 4, 5, 8]
    @test ordered.merge == h.merge
end


text:
"""
    EGARCH{o, p, q, T<:AbstractFloat} <: UnivariateVolatilitySpec{T}
"""
struct EGARCH{o, p, q, T<:AbstractFloat} <: UnivariateVolatilitySpec{T}
    coefs::Vector{T}
    function EGARCH{o, p, q, T}(coefs::Vector{T}) where {o, p, q, T}
        length(coefs) == nparams(EGARCH{o, p, q})  || throw(NumParamError(nparams(EGARCH{o, p, q}), length(coefs)))
        new{o, p, q, T}(coefs)
    end
end

"""
    EGARCH{o, p, q}(coefs) -> UnivariateVolatilitySpec

Construct an EGARCH specification with the given parameters.

# Example:
```jldoctest
julia> EGARCH{1, 1, 1}([-0.1, .1, .9, .04])
EGARCH{1,1,1} specification.

─────────────────────────────────
                ω   γ₁   β₁    α₁
─────────────────────────────────
Parameters:  -0.1  0.1  0.9  0.04
─────────────────────────────────
```
"""
EGARCH{o, p, q}(coefs::Vector{T}) where {o, p, q, T}  = EGARCH{o, p, q, T}(coefs)

@inline nparams(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = o+p+q+1
@inline nparams(::Type{<:EGARCH{o, p, q}}, subset) where {o, p, q} = isempty(subset) ? 1 : sum(subset) + 1

@inline presample(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = max(o, p, q)

Base.@propagate_inbounds @inline function update!(
            ht, lht, zt, at, ::Type{<:EGARCH{o, p ,q}}, garchcoefs
            ) where {o, p, q}
    mlht = garchcoefs[1]
    for i = 1:o
        mlht += garchcoefs[i+1]*zt[end-i+1]
    end
    for i = 1:p
        mlht += garchcoefs[i+1+o]*lht[end-i+1]
    end
    for i = 1:q
        mlht += garchcoefs[i+1+o+p]*(abs(zt[end-i+1]) - sqrt2invpi)
    end
    push!(lht, mlht)
    push!(ht, exp(mlht))
    return nothing
end

@inline function uncond(::Type{<:EGARCH{o, p, q}}, coefs::Vector{T}) where {o, p, q, T}
    eg = one(T)
    for i=1:max(o, q)
        γ = (i<=o ? coefs[1+i] : zero(T))
        α = (i<=q ? coefs[o+p+1+i] : zero(T))
        eg *= exp(-α*sqrt2invpi) * (exp(.5*(γ+α)^2)*normcdf(γ+α) + exp(.5*(γ-α)^2)*normcdf(α-γ))
    end
    h0 = (exp(coefs[1])*eg)^(1/(1-sum(coefs[o+2:o+p+1])))
end

function startingvals(spec::Type{<:EGARCH{o, p, q}}, data::Array{T}) where {o, p, q, T}
    x0 = zeros(T, o+p+q+1)
    x0[1]=1
    x0[2:o+1] .= 0
    x0[o+2:o+p+1] .= 0.9/p
    x0[o+p+2:end] .= 0.05/q
    x0[1] = var(data)/uncond(spec, x0)
    return x0
end

function startingvals(TT::Type{<:EGARCH}, data::Array{T} , subset::Tuple) where {T}
	o, p, q = subsettuple(TT, subsetmask(TT, subset)) # defend against (p, q) instead of (o, p, q)
	x0 = zeros(T, o+p+q+1)
    x0[2:o+1] .= 0.04/o
    x0[o+2:o+p+1] .= 0.9/p
    x0[o+p+2:end] .= o>0 ? 0.01/q : 0.05/q
    x0[1] = var(data)*(one(T)-sum(x0[2:o+1])/2-sum(x0[o+2:end]))
	mask = subsetmask(TT, subset)
	x0long = zeros(T, length(mask))
	x0long[mask] .= x0
    return x0long
end

function constraints(::Type{<:EGARCH{o, p,q}}, ::Type{T}) where {o, p, q, T}
    lower = zeros(T, o+p+q+1)
    upper = zeros(T, o+p+q+1)
    lower .=  T(-Inf)
    upper .= T(Inf)
    lower[1] = T(-Inf)
    lower[o+2:o+p+1] .= zero(T)
    upper[o+2:o+p+1] .= one(T)
    return lower, upper
end

function coefnames(::Type{<:EGARCH{o, p, q}}) where {o, p, q}
    names = Array{String, 1}(undef, o+p+q+1)
    names[1] = "ω"
    names[2:o+1] .= (i -> "γ"*subscript(i)).([1:o...])
    names[o+2:o+p+1] .= (i -> "β"*subscript(i)).([1:p...])
    names[o+p+2:o+p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
    return names
end

@inline function subsetmask(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subs) where {o, p, q, T}
	ind = falses(nparams(VS_large))
	subset = zeros(Int, 3)
	subset[4-length(subs):end] .= subs
	ind[1] = true
	os = subset[1]
	ps = subset[2]
	qs = subset[3]
	@assert os <= o
	@assert ps <= p
	@assert qs <= q
	ind[2:2+os-1] .= true
	ind[2+o:2+o+ps-1] .= true
	ind[2+o+p:2+o+p+qs-1] .= true
	ind
end

@inline function subsettuple(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subsetmask) where {o, p, q, T}
	os = 0
	ps = 0
	qs = 0
	@inbounds @simd ivdep for i = 2 : o + 1
		os += subsetmask[i]
	end
	@inbounds @simd ivdep for i = o + 2 : o + p + 1
		ps += subsetmask[i]
	end
	@inbounds @simd ivdep for i = o + p + 2 : o + p + q + 1
		qs += subsetmask[i]
	end
	(os, ps, qs)
end


text:
using BinDeps
using Compat

@BinDeps.setup

function validate_clp(name,handle)
    try
        p = Libdl.dlsym(handle, :Clp_VersionMajor)
        return p != C_NULL
    catch
        return false
    end
end

function validate_cbc(name,handle)
    try
        # Pre 2.8.12 doesn't have this defined
        p = Libdl.dlsym(handle, :Cbc_setInitialSolution)
        return p != C_NULL
    catch
        return false
    end
end

if is_unix()
    libclp = library_dependency("libclp",aliases=["libClp"], validate=validate_clp)
    libcbcsolver = library_dependency("libcbcsolver",aliases=["libCbcSolver"], validate=validate_cbc)
end
if is_windows()
    using WinRPM
    libclp = library_dependency("libclp",aliases=["libClp-1"], validate=validate_clp)
    libcbcsolver = library_dependency("libcbcsolver",aliases=["libCbcSolver-3"], validate=validate_cbc)
    provides(WinRPM.RPM, "Cbc", [libclp,libcbcsolver], os = :Windows)
end

cbcname = "Cbc-2.9.8"

provides(Sources, URI("http://www.coin-or.org/download/source/Cbc/$cbcname.tgz"),
    [libclp,libcbcsolver], os = :Unix)

if is_apple()
    using Homebrew
    if Homebrew.installed("coinmp") # coinmp package is old and conflicts with cbc package
        Homebrew.rm("coinmp")
    end
    provides( Homebrew.HB, "cbc", [libclp, libcbcsolver], os = :Darwin )
end

prefix=joinpath(BinDeps.depsdir(libclp),"usr")
patchdir=BinDeps.depsdir(libclp)
srcdir = joinpath(BinDeps.depsdir(libclp),"src",cbcname)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libclp)
        @build_steps begin
            ChangeDirectory(srcdir)
            `./configure --prefix=$prefix --without-blas --without-lapack --enable-cbc-parallel`
            `make install`
        end
    end),[libclp,libcbcsolver], os = :Unix)

@BinDeps.install Dict(:libclp => :libclp, :libcbcsolver => :libcbcsolver)


text:
# Starting of with Graph.jl implementation. There are a lot of overlap with metagraphs and some data can be moved there.
"""
$(TYPEDEF)
Encapsulation structure for a DFGNode (Variable or Factor) in LightGraphs.jl graph.
"""
mutable struct LightGraphsNode
    index::Int
    dfgNode::DFGNode
end

const LFGType = MetaGraph{Int64,Float64}

mutable struct MetaGraphsDFG{T <: AbstractParams} <: AbstractDFG
    g::LFGType
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    #NOTE Removed nodeCounter
    # nodeCounter::Int64
    #NOTE using matagraphs labels
    # labelDict::Dict{Symbol, Int64}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
end

#TODO? do we not want props such as userId, robotId, sessionId, etc...
function MetaGraphsDFG{T}(g::LFGType=MetaGraph();
                           description::String="LightGraphs.jl implementation",
                           userId::String="User ID",
                           robotId::String="Robot ID",
                           sessionId::String="Session ID",
                           userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           params::T=NoSolverParams()) where T <: AbstractParams
    set_prop!(g, :description, description)
    set_prop!(g, :userId, userId)
    set_prop!(g, :robotId, robotId)
    set_prop!(g, :sessionId, sessionId)
    set_indexing_prop!(g, :label)
    MetaGraphsDFG{T}(g, description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], params)
end

Base.propertynames(x::MetaGraphsDFG, private::Bool=false) =
    (:g, :description, :userId, :robotId, :sessionId, :nodeCounter, :labelDict, :addHistory, :solverParams)
        # (private ? fieldnames(typeof(x)) : ())...)

Base.getproperty(x::MetaGraphsDFG,f::Symbol) = begin
    if f == :nodeCounter
        @warn "Field nodeCounter depreciated. returning number of nodes"
        nv(x.g)
    elseif f == :labelDict
        @warn "Field labelDict depreciated. Consider using exists(dfg,label) or getLabelDict(dfg) instead. Returning internals copy"
        copy(x.g.metaindex[:label])
    else
        getfield(x,f)
    end
end

"""
$(TYPEDEF)
Basic Structure for just the LightGraphs object and keys
"""
struct LightGraphsSkeleton
    G::LightGraphs.SimpleGraph
    labels::Dict{Symbol, Int}
end

LightGraphsSkeleton(dfg::MetaGraphsDFG) = LightGraphsSkeleton(dfg.g.graph, dfg.g.metaindex[:label])


text:
using SkewDist
using Base.Test

function rand_pd(d::Int)
    A = randn(d,d)
    A'A
end

function rand_cor(d::Int)
    A = randn(d,d)
    Σ = A'A
    cov2cor(Σ)
end

include("test_MvSkewNormal.jl")
include("test_MvSkewTDist.jl")
include("test_operators.jl")
include("test_fit_skew.jl")
include("test_skewtdist.jl")


text:
# # constructor
_BlockDiagonal_pullback(Δ::Tangent, nrows, ncols) = (NoTangent(), Δ.blocks)
_BlockDiagonal_pullback(Δ::BlockDiagonal, nrows, ncols) = (NoTangent(), Δ.blocks)
function _BlockDiagonal_pullback(Δ::AbstractThunk, nrows, ncols)
    _BlockDiagonal_pullback(unthunk(Δ), nrows, ncols)
end
function _BlockDiagonal_pullback(Δ::AbstractArray, nrows, ncols)
    row_idxs = cumsum(nrows) .- nrows .+ 1
    col_idxs = cumsum(ncols) .- ncols .+ 1
    Δblocks = map(eachindex(nrows)) do n
        block_rows = row_idxs[n]:(row_idxs[n] + nrows[n] - 1)
        block_cols = col_idxs[n]:(col_idxs[n] + ncols[n] - 1)
        return Δ[block_rows, block_cols]
    end
    return NoTangent(), Δblocks
end

function ChainRulesCore.rrule(::Type{<:BlockDiagonal}, blocks::Vector{V}) where {V}
    nrows = size.(blocks, 1)
    ncols = size.(blocks, 2)
    function BlockDiagonal_pullback(Δ)
        return _BlockDiagonal_pullback(Δ, nrows, ncols)
    end
    return BlockDiagonal(blocks), BlockDiagonal_pullback
end

# densification
function _densification_pullback(Ȳ::Matrix, T, nrows, ncols)
    row_idxs = cumsum(nrows) .- nrows .+ 1
    col_idxs = cumsum(ncols) .- ncols .+ 1

    Δblocks = map(eachindex(nrows)) do n
        block_rows = row_idxs[n]:(row_idxs[n] + nrows[n] - 1)
        block_cols = col_idxs[n]:(col_idxs[n] + ncols[n] - 1)
        return Ȳ[block_rows, block_cols]
    end
    return (NoTangent(), Tangent{T}(blocks=Δblocks))
end
function _densification_pullback(Ȳ::AbstractThunk, T, nrows, ncols)
    return _densification_pullback(unthunk(Ȳ), T, nrows, ncols)
end
function ChainRulesCore.rrule(::Type{<:Base.Matrix}, B::T) where {T<:BlockDiagonal}
    nrows = size.(B.blocks, 1)
    ncols = size.(B.blocks, 2)
    densification_pullback(ȳ) = _densification_pullback(ȳ, T, nrows, ncols)
    return Matrix(B), densification_pullback
end

# multiplication
function ChainRulesCore.rrule(
        ::typeof(*),
        bm::BlockDiagonal{T, V},
        v::StridedVector{T}
    ) where {T<:Union{Real, Complex}, V<:Matrix{T}}

    y = bm * v

    # needed for computing Δ * v' blockwise
    nrows = size.(bm.blocks, 1)
    ncols = size.(bm.blocks, 2)
    row_idxs = cumsum(nrows) .- nrows .+ 1
    col_idxs = cumsum(ncols) .- ncols .+ 1

    function bm_vector_mul_pullback(Δy)
        ȳ = unthunk(Δy)
        Δblocks = map(eachindex(nrows)) do i
            block_rows = row_idxs[i]:(row_idxs[i] + nrows[i] - 1)
            block_cols = col_idxs[i]:(col_idxs[i] + ncols[i] - 1)
            return InplaceableThunk(
                X̄ -> mul!(X̄, ȳ[block_rows], v[block_cols]', true, true),
                @thunk(ȳ[block_rows] * v[block_cols]'),
            )
        end

        b̄m = Tangent{BlockDiagonal{T, V}}(;blocks=Δblocks)
        v̄ = InplaceableThunk(X̄ -> mul!(X̄, bm', ȳ, true, true), @thunk(bm' * ȳ))
        return NoTangent(), b̄m, v̄
    end
    return y, bm_vector_mul_pullback
end

function ProjectTo(b::BlockDiagonal)
    blocks = map(ProjectTo, b.blocks)
    return ProjectTo{BlockDiagonal}(; blocks=blocks, blocksizes=blocksizes(b))
end

function (project::ProjectTo{BlockDiagonal})(dx::AbstractArray)
    # prepare to index into the dense array
    nrows = first.(project.blocksizes)
    ncols = last.(project.blocksizes)
    row_idxs = cumsum(nrows) .- nrows .+ 1
    col_idxs = cumsum(ncols) .- ncols .+ 1
    # project each block individually
    blocks = map(eachindex(nrows)) do i
        block_rows = row_idxs[i]:(row_idxs[i] + nrows[i] - 1)
        block_cols = col_idxs[i]:(col_idxs[i] + ncols[i] - 1)
        project.blocks[i](dx[block_rows, block_cols])
    end
    return BlockDiagonal(blocks)
end


text:
@testset "block_node" begin
    # We should identity-map wrapping of the same block.
    x = block_node(b1)
    y = block_node(b1)
    @test x.op == y.op
    @test x == y
    @test x === y

    # Similar blocks that are non-identical should be identity mapped.
    x = block_node(Block{Int64}())
    y = block_node(Block{Int64}())
    @test x.op == y.op
    @test x == y
    @test x === y
end

@testset "iterdates" begin
    @testset "midnight" begin
        n = iterdates()
        @test n === iterdates(Time(0))
        @test value_type(n) == DateTime

        expected_times = collect(_T_START:Day(1):(_T_END - Day(1)))
        @test _eval(n) == Block(expected_times, expected_times)
    end

    @testset "non-midnight" begin
        n = iterdates(Time(1))
        @test n === iterdates(Time(1))
        @test value_type(n) == DateTime

        first = Date(_T_START) + Time(1)
        last = Date(_T_END) - Day(1) + Time(1)
        expected_times = collect(first:Day(1):last)
        @test _eval(n) == Block(expected_times, expected_times)
    end

    @testset "non-UTC" begin
        timezone = tz"America/Chicago"

        n = iterdates(Time(0), timezone)
        @test n === iterdates(Time(0), timezone)
        @test value_type(n) == DateTime

        # Evaluate over a sufficient range to get some DST changes.
        block = _evaluate(n, DateTime(2020), DateTime(2021))
        @test block.times === block.values
        zdts = ZonedDateTime.(block.times, tz"UTC")
        local_time_of_day = Time.(astimezone.(zdts, timezone))
        @test all(local_time_of_day .== Time(0))
    end

    @testset "non-existent time" begin
        timezone = tz"America/Chicago"
        n = iterdates(Time(2), timezone)
        @test_throws NonExistentTimeError _evaluate(n, DateTime(2020), DateTime(2021))
    end
end

@testset "pulse" begin
    @testset "invalid" begin
        @test_throws ArgumentError pulse(Hour(0))
        @test_throws ArgumentError pulse(Hour(-1))
    end

    @testset "1 day" begin
        n = pulse(Hour(24))
        expected_times = collect(_T_START:Day(1):(_T_END - Day(1)))
        @test _eval(n) == Block(expected_times, expected_times)
    end

    @testset "arbitrary" begin
        # By default, our pulses should be aligned with the default epoch, which is the same
        # as Julia's internal epoch.
        epoch = DateTime(0, 1, 1) - Millisecond(Dates.DATETIMEEPOCH)
        @test epoch.instant.periods == Millisecond(0)

        for delta in (Second(3), Minute(10), Hour(2))
            n = pulse(delta)

            # We expect to be aligned to the Julia epoch. Work in terms of milliseconds
            # since this epoch.
            start = _T_START.instant.periods
            remainder = start % Millisecond(delta)

            pulse_start = remainder == Millisecond(0) ? start : start + (delta - remainder)
            expected_times = [epoch + pulse_start]
            while true
                prev_time = last(expected_times)
                next_time = prev_time + delta
                next_time < _T_END || break
                push!(expected_times, next_time)
            end

            @test _eval(n) == Block(expected_times, expected_times)
        end
    end

    @testset "custom epoch" begin
        delta = Minute(10)
        for offset in Minute(0):Minute(30):Minute(Hour(24))
            epoch = _T_START + offset
            n = pulse(delta; epoch)

            rem_ = (_T_START - epoch) % Millisecond(delta)
            pulse_start = rem_ == Millisecond(0) ? _T_START : _T_START + (delta - rem_)
            expected_times = [pulse_start]
            while true
                prev_time = last(expected_times)
                next_time = prev_time + delta
                next_time < _T_END || break
                push!(expected_times, next_time)
            end
        end
    end

    @testset "equivalents" begin
        @test pulse(Hour(1)) === pulse(Hour(1))
        @test pulse(Hour(1)) === pulse(Minute(60))
        @test pulse(Hour(1)) === pulse(Second(60 * 60))
        @test pulse(Hour(1)) === pulse(Millisecond(60 * 60 * 1000))
        @test pulse(Hour(1); epoch=DateTime(2020)) !=
            pulse(Hour(1); epoch=DateTime(2020, 1, 1, 1, 30))
        # Optimising this in the general case is too hard, even though these two nodes
        # ought to be equivalent
        @test pulse(Hour(1); epoch=DateTime(2020)) != pulse(Hour(1); epoch=DateTime(2021))
    end
end

@testset "tea_file" begin
    mktempdir() do prefix
        # Write some basic data to a file, then read it back.
        path = joinpath(prefix, "moo.tea")
        TeaFiles.write(path, b1)

        n = TimeDag.tea_file(path, :value)
        n2 = TimeDag.tea_file(path, :value)

        @test n2 === n
        @test _eval(n) == b1
    end
end


text:
# Generate some kind of plot to use for the banner image

using ImageProjectiveGeometry
using PyPlot

figure(1)
clf()

# Generate 100 points
nGpts = 100
dx = 700; dy = 700; dz = 100

gx = (2*rand(nGpts)-1) * dx
gy = (2*rand(nGpts)-1) * dy
gz = (2*rand(nGpts)-1) * dz

plot3D(gx,gy,gz,"k*")
hold(true)

# Generate Cameras
f = 4000
rows = 3000;  cols = 4000
ppx = cols/2; ppy = rows/2;
Z = 1000

Ca = []
for X = -500:500:500
    for Y = -800:400:800
        Rc_w = rotx(pi)[1:3,1:3]     # camera points down
        push!(Ca, Camera(P=[X, Y, Z], Rc_w=Rc_w, fx=f, fy=f, ppx=ppx, ppy=ppy,
                         rows=rows, cols=cols))
    end
end

plotcamera(Ca,150)
axis("equal")


text:
using JuMP, VariationalInequality

# https://cdr.lib.unc.edu/indexablecontent/uuid:778ca632-74ca-4858-8c3c-6dcfc7e6e703
# Example 4.2. This example is adapted from the example in (Fukushima, 1986). Let
m = VIPModel()
@variable(m, x[1:3])

@constraint(m, x[1]^2 + 0.4x[2]^2 + 0.6x[3]^2 <= 1)
@constraint(m, 0.6x[1]^2 + 0.4x[2]^2 + x[3]^2 <= 1)
@constraint(m, x[1] + x[2] + x[3] >= sqrt(3))

@mapping(m, F1, 2x[1] + 0.2x[1]^3 - 0.5x[2] + 0.1x[3] - 4)
@mapping(m, F2, -0.5x[1] + x[2] + 0.1x[2]^3 + 0.5)
@mapping(m, F3, 0.5x[1] - 0.2x[2] + 2x[3] - 0.5)

@innerproduct(m, F1, x[1])
@innerproduct(m, F2, x[2])
@innerproduct(m, F3, x[3])
# or
# @innerproduct(m, [F1, F2, F3], x)

solveVIP(m, algorithm=:fixed_point, max_iter=1000, step_size=0.1)
sol1, Fval1, gap1 = save_solution(m)
@assert 0<= gap1 < 1e-6

# The above `solveVIP` sets the value of variables at the solution
clear_values(m)

solveVIP(m, algorithm=:extra_gradient, max_iter=1000, step_size=0.1)
sol2, Fval2, gap2 = save_solution(m)
@assert 0<= gap2 < 1e-6

# sol2, Fval2, gap2 = solveVIP(m, algorithm=:hyperplane, max_iter=1000, step_size=0.1)

@show sol1
@show sol2
@show gap1
@show gap2

# OKAY.
# x∗ = (0.9168, 0.4850, 0.3303)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ImageReference <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    publisher::Any # spec type: Union{ Nothing, String } # spec name: publisher
    offer::Any # spec type: Union{ Nothing, String } # spec name: offer
    sku::Any # spec type: Union{ Nothing, String } # spec name: sku
    version::Any # spec type: Union{ Nothing, String } # spec name: version

    function ImageReference(;id=nothing, publisher=nothing, offer=nothing, sku=nothing, version=nothing)
        o = new()
        validate_property(ImageReference, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ImageReference, Symbol("publisher"), publisher)
        setfield!(o, Symbol("publisher"), publisher)
        validate_property(ImageReference, Symbol("offer"), offer)
        setfield!(o, Symbol("offer"), offer)
        validate_property(ImageReference, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(ImageReference, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        o
    end
end # type ImageReference

const _property_map_ImageReference = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("publisher")=>Symbol("publisher"), Symbol("offer")=>Symbol("offer"), Symbol("sku")=>Symbol("sku"), Symbol("version")=>Symbol("version"))
const _property_types_ImageReference = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("publisher")=>"String", Symbol("offer")=>"String", Symbol("sku")=>"String", Symbol("version")=>"String")
Base.propertynames(::Type{ ImageReference }) = collect(keys(_property_map_ImageReference))
Swagger.property_type(::Type{ ImageReference }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ImageReference[name]))}
Swagger.field_name(::Type{ ImageReference }, property_name::Symbol) =  _property_map_ImageReference[property_name]

function check_required(o::ImageReference)
    true
end

function validate_property(::Type{ ImageReference }, name::Symbol, val)
end


text:
# spherical quadrature
KitBase.legendre_quadrature(6)
KitBase.octa_quadrature(8)

KitBase.VSpace1D() |> show
KitBase.VSpace2D() |> show
KitBase.VSpace3D() |> show
KitBase.MVSpace1D() |> show
KitBase.MVSpace2D() |> show
KitBase.MVSpace3D() |> show

KitBase.VSpace1D(-5, 5, 16, "newton")
KitBase.VSpace1D(-5, 5, 16, "algebra")
KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "newton")
KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "algebra")
KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "maxwell")
KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, "newton")
KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, "algebra")
KitBase.MVSpace1D(-5, 5, -5, 5, 16, "newton")
KitBase.MVSpace2D(-5, 5, -5, 5, 16, -5, 5, -5, 5, 16, "newton")
KitBase.MVSpace3D(-5, 5, -5, 5, 8, -5, 5, -5, 5, 8, -5, 5, -5, 5, 8, "newton")
KitBase.UnstructVSpace(-1, 1, 16, rand(16), ones(16))


text:
import Base.TCPServer

type PokeServer
    socket::TCPServer
    cache::SortedDict{Vector{UInt8}, Vector{UInt8}}
    cache_size::Int64
    dump_path::UTF8String
end

function startserver(port::Int, dump_path::AbstractString)
    server_socket = listen(port)
    pserv = PokeServer(server_socket, new_cache(), 0, dump_path)
    task = @async begin
        try
            while true
                sock = accept(server_socket)
                serve_conn(pserv, sock)
            end
        catch ex
            if !isa(ex, InterruptException)
                throw(e)
            end
        finally
            close(server_socket)
        end
    end
    return task, pserv  # TODO: return pserv instead
end

function stopserver(server_task::Task)
    Base.throwto(server_task, InterruptException())
end

function serve_conn(poke_server::PokeServer, sock::TCPSocket)
    while true
        req = readobj(sock, PokeRequest)
        println(req)
        for (k, v) in req.data
            cache[k] = v
        end
        # TODO: check if size of cache is large enough to start merge
        write(RequestResponse(req.id, 0, ""))
    end
end

function Base.isless(x::Vector{UInt8}, y::Vector{UInt8})
    return isless(bytestring(x), bytestring(y))
end

function merge()

end

function mergedump(dump_path::AbstractString, cache::Associative)
    if isfile(dump_path)
        # create

    else
        # merge
        dumpf = open(dump_path)
        temp_path, tempf = mktemp()
        try

        finally
            close(dumpf)
            close(tempf)
            rm(tempf)
        end
    end
end


text:
# Restricted for a single non zero block
struct BlockFieldArrayCoo{T,N,A,X} <: AbstractBlockArray{T,N}
  axes::X
  blockids::Vector{NTuple{N,Int}}
  block::A
  function BlockFieldArrayCoo(
    _axes::NTuple{N},
    blockids::Vector{NTuple{N,Int}},
    block::A) where {T,N,A<:AbstractArray{T,N}}

    @assert length(blockids) == 1
    #I = first(blockids)
    #@check blocks_equal(axes(block),map(local_range,_axes,I)) "The given block and axes are incompatible."

    X = typeof(_axes)
    new{T,N,A,X}(_axes,blockids,block)
  end
end

struct BlockFieldArrayCooMap{N} <: Map
  blocksize::NTuple{N,Int}
  blockids::Vector{NTuple{N,Int}}
  ptrs::Array{Int,N}
  function BlockFieldArrayCooMap(blocksize::NTuple{N,Int}, blockids::Vector{NTuple{N,Int}}) where N
    ptrs = fill(-1,blocksize)
    for (p,I) in enumerate(blockids)
      ptrs[I...] = p
      for i in 1:N
        @check blocksize[i] >= I[i]
      end
    end
    new{N}(blocksize,blockids,ptrs)
  end
end

@inline function lazy_map(k::BlockFieldArrayCooMap,T::Type,f::AbstractArray...)
  s = Arrays._common_size(f...)
  N = length(s)
  LazyArray(T,Val(N),Fill(k,s),f...)
end

BlockArrays.blocksize(a::BlockFieldArrayCooMap) = a.blocksize
is_zero_block(a::BlockFieldArrayCooMap,i::Integer) = a.ptrs[i] < 0
is_zero_block(a::BlockFieldArrayCooMap{N},i::Vararg{Integer,N}) where N = a.ptrs[i...] < 0
is_zero_block(a::BlockFieldArrayCooMap{N},i::Vararg{Block,N}) where N = is_zero_block(a,map(Int,i)...)
is_zero_block(a::BlockFieldArrayCooMap,i::Block) = is_zero_block(a,convert(Tuple,i)...)
is_zero_block(a::BlockFieldArrayCooMap,i::CartesianIndex) = is_zero_block(a,Tuple(i)...)

function evaluate!(cache,k::BlockFieldArrayCooMap,axes,block)
  @check map(i->first(blocksize(i)),axes) == k.blocksize "The given axes are not compatible with the given BlockFieldArrayCooMap"
  BlockFieldArrayCoo(axes,k.blockids,block)
end

testitem(a::BlockFieldArrayCoo) = testitem(a.block)

# Specific API

function is_zero_block(a::BlockFieldArrayCoo{T,N},i::Vararg{Integer,N}) where {T,N}
  i != first(a.blockids)
end

# AbstractBlockArray

@inline function BlockArrays.getblock(a::BlockFieldArrayCoo{T,N}, i::Vararg{Integer, N}) where {T,N}
  if i == first(a.blockids)
    a.block
  else
    @notimplemented "Cannot get a zero block from a BlockFieldArrayCoo"
  end
end

# AbstractArray

Base.size(a::BlockFieldArrayCoo) = map(length,Base.axes(a))
Base.axes(a::BlockFieldArrayCoo) = a.axes
Base.IndexStyle(::Type{<:BlockFieldArrayCoo}) = IndexCartesian()

function Base.getindex(a::BlockFieldArrayCoo{T,N},i::Vararg{Integer,N}) where {T,N}
  s = map(findblockindex,a.axes,i)
  ai = a[s...]
  ai
end

# Evaluation

function return_cache(a::BlockFieldArrayCoo,x::Point)
  fc = return_cache(a.block,x)
  fx = return_value(a.block,x)
  bsize = map(i->first(blocksize(i)),a.axes)
  k = BlockArrayCooMap(bsize,a.blockids)
  cr = return_cache(k,a.axes,fx)
  (k,fc,cr)
end

@inline function evaluate!(cache,a::BlockFieldArrayCoo,x::Point)
  k,fc,cr = cache
  fx = evaluate!(fc,a.block,x)
  evaluate!(cr,k,a.axes,fx)
end

function return_cache(a::BlockFieldArrayCoo,x::AbstractVector{<:Point})
  fc = return_cache(a.block,x)
  fx = return_value(a.block,x)
  pr = similar_range(first(a.axes),length(x))
  axs = (pr,a.axes...)
  blockids = map(i->(1,i...),a.blockids)
  bsize = map(i->first(blocksize(i)),axs)
  k = BlockArrayCooMap(bsize,blockids)
  cr = return_cache(k,axs,fx)
  (k,fc,axs,cr)
end

@inline function evaluate!(cache,a::BlockFieldArrayCoo,x::AbstractVector{<:Point})
  k,fc,axs,cr = cache
  fx = evaluate!(fc,a.block,x)
  evaluate!(cr,k,axs,fx)
end

# Gradient

function evaluate!(cache,k::Broadcasting{typeof(∇)},a::BlockFieldArrayCoo)
  g = k(a.block)
  BlockFieldArrayCoo(a.axes,a.blockids,g)
end

function evaluate!(cache,k::Broadcasting{typeof(∇∇)},a::BlockFieldArrayCoo)
  g = k(a.block)
  BlockFieldArrayCoo(a.axes,a.blockids,g)
end

# Transpose

function Base.transpose(a::BlockFieldArrayCoo{T,1} where T)
  r = similar_range(first(axes(a)),1)
  axs = (r,axes(a)...)
  blockids = map(i->(1,i...),a.blockids)
  BlockFieldArrayCoo(axs,blockids,transpose(a.block))
end

# Global optimizations

# Evaluation

function lazy_map(
  ::typeof(evaluate),
  a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}},
  x::AbstractArray{<:Point})

  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_aix = lazy_map(evaluate,cell_ai,x)
  lazy_map(BlockArrayCooMap(m.blocksize,m.blockids),cell_axs,cell_aix)
end

function lazy_map(
  ::typeof(evaluate),
  a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}},
  cell_x::AbstractArray{<:AbstractVector{<:Point}})

  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_aix = lazy_map(evaluate,cell_ai,cell_x)

  function result_axes_on_point_vector(axs,x)
    pr = similar_range(first(axs),length(x))
    axs = (pr,axs...)
  end

  cell_axs_new = lazy_map(result_axes_on_point_vector,cell_axs,cell_x)

  bsize_new = (1,m.blocksize...)
  bids_new = map(i->(1,i...),m.blockids)

  r = lazy_map(BlockArrayCooMap(bsize_new,bids_new),cell_axs_new, cell_aix)
end

# Gradient

function lazy_map(
  k::Broadcasting{typeof(∇)}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}})

  m = a.maps.value
  cell_axs, cell_ai = a.args

  cell_gi = lazy_map(k,cell_ai)
  lazy_map(m,cell_axs, cell_gi)
end

function lazy_map(
  k::Broadcasting{typeof(∇∇)}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}})

  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_gi = lazy_map(k,cell_ai)
  lazy_map(m,cell_axs,cell_gi)
end

function lazy_map(::typeof(axes),a::LazyArray{<:Fill{Broadcasting{typeof(∇)}}})
  b = a.args[1]
  lazy_map(axes,b)
end

function lazy_map(::typeof(axes),a::LazyArray{<:Fill{Broadcasting{typeof(∇∇)}}})
  b = a.args[1]
  lazy_map(axes,b)
end

# Composition

function lazy_map(
  k::Broadcasting{typeof(∘)}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}},b::AbstractArray{<:Field})

  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_gi = lazy_map(k,cell_ai,b)
  lazy_map(m,cell_axs,cell_gi)
end

function lazy_map(::typeof(axes),a::LazyArray{<:Fill{Broadcasting{typeof(∘)}}})
  b = a.args[1]
  lazy_map(axes,b)
end

# Transpose

function lazy_map(
  k::typeof(transpose), a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}})

  m = a.maps.value
  cell_axs, cell_ai = a.args

  cell_ait = lazy_map(transpose,cell_ai)

  cell_axs_new = lazy_map(cell_axs) do axs
    @check length(axs) == 1
    r = similar_range(first(axs),1)
    (r,axs...)
  end

  bsize_new = (1,m.blocksize...)
  bids_new = map(i->(1,i...),m.blockids)

  lazy_map(BlockFieldArrayCooMap(bsize_new,bids_new),cell_axs_new, cell_ait)
end

# Operations before evaluating

function lazy_map(
  k::Broadcasting{<:Operation}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}})
  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_aif = lazy_map(k,cell_ai)
  lazy_map(m,cell_axs,cell_aif)
end

function lazy_map(
  k::Broadcasting{<:Operation}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}},f::AbstractArray{<:Field})
  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_aif = lazy_map(k,cell_ai,f)
  lazy_map(m,cell_axs,cell_aif)
end

function lazy_map(
  k::Broadcasting{<:Operation}, f::AbstractArray{<:Field}, a::LazyArray{<:Fill{<:BlockFieldArrayCooMap}})
  m = a.maps.value
  cell_axs, cell_ai = a.args
  cell_aif = lazy_map(k,f,cell_ai)
  lazy_map(m,cell_axs,cell_aif)
end

# Operations on values

function _get_axes_and_blocks(f)
  f[1], f[2:end]
end
#
# Unary operations
# Assumption: op is a scaling of a
function lazy_map(
  k::BroadcastingFieldOpMap, a::LazyArray{<:Fill{<:BlockArrayCooMap}})
  m = a.maps.value
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  cell_blocks_new = map(b->lazy_map(k,b),cell_blocks)
  lazy_map(m,cell_axs,cell_blocks_new...)
end

# Binary test/field or trial/field

# Assumption: op is a scaling of a
function _op_basis_vs_field(k,a,f)
  @check ndims(eltype(a)) > ndims(eltype(f))
  m = a.maps.value
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  cell_blocks_new = map(b->lazy_map(k,b,f),cell_blocks)
  lazy_map(m,cell_axs,cell_blocks_new...)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{1}}},
  f::AbstractArray{<:Number})
  _op_basis_vs_field(k,a,f)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{2}}},
  f::AbstractArray{<:Number})
  _op_basis_vs_field(k,a,f)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{2}}},
  f::AbstractArray{<:AbstractVector})
  _op_basis_vs_field(k,a,f)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{3}}},
  f::AbstractArray{<:AbstractVector})
  _op_basis_vs_field(k,a,f)
end

# Assumption: op is a scaling of a
function _op_field_vs_basis(k,f,a)
  @check ndims(eltype(a)) > ndims(eltype(f))
  m = a.maps.value
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  cell_blocks_new = map(b->lazy_map(k,f,b),cell_blocks)
  lazy_map(m,cell_axs,cell_blocks_new...)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  f::AbstractArray{<:Number},
  a::LazyArray{<:Fill{BlockArrayCooMap{1}}})
  _op_field_vs_basis(k,f,a)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  f::AbstractArray{<:Number},
  a::LazyArray{<:Fill{BlockArrayCooMap{2}}})
  _op_field_vs_basis(k,f,a)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  f::AbstractArray{<:AbstractVector},
  a::LazyArray{<:Fill{BlockArrayCooMap{2}}})
  _op_field_vs_basis(k,f,a)
end

function lazy_map(
  k::BroadcastingFieldOpMap,
  f::AbstractArray{<:AbstractVector},
  a::LazyArray{<:Fill{BlockArrayCooMap{3}}})
  _op_field_vs_basis(k,f,a)
end

# Binary test/test or trial/trial
# Assumption: op is a linear combination of a and b
function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{N}}},
  b::LazyArray{<:Fill{BlockArrayCooMap{N}}}) where N

  ma = a.maps.value
  mb = b.maps.value
  @assert blocksize(ma) == blocksize(mb)

  cell_axs, = a.args

  blocks = []
  blockids = eltype(ma.blockids)[]
  for I in eachblockid(ma)
    if is_nonzero_block(ma,I) || is_nonzero_block(mb,I)
      if is_nonzero_block(ma,I) && is_nonzero_block(mb,I)
        aI = _get_cell_block(ma,a,I)
        bI = _get_cell_block(mb,b,I)
      elseif is_nonzero_block(ma,I)
        aI = _get_cell_block(ma,a,I)
        bI = _lazy_array_of_zero_blocks(_block_type(eltype(b)),cell_axs,I,aI)
      else
        bI = _get_cell_block(mb,b,I)
        aI = _lazy_array_of_zero_blocks(_block_type(eltype(a)),cell_axs,I,bI)
      end
      block = lazy_map(k,aI,bI)
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end

  @assert length(blocks) > 0
  m = BlockArrayCooMap(blocksize(ma),blockids)
  lazy_map(m,cell_axs,blocks...)
end

function _block_type(::Type{<:BlockArrayCoo{T,N,A}}) where {T,N,A}
  A
end

function _get_cell_block(
  m::BlockArrayCooMap,
  a::AbstractArray{<:BlockArrayCoo{T,N,A}},
  b::Block) where {T,N,A}
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  I = convert(Tuple,b)
  p = m.ptrs[I...]
  @assert p > 0
  cell_blocks[p]
end

function _lazy_array_of_zero_blocks(::Type{A},cell_axs,I,aI) where A
  lazy_map(cell_axs) do axs
    laxs = map( local_range, axs, I.n)
    Arrays._zero_block(A,laxs)
  end
end

function _lazy_array_of_zero_blocks(
  ::Type{A},cell_axs,I,aI::LazyArray{<:Fill{BlockArrayCooMap{N}}}) where {A,N}

  cell_axsI, = aI.args
  maI = aI.maps.value
  blocks = []
  for J in eachblockid(maI)
    if is_nonzero_block(maI,J)
      block = _lazy_array_of_zero_blocks(_block_type(A),cell_axsI,J,nothing)
      push!(blocks,block)
    end
  end

  lazy_map(maI,cell_axsI,blocks...)
end

# Binary test/trial
# Assumption: op is a product of a and b

#function lazy_map(
#  k::BroadcastingFieldOpMap,
#  a::LazyArray{<:Fill{BlockArrayCooMap{1}}},
#  b::LazyArray{<:Fill{BlockArrayCooMap{2}}})
#
#end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{2}}},
  b::LazyArray{<:Fill{BlockArrayCooMap{3}}})

  cell_axs_a, = a.args
  cell_axs_b, = b.args
  cell_axs = lazy_map((a1,a2) -> (a1[1],a1[2],a2[3]),cell_axs_a,cell_axs_b)
  ma = a.maps.value
  mb = b.maps.value

  blocks = []
  blockids = NTuple{3,Int}[]
  nfield1 = size(ma.ptrs,2)
  nfield2 = size(mb.ptrs,3)
  for f1 in 1:nfield1
    I1 = Block(1,f1)
    for f2 in 1:nfield2
      I2 = Block(1,1,f2)
      if is_nonzero_block(ma,I1) && is_nonzero_block(mb,I2)
        aI1 = _get_cell_block(ma,a,I1)
        bI2 = _get_cell_block(mb,b,I2)
        block = lazy_map(k,aI1,bI2)
        push!(blocks,block)
        push!(blockids,(1,f1,f2))
      end
    end
  end

  @assert length(blocks) > 0
  bs = (ma.blocksize[1],ma.blocksize[2],mb.blocksize[3])
  m = BlockArrayCooMap(bs,blockids)
  lazy_map(m,cell_axs,blocks...)
end

#function lazy_map(
#  k::BroadcastingFieldOpMap,
#  a::LazyArray{<:Fill{BlockArrayCooMap{2}}},
#  b::LazyArray{<:Fill{BlockArrayCooMap{1}}})
#
#end

function lazy_map(
  k::BroadcastingFieldOpMap,
  a::LazyArray{<:Fill{BlockArrayCooMap{3}}},
  b::LazyArray{<:Fill{BlockArrayCooMap{2}}})

  cell_axs_a, = a.args
  cell_axs_b, = b.args
  cell_axs = lazy_map((a1,a2) -> (a2[1],a2[2],a1[3]),cell_axs_a,cell_axs_b)
  ma = a.maps.value
  mb = b.maps.value

  blocks = []
  blockids = NTuple{3,Int}[]
  nfield1 = size(mb.ptrs,2)
  nfield2 = size(ma.ptrs,3)
  for f1 in 1:nfield1
    I1 = Block(1,f1)
    for f2 in 1:nfield2
      I2 = Block(1,1,f2)
      if is_nonzero_block(mb,I1) && is_nonzero_block(ma,I2)
        bI1 = _get_cell_block(mb,b,I1)
        aI2 = _get_cell_block(ma,a,I2)
        block = lazy_map(k,bI1,aI2)
        push!(blocks,block)
        push!(blockids,(1,f1,f2))
      end
    end
  end

  @assert length(blocks) > 0
  bs = (mb.blocksize[1],mb.blocksize[2],ma.blocksize[3])
  m = BlockArrayCooMap(bs,blockids)
  lazy_map(m,cell_axs,blocks...)
end

# Integration of elem vectors
function lazy_map(k::IntegrationMap,a::LazyArray{<:Fill{BlockArrayCooMap{2}}},w::AbstractArray,j::AbstractArray)
  ma = a.maps.value
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  cell_axs_new = lazy_map(a->(a[2],),cell_axs)
  blocks = map(block->lazy_map(k,block,w,j),cell_blocks)
  blockids = [ (ids[2],) for ids in ma.blockids ]
  m = BlockArrayCooMap(ma.blocksize[2:end],blockids)
  lazy_map(m,cell_axs_new,blocks...)
end

## Integration of elem matrices
function lazy_map(k::IntegrationMap,a::LazyArray{<:Fill{BlockArrayCooMap{3}}},w::AbstractArray,j::AbstractArray)
  ma = a.maps.value
  cell_axs, cell_blocks = _get_axes_and_blocks(a.args)
  cell_axs_new = lazy_map(a->(a[2],a[3]),cell_axs)
  blocks = map(block->lazy_map(k,block,w,j),cell_blocks)
  blockids = [ (ids[2],ids[3]) for ids in ma.blockids ]
  m = BlockArrayCooMap(ma.blocksize[2:end],blockids)
  lazy_map(m,cell_axs_new,blocks...)
end
#function lazy_map(k::IntMap,f::VectorOfBlockArrayCoo{T,3} where T,w::AbstractArray,j::AbstractArray)
#  ax = lazy_map(a->(a[2],a[3]),f.axes)
#  blocks = map(block->lazy_map(k,block,w,j),f.blocks)
#  blockids = [ (ids[2], ids[3]) for ids in f.blockids ]
#  VectorOfBlockArrayCoo(blocks,blockids,ax)
#end


text:
__precompile__()

module SimpleDiffEq

using Reexport, MuladdMacro
@reexport using DiffEqBase
using StaticArrays
using RecursiveArrayTools

@inline _copy(a::SArray) = a
@inline _copy(a) = copy(a)

include("functionmap.jl")
include("euler_maruyama.jl")
include("rk4/rk4.jl")
include("rk4/gpurk4.jl")
include("rk4/looprk4.jl")
include("tsit5/atsit5_cache.jl")
include("tsit5/tsit5.jl")
include("tsit5/atsit5.jl")
include("tsit5/gpuatsit5.jl")

end # module


text:
import CLFFT

# figure out a gc safe way to store plans.
# weak refs won't work, since the caching should keep them alive.
# But at the end, we need to free all of these, otherwise CLFFT will crash
# at closing time.
# An atexit hook here, which will empty the dictionary seems to introduce racing
# conditions.
#const plan_dict = Dict()
import Base: *, plan_ifft!, plan_fft!, plan_fft, plan_ifft, size, plan_bfft, plan_bfft!

struct CLFFTPlan{Direction, Inplace, T, N} <: Base.FFTW.FFTWPlan{T, Direction, Inplace}
    plan::CLFFT.Plan{T}
    function CLFFTPlan{Direction, Inplace}(A::CLArray{T, N}) where {T, N, Direction, Inplace}
        ctx = context(A)
        p = CLFFT.Plan(T, ctx.context, size(A))
        CLFFT.set_layout!(p, :interleaved, :interleaved)
        if Inplace
            CLFFT.set_result!(p, :inplace)
        else
            CLFFT.set_result!(p, :outofplace)
        end
        CLFFT.set_scaling_factor!(p, Direction, 1f0)
        CLFFT.bake!(p, ctx.queue)
        new{Direction, Inplace, T, N}(p)
    end
end
size(x::CLFFTPlan) = (CLFFT.lengths(x.plan)...,)

# ignore flags, but have them to make it base compatible.
# TODO can we actually implement the flags?
function plan_fft(A::CLArray; flags = nothing, timelimit = Inf)
    CLFFTPlan{:forward, false}(A)
end
function plan_fft!(A::CLArray; flags = nothing, timelimit = Inf)
    CLFFTPlan{:forward, true}(A)
end
function plan_bfft(A::CLArray, region; flags = nothing, timelimit = Inf)
    CLFFTPlan{:backward, false}(A)
end
function plan_bfft!(A::CLArray, region; flags = nothing, timelimit = Inf)
    CLFFTPlan{:backward, true}(A)
end

const _queue_ref = Vector{cl.CmdQueue}(1)
function *(plan::CLFFTPlan{Direction, true, T, N}, A::CLArray{T, N}) where {T, N, Direction}
    _queue_ref[] = context(A).queue
    CLFFT.enqueue_transform(plan.plan, Direction, _queue_ref, buffer(A), nothing)
    A
end
function *(plan::CLFFTPlan{Direction, false, T, N}, A::CLArray{T, N}) where {T, N, Direction}
    _queue_ref[] = context(A).queue
    y = typeof(A)(size(plan))
    CLFFT.enqueue_transform(plan.plan, Direction, _queue_ref, buffer(A), buffer(y))
    y
end