text:
"""
Compare our results with FJ model (and upper bound)

If we directly modify k nodes' innate opinions, each one modify with at most epsilon,
how much will it influence the final measures?

Ref:
    Gaitonde, Jason, Jon Kleinberg, and Eva Tardos.
    "Adversarial perturbations of opinion dynamics in networks." Proceedings of the 21st ACM Conference on Economics and Computation. 2020.

    Chen, Mayee F., and Miklos Z. Racz.
    "Network disruption: maximizing disagreement and polarization in social networks." arXiv preprint arXiv:2003.08377 (2020).
"""

include("../src/Graph.jl")
include("../src/Tools.jl")
include("../src/Algorithm.jl")
include("../src/Sampling.jl")
include("../src/Methods.jl")

using LinearAlgebra
using StatsBase
using JSON
using JLD2

epsilon = 0.1

jsdict = Dict{String, Any}()
tags = []

# Read Graph
buf = split(ARGS[1], ',')
fileName = string("../data/all/", buf[1], ".jld2")

(name, n, m, V, E) = load(fileName, "name", "n", "m", "V", "E")
G = Graph(n, m, V, E)

merge!(jsdict, Dict("name"=>name, "n"=>n, "m"=>m))

# Before Making changes
if ARGS[2] == "1"
    # Uniform distribution
    push!(tags, "Uniform")
    s = load(fileName, "s_Uni")
elseif ARGS[2] == "2"
    # Exponential distribution
    push!(tags, "Exponential")
    s = load(fileName, "s_Exp")
elseif ARGS[2] == "3"
    # Power-law distribution
    push!(tags, "Power-law")
    s = load(fileName, "s_Pow")
end

# Different parameters to select k
if ARGS[3] == "1"
    push!(tags, "0.005")
    k = load(fileName, "k1")
elseif ARGS[3] == "2"
    push!(tags, "0.015")
    k = load(fileName, "k2")
elseif ARGS[3] == "3"
    push!(tags, "0.02")
    k = load(fileName, "k3")
end

merge!(jsdict, Dict("k"=>k))
merge!(jsdict, Dict("s_dist"=>tags[1], "k_par"=>tags[2]))

name = ["z_sum", "aci", "ad", "ap", "aidc"]

##################### Init

# The initial of different measures
_, init_aci, init_ad, init_ap, init_aidc = Approx(G, s)
init_result = [sum(s), init_aci, init_ad, init_ap, init_aidc, init_aidc] # the last one is a upper bound

init_labels = "init_".*name
for i in 1:length(name)
    merge!(jsdict, Dict(init_labels[i]=>init_result[i]))
end

#################### Changes
# Precalculation

L = Symmetric(getL(G))

W = getW(G)
W = (W + W') / 2
#println(issymmetric(W))

Mp = W*W - ones(n, n)/n
Mp = (Mp' + Mp)/2

Md = W*L*W
Md = (Md' + Md)/2

Mic = W*L*L*W
Mic = (Mic' + Mic)/2

Midc = W

s_sqrt = sqrt(sum([i^2 for i in s]))

g = zeros(n)

for (ID, u, v, w) in G.E
    g[v] += 1
    g[u] += 1
end

d_max = maximum(g) / 2

########################################### Upper bound

# Polarization
maxEigen = eigmax(Mp)
uppChange = maxEigen * epsilon^2 * k + min(2*k*epsilon, 2*sqrt(k)*epsilon*maxEigen*s_sqrt)
merge!(jsdict, Dict("FJ_upp_ap"=>uppChange))
#println("polarization ", uppChange)

# Disagreement
maxEigen = eigmax(Md)
uppChange = maxEigen * epsilon^2 * k + min(2*k*epsilon*d_max, 2*sqrt(k)*epsilon*maxEigen*s_sqrt)
merge!(jsdict, Dict("FJ_upp_ad"=>uppChange))
#println("Disagreement ", uppChange)

# Internal conflict
maxEigen = eigmax(Mic)
uppChange = maxEigen * epsilon^2 * k + 2*sqrt(k)*epsilon*maxEigen*s_sqrt
merge!(jsdict, Dict("FJ_upp_aci"=>uppChange))
#println("IC ", uppChange)

# Disagreement-Controversy
maxEigen = eigmax(Midc)
uppChange = maxEigen * epsilon^2 * k + 2*sqrt(k)*epsilon*maxEigen*s_sqrt
merge!(jsdict, Dict("FJ_upp_aidc"=>uppChange))
#println("IDC ", uppChange)

########################################### Algorithm
function greedyFJ(Mx, s, n, k, epsilon)
    s_new = deepcopy(s)

    Lvec = 2*Mx*s
    cSet = Set{Int}()
    AllSet = BitSet(1:n)
    epVec = zeros(n)

    TotalGain = 0

    for i = 1:k
        MaxGain = 0
        MaxIdx = 0
        MaxFlag = 0
        for j in AllSet
            TmpGain1 = Lvec[j]*epsilon + sum([Mx[t, j]*epVec[t] + Mx[j, t]*epVec[t] for t in cSet])*(epsilon) + Mx[j, j]*epsilon^2
            TmpGain2 = -Lvec[j]*epsilon + sum([Mx[t, j]*epVec[t] + Mx[j, t]*epVec[t] for t in cSet])*(-epsilon) + Mx[j, j]*epsilon^2
            if  TmpGain1 > MaxGain
                MaxIdx = j
                MaxGain = TmpGain1
                MaxFlag = 1
            end
            if TmpGain2 > MaxGain
                MaxIdx = j
                MaxGain = TmpGain2
                MaxFlag = 2
            end
        end
        if MaxIdx == 0
            break
        else
            if MaxFlag == 1
                epVec[MaxIdx] = epsilon
                s_new[MaxIdx] = s[MaxIdx] + epsilon
            elseif MaxFlag == 2
                epVec[MaxIdx] = epsilon
                s_new[MaxIdx] = s[MaxIdx] - epsilon
            end
            push!(cSet, MaxIdx)
            delete!(AllSet, MaxIdx)
            TotalGain += MaxGain
        end
    end

    #println("chosen nodes are ", cSet)
    #println("changes is: ", s_new'*Mx*s_new - s'*Mx*s, " or ", TotalGain, " or ", epVec'*Lvec + epVec'*Mx*epVec)
    return TotalGain
end

# Polarization
Change = greedyFJ(Mp, s, n, k, epsilon)
merge!(jsdict, Dict("FJ_ap"=>Change))
#println("polarization ", uppChange)

# Disagreement
Change = greedyFJ(Md, s, n, k, epsilon)
merge!(jsdict, Dict("FJ_ad"=>Change))
#println("Disagreement ", uppChange)

# Internal conflict
Change = greedyFJ(Mic, s, n, k, epsilon)
merge!(jsdict, Dict("FJ_aci"=>Change))
#println("IC ", uppChange)

# Disagreement-Controversy
Change = greedyFJ(Midc, s, n, k, epsilon)
merge!(jsdict, Dict("FJ_aidc"=>Change))
#println("IDC ", uppChange)

# zsum
merge!(jsdict, Dict("FJ_z_sum"=>k*epsilon))

#println("zsum ", k*epsilon)

js = JSON.json(jsdict)

println(jsdict)

open("out/fj_compare.json", "a+") do ff
    write(ff, js)
    write(ff, "\n")
end


text:
# QuEST_jl/src/base/debugging.jl
#

function copyStateFromGPU(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:copyStateFromGPU, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function copyStateToGPU(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:copyStateToGPU, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function getEnvironmentString(env ::QuEST_Types.QuESTEnv, qureg ::QuEST_Types.Qureg) ::String
    cstr = Vector{Cchar}(undef,232)
    ccall(:getEnvironmentString, Cvoid, (QuEST_Types.QuESTEnv, QuEST_Types.Qureg, Ptr{Cchar}), env, qureg, cstr)

    return unsafe_string(pointer(cstr))
end

function initDebugState(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:initDebugState, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

#
#void 	invalidQuESTInputError (const char *errMsg, const char *errFunc)
#An internal function called when invalid arguments are passed to a QuEST API call, which the user can optionally override by redefining.  More...
#

function reportPauliHamil(hamil     ::QuEST_Types.PauliHamil)   ::Nothing
    ccall(:reportPauliHamil, Cvoid, (QuEST_Types.PauliHamil,), hamil)
    return nothing
end

function reportQuESTEnv(env ::QuEST_Types.QuESTEnv) ::Nothing
    ccall(:reportQuESTEnv, Cvoid, (QuEST_Types.QuESTEnv,), env)
    return nothing
end

function reportQuregParams(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:reportQuregParams, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function reportState(qureg ::QuEST_Types.Qureg) ::Nothing
    ccall(:reportState, Cvoid, (QuEST_Types.Qureg,), qureg)
    return nothing
end

function reportStateToScreen(qureg      ::QuEST_Types.Qureg,
                             env        ::QuEST_Types.QuESTEnv,
                             reportRank ::Integer)              ::Nothing
    ccall(:reportStateToScreen,
          Cvoid,
          (QuEST_Types.Qureg, QuEST_Types.QuESTEnv, Cint),
          qureg,              env,                  reportRank)
    nothing
end

function seedQuEST(seedarray ::Base.Vector{Culong}) ::Nothing

    @assert  ! isempty(seedarray)

    ccall(:seedQuEST,
          Cvoid, (Ptr{Culong}, Cint),
          seedarray,
          Cint(length(seedarray)))
    return nothing
end

function seedQuESTDefault() ::Nothing
    ccall(:seedQuESTDefault, Cvoid, () )
    return nothing
end

function syncQuESTEnv(env ::QuEST_Types.QuESTEnv) ::Nothing
    ccall(:syncQuESTEnv, Cvoid,
          (QuEST_Types.QuESTEnv,),
          env)
    return nothing
end

function syncQuESTSuccess(successCode       ::Integer)       ::Cint
    return ccall(:syncQuESTSuccess, Cint,
                 (Cint,),
                 successCode)
end

#EOF

ext: jl


text:
module NeuralNetworks
include("types.jl")
include("Shape.jl")
#include("SpatialConv.jl")

# weighted layers
include("Linear.jl")

# nonlinearities
include("ReLU.jl")
include("Sigmoid.jl")

# containers
include("Serial.jl")

# miscellaneous
include("View.jl")

# criterions
include("MSECriterion.jl")

end


text:
@testset "Parallel.Radiality" begin
    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")
    c = vec(readdlm(joinpath(testdir, "testdata", "graph-50-500-rc.txt"), ','))
    for g in testdigraphs(gint)
        zd = @inferred(Parallel.radiality_centrality(g; parallel=:distributed))
        @test zd == c
        zt = @inferred(Parallel.radiality_centrality(g; parallel=:threads))
        @test zt == c
    end

    g1 = CycleGraph(4)
    add_vertex!(g1)
    add_edge!(g1, 4, 5)
    for g in testgraphs(g1)
        zd = @inferred(Parallel.radiality_centrality(g; parallel=:distributed))
        @test zd ≈ [5 // 6, 3 // 4, 5 // 6, 11 // 12, 2 // 3]
        zt = @inferred(Parallel.radiality_centrality(g; parallel=:threads))
        @test zt ≈ [5 // 6, 3 // 4, 5 // 6, 11 // 12, 2 // 3]
    end
end


text:
struct Kennet <: AbstractSimulator
    wavelet::Vector{Float64}
    dt::Float64
    mopt::Int
    fs::Int
    varnames::Vector{Symbol}
end

function Kennet(;wavelet::Vector{Float64}=[], dt::Float64=1e-3,
        mopt::Int=2, fs::Int=0, varnames::Vector{Symbol}=[:rho, :vp, :dh])
    @assert !isempty(wavelet) "No wavelet provided."
    @assert dt > 0 "Sampling rate dt must be positive."
    @assert mopt >= 0 "Multiple order mopt must be non-negative."

    Kennet(wavelet, dt, mopt, fs, varnames)
end

function (sim::Kennet)(lyr::Dict{Symbol,Vector{Float64}})
    ρ, vₚ, d = (lyr[key] for key in sim.varnames)
    simkennet(ρ, vₚ, d, sim.wavelet, sim.dt, sim.mopt, sim.fs)
end

(sim::Kennet)(p::AbstractPattern) = sim(sample(p))
(sim::Kennet)(ps::Vector{<:AbstractPattern}, cut::UnitRange{Int}=1:length(ps)) = sim(sample.(ps), cut)

function (sim::Kennet)(lyr::Vector{Dict{Symbol,Vector{Float64}}}, cut::UnitRange{Int}=1:length(lyr))
    s = sim(stack(lyr))

    if first(cut) > 0
        t = d2t(lyr)
        len = length(first(values(first(lyr))))
        start = (first(cut) - 1) * len + 1
        stop = last(cut) * len
        # @show t[start], t[stop]
        t_s = sim.dt * (0:(length(s) - 1))
        range = (t_s .>= t[start]) .& (t_s .<= t[stop])
        s = s[range]
    end

    s
end

function (sim::Kennet)(pset::PatternSet; mask::Vector{Int}=[0, 0, 0], cut::UnitRange{Int}=1:length(mask))
    n = length(mask)
    0 <= first(cut) <= last(cut) <= n || error("Attempt to cut $n-element stack at index $cut")

    # sample n patterns from the set
    idx = map(x -> x == 0 ? rand(1:length(pset)) : x, mask)
    patterns = pset.patterns[idx]

    # compute synthetic seismogram
    s = sim(patterns, cut)

    s
end

function (sim::Kennet)(pset::PatternSet, masks::Vector{Vector{Int}}, cut::UnitRange{Int}=1:length(first(masks)))
    nworkers() == 1 && (return [sim(pset, mask, cut) for mask in masks])
    println("starting simulation on $(nworkers()) workers")

    f = mask -> sim(pset, mask, cut)
    return pmap(f, masks)
end

function simkennet(ρ::Vector{T}, vₚ::Vector{T}, d::Vector{T}, wavelet::Vector{T}, dt::T=1e-3, mopt::Int64=2, fs::Int64=0) where {T<:Real}
    nlr = length(ρ)
    n = length(wavelet)

    om = (2π/(n*dt)) .* (0:(n/2 - 1.)) |> collect
    freq = om ./ (2π)

    p₀ = ifft(wavelet)
    w₀ = (1 ./ (ρ[1]*vₚ[1])) .* p₀
    rdhat = zeros(size(om))
    # tdhat = ones(size(om))

    deno = @. ρ[2:end]*vₚ[2:end]+ ρ[1:end-1]*vₚ[1:end-1]
    rd = .- diff(ρ.*vₚ) ./ deno
    td = @. 2*sqrt(ρ[2:end]*vₚ[2:end]*ρ[1:end-1]*vₚ[1:end-1])/deno

    if fs == 1
        pushfirst!(rd, -1)
        pushfirst!(td, 1)
    else
        pushfirst!(rd, 0)
        pushfirst!(td, 1)
    end

    ru = .- rd
    # tu = td

    for j=nlr:-1:1
        ed = exp.((d[j]/vₚ[j])*im*om)

        if mopt == 0
            reverb = ones(size(om))
        elseif mopt == 1
            reverb = 1. .+ ru[j] .*ed .*rdhat .*ed
        else
            reverb = 1. ./ (1. .- ru[j] .*ed .*rdhat .*ed)
        end

        rdhat = rd[j] .+ td[j] .*ed .*rdhat .*ed .*reverb .*td[j]
        #tdhat = tdhat .*ed .*reverb .*td[j]

        if mopt == 1
            dx=findall(real(rdhat) .> 1.); rdhat[dx] = 1. .+ im*imag(rdhat[dx]);
            dx=findall(imag(rdhat) .> 1.); rdhat[dx] = real(rdhat[dx]) .+ 1.0im;
            # dx=findall(real(tdhat) .> 1); tdhat[dx] = 1. + im*imag(tdhat[dx]);
            # dx=findall(imag(tdhat) .> 1); tdhat[dx] = real(tdhat[dx]) + im;
        end
    end

    # tf = hcat(freq, rdhat, tdhat)
    # pz = tdhat .* p₀[1:length(om)]
    wz = rdhat .* p₀[1:length(om)]

    # pz = vcat(pz, 0, conj(reverse(pz[2:end], dims=1)))
    wz = vcat(0, wz[2:end], 0, conj(reverse(wz[2:end], dims=1)))

    # pz = real(fft(pz))
    wz = real(fft(wz))
end

function Base.show(io::IO, sim::Kennet)
    variables = join(sim.varnames, ", ")
    print(io, "Kennet($variables)")
end

function Base.show(io::IO, ::MIME"text/plain", sim::Kennet)
    println(io, "Kennet")
    println(io, "  Variables:         ", sim.varnames)
    println(io, "  Sampling rate Δt:  ", sim.dt, " s")
    println(io, "  Multiples order:   ", sim.mopt)
    println(io, "  Wavelet length:    ", length(sim.wavelet), " points")
end


text:
function strategy_static_weights(returns::Returns, ratio::Vector{Float64})::Weights
    return TimeSeries.map((timestamp, value) -> (timestamp, ratio), returns)
end

function strategy_static_equal_weights(returns::Returns)::Weights
    return TimeSeries.map((timestamp, value) -> (timestamp, repeat([1 / length(colnames(returns))], length(colnames(returns)))), returns)
end


text:
# using NNlib, Test

maxpool_answer_dict = Dict(
    1 => Dict(
        "y"          => [2, 4.],
        "y_nostride" => [2, 3, 4, 5.],
        "y_pad"      => [1, 3, 5.],

        "dx"          => [0, 2, 0, 4, 0.],
        "dx_nostride" => [0, 2, 3, 4, 5.],
        "dx_pad"      => [1, 0, 3, 0, 5.],
    ),
    2 => Dict(
        "y" => [
            7 17.;
            9 19.
        ],
        "y_nostride" => [
            7  12 17;
            8  13 18;
            9  14 19;
            10 15 20.
        ],
        "y_pad" => [
            1  11 16;
            3  13 18;
            5  15 20.
        ],

        "dx" => [
            0 0 0  0;
            0 7 0 17;
            0 0 0  0;
            0 9 0 19;
            0 0 0  0.
        ],
        "dx_nostride" => [
            0  0  0  0;
            0  7 12 17;
            0  8 13 18;
            0  9 14 19;
            0 10 15 20.
        ],
        "dx_pad"      => [
            1 0 11 16;
            0 0  0  0;
            3 0 13 18;
            0 0  0  0;
            5 0 15 20.
        ],
    ),
    3 => Dict(
        "y" => reshape([
            27, 29,
            37, 39.
        ], (2, 2, 1)),
        "y_nostride" => reshape([
            27, 28, 29, 30,
            32, 33, 34, 35,
            37, 38, 39, 40,

            47, 48, 49, 50,
            52, 53, 54, 55,
            57, 58, 59, 60.
        ], (4, 3, 2)),
        "y_pad" => reshape([
             1,  3, 5,
            11, 13, 15,
            16, 18, 20,

            41, 43, 45,
            51, 53, 55,
            56, 58, 60.
        ], (3, 3, 2)),

        "dx" => reshape([
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            0,  0, 0,  0, 0,
            0, 27, 0, 29, 0,
            0,  0, 0,  0, 0,
            0, 37, 0, 39, 0,

            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0.
        ], (5, 4, 3)),
        "dx_nostride" => reshape([
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            0,  0,  0,  0,  0,
            0, 27, 28, 29, 30,
            0, 32, 33, 34, 35,
            0, 37, 38, 39, 40,

            0,  0,  0,  0,  0,
            0, 47, 48, 49, 50,
            0, 52, 53, 54, 55,
            0, 57, 58, 59, 60.
        ], (5, 4, 3)),
        "dx_pad" => reshape([
             1, 0,  3, 0,  5,
             0, 0,  0, 0,  0,
            11, 0, 13, 0, 15,
            16, 0, 18, 0, 20,

            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,
            0, 0, 0, 0, 0,

            41, 0, 43, 0, 45,
             0, 0,  0, 0,  0,
            51, 0, 53, 0, 55,
            56, 0, 58, 0, 60.
        ], (5, 4, 3)),
    )
)

meanpool_answer_dict = Dict(
    1 => Dict(
        "y"          => [1.5, 3.5],
        "y_nostride" => [1.5, 2.5, 3.5, 4.5],
        "y_pad"      => [0.5, 2.5, 4.5],

        "dx"          => [0.75, 0.75, 1.75, 1.75,  0.0],
        "dx_nostride" => [0.75,  2.0,  3.0,  4.0, 2.25],
        "dx_pad"      => [0.25, 1.25, 1.25, 2.25, 2.25],
    ),
    2 => Dict(
        "y" => [
            4.0 14.0;
            6.0 16.0
        ],
        "y_nostride" => [
            4.0  9.0 14.0
            5.0 10.0 15.0
            6.0 11.0 16.0
            7.0 12.0 17.0
        ],
        "y_pad" => [
            0.25  4.25 4.0
            1.25 10.0  8.75
            2.25 12.0  9.75
        ],

        "dx" => [
            1.0 1.0 3.5 3.5;
            1.0 1.0 3.5 3.5;
            1.5 1.5 4.0 4.0;
            1.5 1.5 4.0 4.0;
            0.0 0.0 0.0 0.0
        ],
        "dx_nostride" => [
            1.0  3.25  5.75 3.5;
            2.25 7.0  12.0  7.25;
            2.75 8.0  13.0  7.75;
            3.25 9.0  14.0  8.25;
            1.75 4.75  7.25 4.25
        ],
        "dx_pad"      => [
            0.0625 1.0625 1.0625 1.0;
            0.3125 2.5    2.5    2.1875;
            0.3125 2.5    2.5    2.1875;
            0.5625 3.0    3.0    2.4375;
            0.5625 3.0    3.0    2.4375
        ],
    ),
    3 => Dict(
        "y" => reshape([
            14.0, 16.0,
            24.0, 26.0
        ], (2, 2, 1)),
        "y_nostride" => reshape([
            14.0, 15.0, 16.0, 17.0,
            19.0, 20.0, 21.0, 22.0,
            24.0, 25.0, 26.0, 27.0,

            34.0, 35.0, 36.0, 37.0,
            39.0, 40.0, 41.0, 42.0,
            44.0, 45.0, 46.0, 47.0
        ], (4, 3, 2)),
        "y_pad" => reshape([
            0.125, 0.625, 1.125,
            2.125, 5.0,   6.0,
            2.0,   4.375, 4.875,

             7.75, 16.25, 17.25,
            19.25, 40.0,  42.0,
            11.5,  23.75, 24.75,
        ], (3, 3, 2)),

        "dx" => reshape([
            1.75, 1.75, 2.0, 2.0, 0.0,
            1.75, 1.75, 2.0, 2.0, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,

            1.75, 1.75, 2.0, 2.0, 0.0,
            1.75, 1.75, 2.0, 2.0, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,
            3.0, 3.0, 3.25, 3.25, 0.0,

            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0,
        ], (5, 4, 3)),
        "dx_nostride" => reshape([
            1.75,   3.625,  3.875,  4.125, 2.125,
            4.125,  8.5,    9.0,    9.5,   4.875,
            5.375, 11.0,   11.5,   12.0,   6.125,
            3.0,    6.125,  6.375,  6.625, 3.375,

             6.0,  12.25, 12.75, 13.25,  6.75,
            13.25, 27.0,  28.0,  29.0,  14.75,
            15.75, 32.0,  33.0,  34.0,  17.25,
             8.5,  17.25, 17.75, 18.25,  9.25,

             4.25,   8.625,  8.875,  9.125,  4.625,
             9.125, 18.5,   19.0,   19.5,    9.875,
            10.375, 21.0,   21.5,   22.0,   11.125,
             5.5,   11.125, 11.375, 11.625,  5.875
        ], (5, 4, 3)),
        "dx_pad" => reshape([
            0.015625, 0.078125, 0.078125, 0.140625, 0.140625,
            0.265625, 0.625,    0.625,    0.75,     0.75,
            0.265625, 0.625,    0.625,    0.75,     0.75,
            0.25,     0.546875, 0.546875, 0.609375, 0.609375,

            0.96875, 2.03125, 2.03125, 2.15625, 2.15625,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            1.4375,  2.96875, 2.96875, 3.09375, 3.09375,

            0.96875, 2.03125, 2.03125, 2.15625, 2.15625,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            2.40625, 5.0,     5.0,     5.25,    5.25,
            1.4375,  2.96875, 2.96875, 3.09375, 3.09375
        ], (5, 4, 3)),
    )
)

for rank in (1, 2, 3)
    @testset "pool$(rank)d" begin
        for (pool, ∇pool, answer_dict) in (
                # Main API name
                (maxpool, ∇maxpool, maxpool_answer_dict),
                (meanpool, ∇meanpool, meanpool_answer_dict),

                # _direct name
                (NNlib.maxpool_direct, NNlib.∇maxpool_direct, maxpool_answer_dict),
                (NNlib.meanpool_direct, NNlib.∇meanpool_direct, meanpool_answer_dict),)

            @testset "$(pool)$(rank)d" begin
                y = answer_dict[rank]["y"]
                y_nostride = answer_dict[rank]["y_nostride"]
                y_pad = answer_dict[rank]["y_pad"]
                dx = answer_dict[rank]["dx"]
                dx_nostride = answer_dict[rank]["dx_nostride"]
                dx_pad = answer_dict[rank]["dx_pad"]

                x = reshape(Float64[1:prod(size(dx));], size(dx)..., 1, 1)

                # A "drop channels and batch dimension" helper
                ddims(x) = dropdims(x, dims=(rank + 1, rank + 2))

                # Let's ensure that a 1x1x1 pooling kernel always just returns `x`
                @test pool(x, PoolDims(x, 1)) == x

                # Test vanilla pooling
                pdims = PoolDims(x, 2)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx

                # Strided pooling
                pdims = PoolDims(x, 2; stride=1)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y_nostride
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx_nostride

                # Padded pooling
                pdims = PoolDims(x, 2; padding=1)
                y_hat = pool(x, pdims)
                @test ddims(y_hat) == y_pad
                @test ddims(∇pool(y_hat, y_hat, x, pdims)) == dx_pad
            end
        end
    end
end

@testset "Pooling - Check Sizes" begin
    x = rand(10, 10, 3, 10)
    @test size(maxpool(x, (2, 2))) == (5, 5, 3, 10)
    @test size(maxpool(x, (2, 2); pad=(1, 1), stride=(2, 2))) == (6, 6, 3, 10)
    @test size(meanpool(x, (2, 2))) == (5, 5, 3, 10)
    @test size(meanpool(x, (2, 2); pad=(1, 1), stride=(2, 2))) == (6, 6, 3, 10)
end

# Add another test for 2d maxpool that uses an odd-length size:
@testset "Issue #133" begin
    x = reshape([(1.:9.)...], 3, 3, 1, 1)
    pdims = PoolDims(size(x), (2, 2), padding=(1, 1), stride=(2, 2))
    y = maxpool(x, pdims)

    dy = y .* 0 .+ 1
    dx = ∇maxpool(dy, y, x, pdims)
    @test dx[:,:,1,1] == [1.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 1.0]
end

# test "true" strided case, see https://github.com/FluxML/NNlib.jl/issues/205

# obtained with
# using FiniteDifferences
maxpool_answer_nature = Dict(
    "rank1" => Dict(
        # kernel size 2, stride 1, pad 0
        "k2s1p0" => (size = (2,),
            stride = 1,
            pad = 0,

            x = reshape([
                0.0299635,  0.233456,  0.596161,   0.161514,  0.0094027
            ], 5, 1, 1), # width, channel, batch_size

            dx_maxpool = reshape([
                 0.0, 1.0, 2.0, 1.0, 0.0
            ], 5, 1, 1),

            dx_meanpool = reshape([
                 0.5, 1.0, 1.0, 1.0, 0.5
            ], 5, 1, 1),),
        "k2s1p1" => (size = (2,),
            stride = 1,
            pad = 1,

            x = reshape([
                0.0299635,  0.233456,  0.596161,   0.161514,  0.0094027
            ], 5, 1, 1),

            dx_maxpool = reshape([
                 1.0, 1.0, 2.0, 1.0, 1.0
            ], 5, 1, 1),

            dx_meanpool = reshape([
                 1.0, 1.0, 1.0, 1.0, 1.0
            ], 5, 1, 1),),
        "k3s1p1" => (size = (3,),
            stride = 1,
            pad = 1,

            x = reshape([
                0.0299635,  0.233456,  0.596161,   0.161514,  0.0094027
            ], 5, 1, 1),

            dx_maxpool = reshape([
                 0.0, 1.0, 3.0, 1.0, 0.0
            ], 5, 1, 1),

            dx_meanpool = reshape([
                 0.6666666666, 1.0, 1.0, 1.0, 0.6666666666
            ], 5, 1, 1),),
        "k3s2p1" => (size = (3,),
            stride = 2,
            pad = 1,

            x = reshape([
                0.0299635,  0.233456,  0.596161,   0.161514,  0.0094027
            ], 5, 1, 1),

            dx_maxpool = reshape([
                 0.0, 1.0, 1.0, 1.0, 0.0
            ], 5, 1, 1),

            dx_meanpool = reshape([
                 0.333333333,
                 0.666666666,
                 0.333333333,
                 0.666666666,
                 0.333333333,
            ], 5, 1, 1),)
    ),
    "rank2" => Dict(
        # kernel size 2, stride 1, pad 0
        "k2s1p0" => (size = (2, 2),
            stride = 1,
            pad = 0,

            x = reshape([
                0.0299635  0.233456  0.596161   0.161514  0.0094027
                0.389984   0.235158  0.579525   0.301893  0.561358
                0.0830242  0.483759  0.914904   0.253871  0.820061
                0.425287   0.53451   0.0405225  0.729861  0.403925
                0.473724   0.571418  0.558427   0.552183  0.561624
            ], 5, 5, 1, 1),

            dx_maxpool = reshape([
                0.0  0.0  2.0  0.0  0.0
                1.0  0.0  0.0  0.0  1.0
                0.0  1.0  4.0  0.0  2.0
                0.0  1.0  0.0  2.0  0.0
                0.0  2.0  0.0  0.0  0.0
            ], 5, 5, 1, 1),

            dx_meanpool = reshape([
                0.25  0.5  0.5  0.5  0.25
                0.5   1.0  1.0  1.0  0.5
                0.5   1.0  1.0  1.0  0.5
                0.5   1.0  1.0  1.0  0.5
                0.25  0.5  0.5  0.5  0.25
            ], 5, 5, 1, 1)),
        "k2s1p1" => (size = (2, 2),
            stride = 1,
            pad = 1,

            x = reshape([
                0.0299635  0.233456  0.596161   0.161514  0.0094027
                0.389984   0.235158  0.579525   0.301893  0.561358
                0.0830242  0.483759  0.914904   0.253871  0.820061
                0.425287   0.53451   0.0405225  0.729861  0.403925
                0.473724   0.571418  0.558427   0.552183  0.561624
            ], 5, 5, 1, 1),

            dx_maxpool = reshape([
                1.0  1.0  4.0  1.0  1.0
                3.0  0.0  0.0  0.0  2.0
                0.0  1.0  4.0  0.0  4.0
                1.0  1.0  0.0  2.0  0.0
                2.0  4.0  1.0  0.0  3.0
            ], 5, 5, 1, 1),

            dx_meanpool = reshape([
                1.0  1.0  1.0  1.0  1.0
                1.0  1.0  1.0  1.0  1.0
                1.0  1.0  1.0  1.0  1.0
                1.0  1.0  1.0  1.0  1.0
                1.0  1.0  1.0  1.0  1.0
            ], 5, 5, 1, 1)),
        "k3s1p1" => (size = (3, 3),
            stride = 1,
            pad = 1,

            x = reshape([
                0.0299635  0.233456  0.596161   0.161514  0.0094027
                0.389984   0.235158  0.579525   0.301893  0.561358
                0.0830242  0.483759  0.914904   0.253871  0.820061
                0.425287   0.53451   0.0405225  0.729861  0.403925
                0.473724   0.571418  0.558427   0.552183  0.561624
            ], 5, 5, 1, 1),

            dx_maxpool = reshape([
                0.0  0.0  3.0  0.0  0.0
                1.0  0.0  0.0  0.0  1.0
                0.0  1.0  9.0  0.0  3.0
                0.0  1.0  0.0  3.0  0.0
                0.0  3.0  0.0  0.0  0.0
            ], 5, 5, 1, 1),

            dx_meanpool = reshape([
                0.444444  0.666667  0.666667  0.666667  0.444444
                0.666667  1.0       1.0       1.0       0.666667
                0.666667  1.0       1.0       1.0       0.666667
                0.666667  1.0       1.0       1.0       0.666667
                0.444444  0.666667  0.666667  0.666667  0.444444
            ], 5, 5, 1, 1)),
        "k3s2p1" => (size = (3, 3),
            stride = 2,
            pad = 1,

            x = reshape([
                0.0299635  0.233456  0.596161   0.161514  0.0094027
                0.389984   0.235158  0.579525   0.301893  0.561358
                0.0830242  0.483759  0.914904   0.253871  0.820061
                0.425287   0.53451   0.0405225  0.729861  0.403925
                0.473724   0.571418  0.558427   0.552183  0.561624
            ], 5, 5, 1, 1),

            dx_maxpool = reshape([
                0.0  0.0  1.0  0.0  0.0
                1.0  0.0  0.0  0.0  1.0
                0.0  0.0  1.0  0.0  1.0
                0.0  1.0  0.0  2.0  0.0
                0.0  1.0  0.0  0.0  0.0
            ], 5, 5, 1, 1),

            dx_meanpool = reshape([
                0.111111  0.222222  0.111111  0.222222  0.111111
                0.222222  0.444444  0.222222  0.444444  0.222222
                0.111111  0.222222  0.111111  0.222222  0.111111
                0.222222  0.444444  0.222222  0.444444  0.222222
                0.111111  0.222222  0.111111  0.222222  0.111111
            ], 5, 5, 1, 1))
    ),
    "rank3" => Dict(
        # kernel size 2, stride 1, pad 0
        "k2s1p0" => (size = (2, 2, 2),
            stride = 1,
            pad = 0,

            x = reshape(cat([
                    0.82584   0.416818  0.92668   0.471931
                    0.798798  0.131608  0.344556  0.79681
                    0.716898  0.320672  0.24453   0.288568
                    0.261484  0.258469  0.121916  0.0685961
                ],
                [
                    0.73934   0.16631    0.525109   0.0223458
                    0.164918  0.790875   0.444085   0.469671
                    0.116848  0.359845   0.0653075  0.804886
                    0.525431  0.0402844  0.846814   0.84876
                ],
                [
                    0.709245  0.325828  0.715952  0.719116
                    0.576722  0.405659  0.770104  0.259131
                    0.640221  0.28811   0.129229  0.97571
                    0.953795  0.1316    0.94538   0.705337
                ],dims=3), 4,4,3,1,1),

            dx_maxpool = reshape(cat([
                     1.0  0.0  2.0  0.0
                     1.0  0.0  0.0  0.0
                     1.0  0.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                ],
                [
                     0.0  0.0  0.0  0.0
                     0.0  5.0  0.0  0.0
                     0.0  0.0  0.0  1.0
                     0.0  0.0  1.0  1.0
                ],
                [
                     0.0  0.0  0.0  0.0
                     0.0  0.0  1.0  0.0
                     0.0  0.0  0.0  2.0
                     1.0  0.0  1.0  0.0
                ],dims=3), 4,4,3,1,1),

            dx_meanpool = reshape(cat([
                     0.125  0.25  0.25  0.125
                     0.25   0.5   0.5   0.25
                     0.25   0.5   0.5   0.25
                     0.125  0.25  0.25  0.125
                ],
                [
                     0.25  0.5  0.5  0.25
                     0.5   1.0  1.0  0.5
                     0.5   1.0  1.0  0.5
                     0.25  0.5  0.5  0.25
                ],
                [
                     0.125  0.25  0.25  0.125
                     0.25   0.5   0.5   0.25
                     0.25   0.5   0.5   0.25
                     0.125  0.25  0.25  0.125
                ],dims=3), 4,4,3,1,1)),
        "k2s1p1" => (size = (2, 2, 2),
            stride = 1,
            pad = 1,

            x = reshape(cat([
                    0.82584   0.416818  0.92668   0.471931
                    0.798798  0.131608  0.344556  0.79681
                    0.716898  0.320672  0.24453   0.288568
                    0.261484  0.258469  0.121916  0.0685961
                ],
                [
                    0.73934   0.16631    0.525109   0.0223458
                    0.164918  0.790875   0.444085   0.469671
                    0.116848  0.359845   0.0653075  0.804886
                    0.525431  0.0402844  0.846814   0.84876
                ],
                [
                    0.709245  0.325828  0.715952  0.719116
                    0.576722  0.405659  0.770104  0.259131
                    0.640221  0.28811   0.129229  0.97571
                    0.953795  0.1316    0.94538   0.705337
                ],dims=3), 4,4,3,1,1),

            dx_maxpool = reshape(cat([
                     8.0  0.0  8.0  2.0
                     4.0  0.0  1.0  4.0
                     4.0  1.0  0.0  2.0
                     2.0  1.0  1.0  1.0
                ],
                [
                     3.0  0.0  0.0  0.0
                     0.0  5.0  0.0  0.0
                     0.0  0.0  0.0  2.0
                     2.0  0.0  2.0  5.0
                ],
                [
                     4.0  0.0  2.0  6.0
                     0.0  0.0  4.0  0.0
                     3.0  0.0  0.0  8.0
                     8.0  0.0  6.0  1.0
                ],dims=3), 4,4,3,1,1),

            dx_meanpool = reshape(cat([
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                ],
                [
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                ],
                [
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                     1.0  1.0  1.0  1.0
                ],dims=3), 4,4,3,1,1)),
        "k3s1p1" => (size = (3, 3, 2),
            stride = 1,
            pad = 1,

            x = reshape(cat([
                    0.82584   0.416818  0.92668   0.471931
                    0.798798  0.131608  0.344556  0.79681
                    0.716898  0.320672  0.24453   0.288568
                    0.261484  0.258469  0.121916  0.0685961
                ],
                [
                    0.73934   0.16631    0.525109   0.0223458
                    0.164918  0.790875   0.444085   0.469671
                    0.116848  0.359845   0.0653075  0.804886
                    0.525431  0.0402844  0.846814   0.84876
                ],
                [
                    0.709245  0.325828  0.715952  0.719116
                    0.576722  0.405659  0.770104  0.259131
                    0.640221  0.28811   0.129229  0.97571
                    0.953795  0.1316    0.94538   0.705337
                ],dims=3), 4,4,3,1,1),

            dx_maxpool = reshape(cat([
                     4.0  0.0  12.0  0.0
                     3.0  0.0   0.0  2.0
                     3.0  1.0   0.0  1.0
                     0.0  0.0   0.0  0.0
                ],
                [
                     0.0  0.0  0.0  0.0
                     0.0  5.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                     0.0  0.0  2.0  4.0
                ],
                [
                     2.0  0.0  0.0   0.0
                     0.0  0.0  5.0   0.0
                     0.0  0.0  0.0  12.0
                     8.0  0.0  0.0   0.0
                ],dims=3), 4,4,3,1,1),

            dx_meanpool = reshape(cat([
                     0.444444  0.666667  0.666667  0.444444
                     0.666667  1.0       1.0       0.666667
                     0.666667  1.0       1.0       0.666667
                     0.444444  0.666667  0.666667  0.444444
                ],
                [
                     0.444444  0.666667  0.666667  0.444444
                     0.666667  1.0       1.0       0.666667
                     0.666667  1.0       1.0       0.666667
                     0.444444  0.666667  0.666667  0.444444
                ],
                [
                     0.444444  0.666667  0.666667  0.444444
                     0.666667  1.0       1.0       0.666667
                     0.666667  1.0       1.0       0.666667
                     0.444444  0.666667  0.666667  0.444444
                ],dims=3), 4,4,3,1,1)),
        "k3s2p1" => (size = (3, 3, 2),
            stride = 2,
            pad = 1,

            x = reshape(cat([
                    0.82584   0.416818  0.92668   0.471931
                    0.798798  0.131608  0.344556  0.79681
                    0.716898  0.320672  0.24453   0.288568
                    0.261484  0.258469  0.121916  0.0685961
                ],
                [
                    0.73934   0.16631    0.525109   0.0223458
                    0.164918  0.790875   0.444085   0.469671
                    0.116848  0.359845   0.0653075  0.804886
                    0.525431  0.0402844  0.846814   0.84876
                ],
                [
                    0.709245  0.325828  0.715952  0.719116
                    0.576722  0.405659  0.770104  0.259131
                    0.640221  0.28811   0.129229  0.97571
                    0.953795  0.1316    0.94538   0.705337
                ],dims=3), 4,4,3,1,1),

            dx_maxpool = reshape(cat([
                     1.0  0.0  1.0  0.0
                     1.0  0.0  0.0  1.0
                     0.0  0.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                ],
                [
                     0.0  0.0  0.0  0.0
                     0.0  2.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                ],
                [
                     0.0  0.0  0.0  0.0
                     0.0  0.0  0.0  0.0
                     0.0  0.0  0.0  1.0
                     1.0  0.0  0.0  0.0
                ],dims=3), 4,4,3,1,1),

            dx_meanpool = reshape(cat([
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.111111   0.222222  0.111111   0.111111
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.0555556  0.111111  0.0555556  0.0555556
                ],
                [
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.111111   0.222222  0.111111   0.111111
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.0555556  0.111111  0.0555556  0.0555556
                ],
                [
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.111111   0.222222  0.111111   0.111111
                     0.0555556  0.111111  0.0555556  0.0555556
                     0.0555556  0.111111  0.0555556  0.0555556
                ],dims=3), 4,4,3,1,1))
    )
)

@testset "more maxpool and meanpool tests" begin
    # issue #205
    function check(config, T)
        # CHECK DEFAULT
        pdims = PoolDims(config.x, config.size; stride=config.stride, padding=config.pad)
        x = T.(config.x)
        y_maxpool = NNlib.maxpool(x, pdims)
        y_meanpool = NNlib.meanpool(x, pdims)
        dy = ones(T, size(y_maxpool)...) # size(y_maxpool) == size(y_meanpool)
        @test isapprox(config.dx_maxpool, NNlib.∇maxpool(dy, y_maxpool, x, pdims), rtol=1e-5)
        @test isapprox(config.dx_meanpool, NNlib.∇meanpool(dy, y_meanpool, x, pdims), rtol=1e-5)
        # CHECK DIRECT
        y_maxpool_dir = NNlib.maxpool_direct(x, pdims)
        y_meanpool_dir = NNlib.meanpool_direct(x, pdims)
        @test y_maxpool_dir ≈ y_maxpool  atol = 1e-6
        @test isapprox(config.dx_maxpool, NNlib.∇maxpool_direct(dy, y_maxpool_dir, x, pdims), rtol=1e-5)
        @test isapprox(config.dx_meanpool, NNlib.∇meanpool_direct(dy, y_meanpool_dir, x, pdims), rtol=1e-5)

        # CHECK NNPACK
        if NNlib.is_nnpack_available() && T == Float32
            if NNlib.nnpack_supported_operation(pdims)
                y_maxpool_nnp = NNlib.maxpool_nnpack(x, pdims)
                @test y_maxpool_nnp ≈ y_maxpool  atol = 1e-6
                # NNPACK maxpool gradient still missing
                # @test isapprox(config.dx_maxpool, NNlib.∇maxpool_nnpack(dy, y_maxpool_nnp, config.x, pdims), rtol=1e-5)
            end
        end
    end

    for (rank_name, config_dict) in maxpool_answer_nature
        for (setting_name, config) in config_dict
            for T in (Float32, Float64)
                check(config, T)
            end
        end
    end

    # issue 210
    x, k = rand(Float32, 5, 2, 1, 3), (2, 1)
    pdims1 = NNlib.PoolDims(x, k, padding=1, stride=1)
    pdims2 = NNlib.PoolDims(x, k, padding=(1, 0, 0, 0), stride=1)
    @test maxpool(x, pdims1) isa Array{Float32,4}
    @test maxpool(x, pdims2) isa Array{Float32,4}

    # issue #229
    x = ones(Float32, 4, 4, 1, 1) .* -1
    pool = meanpool(x, PoolDims(x, 2, padding=1))
    valid = reshape([
    -0.25,  -0.5,  -0.25,
    -0.5,   -1.0,  -0.5,
    -0.25,  -0.5,  -0.25], (3, 3, 1, 1))
    @test all(pool .== valid)

    # if NNlib.is_nnpack_available()
    #     if NNlib.nnpack_supported_operation(pdims1)
    #         @test NNlib.maxpool_nnpack(x, pdims1) isa Array{Float32, 4}
    #     end
    #     if NNlib.nnpack_supported_operation(pdims2)
    #         print("you should not see this")
    #         @test NNlib.maxpool_nnpack(x, pdims2) isa Array{Float32, 4}
    #     end
    # end
end

@testset "AutoDiff: spatial_rank=$spatial_rank" for spatial_rank in (1, 2)
  x = rand(rng, repeat([10], spatial_rank)..., 3, 2)
  pdims = PoolDims(x, 2)
  gradtest(x -> maxpool(x, pdims), x; broken=spatial_rank <= 2)
  gradtest(x -> meanpool(x, pdims), x)
  gradtest(x -> sum(maxpool(x, pdims)), x, broken=spatial_rank == 2)
  gradtest(x -> sum(meanpool(x, pdims)), x)

  #https://github.com/FluxML/NNlib.jl/issues/188
  k = ntuple(_ -> 2, spatial_rank)  # Kernel size of pool in ntuple format
  gradtest(x -> maxpool(x, k), x; broken=spatial_rank <= 2)
  gradtest(x -> meanpool(x, k), x)
  gradtest(x -> sum(maxpool(x, k)), x, broken=spatial_rank == 2)
  gradtest(x -> sum(meanpool(x, k)), x)
end


text:
struct Null end
CachedArrays.isnull(::Null) = true

mutable struct TestBlock
    isnull::Bool
    next::Union{TestBlock,Null}
    previous::Union{TestBlock,Null}
    address::Int
end

TestBlock(isnull = true, address::Integer = 0) = TestBlock(isnull, Null(), Null(), address)

CachedArrays.isnull(x::TestBlock) = x.isnull
CachedArrays.address(x::TestBlock) = x.address
Base.isless(a::TestBlock, b::TestBlock) = CachedArrays.address(a) < CachedArrays.address(b)

@testset "Testing Freelist" begin
    isnull(x) = CachedArrays.isnull(x)
    null = Null()

    @test isnull(TestBlock())

    @testset "Swapping" begin
        #####
        ##### Test Swapping - four cases:
        #####

        # Null - A - B - Null
        a = TestBlock(false, null, null, 0)
        b = TestBlock(false, null, null, 1)

        a.next = b
        b.previous = a
        CachedArrays.swap!(a, b)
        @test !isnull(a)
        @test !isnull(b)
        @test isnull(b.previous)
        @test b.next === a
        @test a.previous === b
        @test isnull(a.next)

        # C - A - B - Null
        a = TestBlock(false, null, null, 0)
        b = TestBlock(false, null, null, 1)
        c = TestBlock(false, null, null, 2)

        c.next = a

        a.next = b
        a.previous = c

        b.previous = a
        CachedArrays.swap!(a, b)

        # tests
        @test isnull(c.previous)
        @test c.next === b

        @test b.previous === c
        @test b.next === a

        @test a.previous === b
        @test isnull(a.next)

        # Null - A - B - D
        a = TestBlock(false, null, null, 0)
        b = TestBlock(false, null, null, 1)
        d = TestBlock(false, null, null, 2)

        a.next = b

        b.previous = a
        b.next = d

        d.previous = b
        CachedArrays.swap!(a, b)

        @test isnull(b.previous)
        @test b.next === a

        @test a.previous === b
        @test a.next === d
        @test d.previous === a
        @test isnull(d.next)

        # C - A - B - D
        c = TestBlock(false, null, null, 3)
        a = TestBlock(false, null, null, 0)
        b = TestBlock(false, null, null, 1)
        d = TestBlock(false, null, null, 2)

        c.next = a
        a.previous = c
        a.next = b
        b.previous = a
        b.next = d
        d.previous = b

        CachedArrays.swap!(a, b)
        @test isnull(c.previous)
        @test c.next === b
        @test b.previous === c
        @test b.next === a
        @test a.previous === b
        @test a.next === d
        @test d.previous === a
        @test isnull(d.next)
    end

    @testset "Freelist" begin
        list = CachedArrays.Freelist{TestBlock}()
        @test length(list) == 0
        @test isempty(list)
        @test issorted(list)

        # Push blocks in reverse order so we can test out sorting.
        block0 = TestBlock(false, 0)
        push!(list, block0)

        @test !isempty(list)
        @test length(list) == 1
        @test issorted(list)
        @test CachedArrays.getbase(list) === block0

        # Pop an empty list
        popped = pop!(list)
        @test popped === block0
        @test isempty(list)

        push!(list, block0)

        # Add another block, with a higher address
        block1 = TestBlock(false, 1)
        block1.previous = block0

        # Make sure `previous` fields get cleared
        @test !isnull(block1.previous)
        push!(list, block1)
        @test isnull(block1.previous)

        @test length(list) == 2
        @test !issorted(list)

        # Do some ordering checks
        blocks = collect(list)
        @test blocks[1] === block1
        @test blocks[2] === block0

        @test isnull(block1.previous)
        @test block1.next === block0
        @test block0.previous === block1
        @test isnull(block0.next)

        # Add yet another block
        block2 = TestBlock(false, 2)
        push!(list, block2)

        @test length(list) == 3
        blocks = collect(list)
        @test blocks[1] === block2
        @test blocks[2] === block1
        @test blocks[3] === block0

        # Try to pop a block off
        block = pop!(list)
        @test block === block2
        @test length(list) == 2

        # Make sure the block1 got nulled out.
        @test isnull(block1.previous)
        push!(list, block2)

        # Now try sorting
        @test length(list) == 3
        @test !issorted(list)
        sort!(list)
        @test length(list) == 3
        @test issorted(list)
        blocks = collect(list)
        @test blocks[1] === block0
        @test blocks[2] === block1
        @test blocks[3] === block2
    end

    @testset "FreelistPtr" begin
        list_ref = Ref(CachedArrays.Block())
        list = CachedArrays.FreelistPtr(
            Base.unsafe_convert(Ptr{CachedArrays.Block}, list_ref)
        )

        # Allocate a collection of blocks in an array to provide stable addresses
        nblocks = 3
        raw_memory = zeros(UInt8, nblocks * CachedArrays.headersize())
        raw_blocks = map(1:nblocks) do i
            return CachedArrays.Block(
                Ptr{Nothing}(pointer(raw_memory, 1 + i * CachedArrays.headersize()))
            )
        end

        block0 = raw_blocks[1]
        block1 = raw_blocks[2]
        block2 = raw_blocks[3]

        @test length(list) == 0
        @test isempty(list)
        @test issorted(list)

        # Push blocks in reverse order so we can test out sorting.
        push!(list, block0)

        @test !isempty(list)
        @test length(list) == 1
        @test issorted(list)
        @test CachedArrays.getbase(list) === block0

        # Pop an empty list
        popped = pop!(list)
        @test popped === block0
        @test isempty(list)

        push!(list, block0)

        # Add another block, with a higher address
        block1.previous = block0

        # Make sure `previous` fields get cleared
        @test !isnull(block1.previous)
        push!(list, block1)
        @test isnull(block1.previous)

        @test length(list) == 2
        @test !issorted(list)

        # Do some ordering checks
        blocks = collect(list)
        @test blocks[1] === block1
        @test blocks[2] === block0

        @test isnull(block1.previous)
        @test block1.next === block0
        @test block0.previous === block1
        @test isnull(block0.next)

        # Add yet another block
        push!(list, block2)

        @test length(list) == 3
        blocks = collect(list)
        @test blocks[1] === block2
        @test blocks[2] === block1
        @test blocks[3] === block0

        # Try to pop a block off
        block = pop!(list)
        @test block === block2
        @test length(list) == 2

        # Make sure the block1 got nulled out.
        @test isnull(block1.previous)
        push!(list, block2)

        # Now try sorting
        @test length(list) == 3
        @test !issorted(list)
        sort!(list)
        @test length(list) == 3
        @test issorted(list)
        blocks = collect(list)
        @test blocks[1] === block0
        @test blocks[2] === block1
        @test blocks[3] === block2
    end
end


text:
using HttpCommon

function readpdf2(path::String)
    jar = "pdfextract-0.2.3.jar"
    pdfstr = readstring(`java -classpath $jar PDFExtractor $path`)
    tokens = String[]
    lines = split(pdfstr, "\n")
    push!(lines, "")
    for line in lines
        if isempty(line)
            push!(tokens, " ")
        else
            items = Vector{String}(split(line,'\t'))
            id = parse(Int, items[1])
            page = parse(Int, items[2])
            items[3] == "TEXT" || continue
            c = items[4]
            length(c) == 1 || (c = "*")
            push!(tokens, c)
        end
    end

    # create sentences
    sents = String[]
    words = String[]
    for t in tokens
        t == " " && isempty(words) && continue
        push!(words, t)
        if t == "."
            isempty(words) && continue
            push!(sents, join(words))
            empty!(words)
        end
    end
    sents
end

function write_xhtml(filename::String, strs::Vector{String})
    strs = map(escapeHTML, strs)
    body = join(map(s -> "
$s
\n", strs))
    xhtml = """
    <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">

    $body

    """
    open(filename, "w") do f
        println(f, xhtml)
    end
end

function main(pdffile::String, dictfile::String)
    sents = readpdf(pdffile)
    dict = open(readlines, dictfile)
    # dict = map(lowercase, dict)

    offset = 0
    count = 1
    for sent in sents
        # sent = lowercase(sent)
        for word in dict
            isempty(word) && continue
            i = searchindex(sent, word)
            i == 0 && continue
            i = offset + 1
            j = i + length(word)
            anno = """
            [$count]
            type = "span"
            position = [$i, $j]
            label = "auto:solvent"
            text = "$word"
            """
            count += 1
        end
        offset += length(sent) + 1
    end

    offset = 0
    word = "cyclohexane"
    count = 1
    for sent in sents
        i = offset + 1
        j = i + 1
        data = """
        [$count]
        type = "span"
        position = [$i, $j]
        label = "auto:solvent"
        text = "cyclohexane"
        """
        println(f, data)
        #println(f, "")
        #i = searchindex(sent, word)
        #if i > 0
            #i += offset
            #j = i + length(word)
            #println("$i, $j")
        #end
        offset += length(sent) + 1
        count += 1
    end
end

filename = "Macromolecules_13_3_653-656_1980.pdf"
main(filename)


text:
abstract type AbstractState end
abstract type AbstractMutableState <: AbstractState end
abstract type AbstractImmutableState <: AbstractState end


text:
using Distributions, ColorSchemes
using Plots
gr()

export fire

"""
cannonball_trajectory(v::Real,θ::Real,t)

this function return the cannonball [x,y] coords at time t for
for initial data v,θ
"""
function cannonball_trajectory(θ::Real, v::Real, width::Real, t)
    g = 9.81
    x,y = 0.5*width .+ v*cos(θ).*t, v*sin(θ).*t - 0.5*g*t.^2
    return (x,y)
end # function cannonball_trajectory

"""
    random_tries(gen_size::Real; max_speed::Real = 20)

this function create random population of size = 'size' for the cannonball
for the velocity 'v' and angle 'θ'
v ∈ [0, typemax(Int64)]
θ ∈ [0, π]
"""
function random_tries(gen_size::Real; max_speed::Real = 20)::Matrix
    generation = Matrix{Real}(undef,gen_size,2)
    for ii in 1:gen_size
        generation[ii,1] = π * rand()
        generation[ii,2] = max_speed * rand()
    end
    return generation
end # function random_tries

"""
hit_coordinate(θ::Real,v::Real,width::Real)

this function find the coordinates of the cannonballs at the bag edge
"""
function hit_coordinate(θ::Real,v::Real,width::Real)
    x = 0.5 * width
    x_hit = width
    if θ > π/2
        x = -x
        x_hit = 0
    end
    t = x/(v*cos(θ))
    y = v*t*sin(θ) - 0.5 * 9.81 * t^2
    y = y < 0 ? 0 : y
    return x_hit, y
 end # function hit_coordinate

"""
    escaped(θ::Real, v::Real, width::Real, height::Real)

this function check if the cannonball escaped the bag
"""
function escaped(θ::Real, v::Real, width::Real, height::Real)
    x_hit, y_hit = hit_coordinate(θ, v, width)
    return (x_hit == 0 || x_hit == width) && y_hit > height
end # function escaped

"""
    cumulative_probabilities(results)

this function build a vector of comulative fittnes results
"""
function cumulative_probabilities(results)
    cprob = Vector{Real}(undef,length(results))
    total = 0.0
    for (ii,res) in enumerate(results)
        total += res[2]
        cprob[ii] = total
    end
    return cprob
end # function cumulative_probabilities

"""
    selection(generation::Matrix,width::Real)

The function appends the y-coordinate where each ball hits the bag
edge to a list, The cumulative_probabilities function then stacks up the running
total of these, giving you your roulette wheel.
"""
function selection(generation::Matrix,width::Real)
    results = [hit_coordinate(θ,v,width) for (θ,v) in zip(generation[:,1],generation[:,2])]
    return cumulative_probabilities(results)
end # function selection

"""
    choose(choices)

The function spins the roulette wheel by picking a random number and
seeing which solution it corresponds to.
"""
function choose(choices)
    p = rand(Uniform(0,choices[end]))
    for (ind,val) in enumerate(choices)
        val >= p && return ind
    end
end # function choose

"""
crossover(generation::Matrix, Width::Real)

Armed with two parents, breed new solutions Not all genetic algorithms
use two parents; some just use one, and there are others that use more.
However, two is a conventional approach.
Notice this makes a new generation of the same size as the last generation.
"""
function crossover(generation::Matrix, width::Real)
    choices = selection(generation,width)
    next_generation = Matrix{Real}(undef,size(generation))
    for (ind,val) in enumerate(eachrow(generation))
        mum = generation[choose(choices),:]
        dad = generation[choose(choices),:] # TODO - why not taking out the pair choosen for mum?
        next_generation[ind,:] = collect(breed(mum,dad))
    end
    return next_generation
end # function crossover

"""
    breed(mum::Vector, dad::Vector)

In order to breed, the information is split—half from one parent and half from
another. Since there are two bits of information, each child will have a velocity
from one parent and an angle from the other. There is only one way to
split this.
In this example, the parents produce a single child.
"""
function breed(mum::Vector, dad::Vector)
    return mum[1],dad[2]
end # function breed

"""
    mutate(generation::Matrix; v_max::Real = 30, v_min::Real = 10)

the function change everything in a population, to keep it general.
the function uses probabilistic mutation, that is you want something to happen
on average, one time out of ten, One approac, is to get a random number between
0 and 1 and do that “something” if you get less than 0.1.
How do you mutate the value? It’s traditional to add or subtract a small
amount or scale (multiply) a little for real numbers. in our case
add a random number to the angle, but only use this if it stays between
0 and 180 degrees. and for the velocity, scale by something random between 0.9
and 1.1, cunningly avoiding the problem of potentially getting zero
"""
function mutate(generation::Matrix; v_max::Real = 30, v_min::Real = 10)
    for (ind,val) in enumerate(eachrow(generation))
        θ = val[1]
        v = val[2]
        if rand() < 0.1
            θ_new = θ + (rand(Uniform(-10,10)) |> deg2rad)
            if 0 < θ_new < π
                θ = θ_new
            end
        end
        if rand() < 0.1
            v *= rand(Uniform(0.9,1.1))
            if v > v_max
                v = v_max
            elseif v < v_min
                v = v_min
            end
        end
        generation[ind,:] = [θ,v]
    end
end # function mutate

"""
    display_epoch!(p::Plots.Plot, generation::Matrix, width::Real, height::Real)

this function plot the current epoch
"""
function display_epoch!(p::Plots.Plot, generation::Matrix, width::Real, height::Real)
    pcolor = rand(ColorSchemes.Accent_5.colors)
    markertype = rand(Plots.supported_markers())
    for (θ,v) in eachrow(generation)
        tend = (2*v*sin(θ))/9.81
        x,y = cannonball_trajectory(θ,v,width,range(0,tend[1],length=200))
        if escaped(θ, v, width, height)
            plot!(p, x, y, color = pcolor)
        else
            plot!(p, x, y, color = :red)
        end
        scatter!(p, x[1:Int(0.1*length(y)):end], y[1:Int(0.1*length(y)):end], marker = markertype, color = pcolor)
    end
end # function display_epoch

"""
    display_epoch(p::Plots.Plot, generation::Matrix, width::Real, height::Real)

this function plot the current epoch
"""
function display_epoch(pin::Plots.Plot, generation::Matrix, width::Real, height::Real)
    p = deepcopy(pin)
    for (θ,v) in eachrow(generation)
        tend = (2*v*sin(θ))/9.81
        x,y = cannonball_trajectory(θ,v,width,range(0,tend[1],length=200))
        if escaped(θ, v, width, height)
            plot!(p, x, y, color = :blue)
            scatter!(p, x[1:Int(0.1*length(y)):end], y[1:Int(0.1*length(y)):end], marker = :circle, color = :blue)
        else
            plot!(p, x, y, color = :red)
            scatter!(p, x[1:Int(0.1*length(y)):end], y[1:Int(0.1*length(y)):end], marker = :circle, color = :red)
        end
    end
    return p
end # function display_epoch

"""
    fire()

the main cannonball function
"""
function fire()
    epochs = 100
    items = 100
    height = 10
    width = 10
    max_cennonball_speed = 20
    box =( x = vec([0 0 width width]), y = vec([height 0 0 height]))
    generation = random_tries(items, max_speed = max_cennonball_speed)
    generation0 = deepcopy(generation) # save to contrast with last epoch
    # p = plot(box.x, box.y, legend = false, lw = 3, aspect_ratio=:equal)
    pbox = plot(box.x, box.y, legend = false, lw = 3)
    titles = ["First generation","Last generation"]
    for ii in 1:epochs
        results = []
        generation = crossover(generation, width)
        mutate(generation)
    end
    genPlot = (generation0,generation)
    pout = Vector{Any}(undef,length(genPlot))
    for (ind,gen) in enumerate(genPlot)
        pout[ind] = display_epoch(pbox, gen, width, height)
        title!(pout[ind],titles[ind])
        xlabel!(pout[ind],"X")
        ylabel!(pout[ind],"Y")
    end
    p = plot(pout[1],pout[2],layout = (2,1))
    display(p)
end # function fire


text:
# Autogenerated wrapper script for Libtask_jll for aarch64-linux-musl
export libtask_julia

JLLWrappers.@generate_wrapper_header("Libtask")
JLLWrappers.@declare_library_product(libtask_julia, "libtask_julia.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtask_julia,
        "lib/libtask_julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
$(TYPEDEF)
"""
struct LinearProblem{uType,isinplace,F,bType,P,K} <: AbstractLinearProblem{bType,isinplace}
    A::F
    b::bType
    u0::uType
    p::P
    kwargs::K
    @add_kwonly function LinearProblem{iip}(A,b,p=NullParameters();u0=nothing,
                                            kwargs...) where iip
        new{typeof(u0),iip,typeof(A),typeof(b),typeof(p),typeof(kwargs)}(
            A,b,u0,p,kwargs
        )
    end
end

function LinearProblem(A,b,args...;kwargs...)
    if A isa AbstractArray
        LinearProblem{true}(DiffEqArrayOperator(A),b,args...;kwargs...)
    else
        LinearProblem{isinplace(A, 4)}(A,b,args...;kwargs...)
    end
end

"""
$(TYPEDEF)
"""
struct NonlinearProblem{uType,isinplace,P,F,K} <: AbstractNonlinearProblem{uType,isinplace}
    f::F
    u0::uType
    p::P
    kwargs::K
    @add_kwonly function NonlinearProblem{iip}(f::AbstractNonlinearFunction{iip},u0,p=NullParameters();kwargs...) where iip
        new{typeof(u0),iip,typeof(p),typeof(f),typeof(kwargs)}(f,u0,p,kwargs)
    end

    """
    $(SIGNATURES)

    Define a steady state problem using the given function.
    `isinplace` optionally sets whether the function is inplace or not.
    This is determined automatically, but not inferred.
    """
    function NonlinearProblem{iip}(f,u0,p=NullParameters()) where iip
      NonlinearProblem(NonlinearFunction{iip}(f),u0,p)
    end
end

"""
$(SIGNATURES)

Define a steady state problem using an instance of
[`AbstractNonlinearFunction`](@ref AbstractNonlinearFunction).
"""
function NonlinearProblem(f::AbstractNonlinearFunction,u0,p=NullParameters();kwargs...)
  NonlinearProblem{isinplace(f)}(f,u0,p;kwargs...)
end

function NonlinearProblem(f,u0,p=NullParameters();kwargs...)
  NonlinearProblem(NonlinearFunction(f),u0,p;kwargs...)
end

"""
$(SIGNATURES)

Define a steady state problem from a standard ODE problem.
"""
NonlinearProblem(prob::AbstractNonlinearProblem) =
      NonlinearProblem{isinplace(prob)}(prob.f,prob.u0,prob.p)

"""
$(TYPEDEF)
"""
struct QuadratureProblem{isinplace,P,F,L,U,K} <: AbstractQuadratureProblem{isinplace}
    f::F
    lb::L
    ub::U
    nout::Int
    p::P
    batch::Int
    kwargs::K
    @add_kwonly function QuadratureProblem{iip}(f,lb,ub,p=NullParameters();
                                                nout=1,
                                                batch = 0, kwargs...) where iip
        new{iip,typeof(p),typeof(f),typeof(lb),
            typeof(ub),typeof(kwargs)}(f,lb,ub,nout,p,batch,kwargs)
    end
end

QuadratureProblem(f,lb,ub,args...;kwargs...) = QuadratureProblem{isinplace(f, 3)}(f,lb,ub,args...;kwargs...)

struct NoAD <: AbstractADType end

struct OptimizationFunction{iip,AD,F,G,H,HV,C,CJ,CH} <: AbstractOptimizationFunction{iip}
    f::F
    adtype::AD
    grad::G
    hess::H
    hv::HV
    cons::C
    cons_j::CJ
    cons_h::CH
end

(f::OptimizationFunction)(args...) = f.f(args...)

OptimizationFunction(args...;kwargs...) = OptimizationFunction{true}(args...;kwargs...)

function OptimizationFunction{iip}(f,adtype::AbstractADType=NoAD();
                     grad=nothing,hess=nothing,hv=nothing,
                     cons=nothing, cons_j=nothing,cons_h=nothing) where iip
    OptimizationFunction{iip,typeof(adtype),typeof(f),typeof(grad),typeof(hess),typeof(hv),
                         typeof(cons),typeof(cons_j),typeof(cons_h)}(
                         f,adtype,grad,hess,hv,cons,cons_j,cons_h)
end

"""
$(TYPEDEF)
"""
struct OptimizationProblem{iip,F,uType,P,B,LC,UC,S,K} <: AbstractOptimizationProblem{isinplace}
    f::F
    u0::uType
    p::P
    lb::B
    ub::B
    lcons::LC
    ucons::UC
    sense::S
    kwargs::K
    @add_kwonly function OptimizationProblem{iip}(f::OptimizationFunction{iip}, u0, p=NullParameters();
                                                  lb = nothing, ub = nothing,
                                                  lcons = nothing, ucons = nothing,
                                                  sense = nothing, kwargs...) where iip
        if xor(lb === nothing, ub === nothing)
            error("If any of `lb` or `ub` is provided, both must be provided.")
        end
        new{iip, typeof(f), typeof(u0), typeof(p),
            typeof(lb), typeof(lcons), typeof(ucons),
            typeof(sense), typeof(kwargs)}(f, u0, p, lb, ub, lcons, ucons, sense, kwargs)
    end
end

OptimizationProblem(f::OptimizationFunction,args...;kwargs...) = OptimizationProblem{isinplace(f)}(f,args...;kwargs...)
OptimizationProblem(f,args...;kwargs...) = OptimizationProblem{true}(OptimizationFunction{true}(f),args...;kwargs...)

isinplace(f::OptimizationFunction{iip}) where iip = iip
isinplace(f::OptimizationProblem{iip}) where iip = iip


text:
export filter_type

import Base: show

"""
    filter_type(type::Type{ T }) where { T }

Creates a `filter_type` operator, which filters items of the source Observable by emitting only
those that match a specified `T` type with a `<:` operator. This operator is a more efficient version of
`filter(v -> v <: T) |> map(T, v -> v)` operators chain.

# Producing

Stream of type `<: Subscribable{ T }` where `T` refers to `type` argument

# Examples
```jldoctest
using Rocket

source = from(Real[ 1, 2.0, 3, 4.0, 5, 6.0 ])
subscribe!(source |> filter_type(Int), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 3
[LogActor] Data: 5
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
filter_type(::Type{ T }) where { T }  = FilterTypeOperator{T}()

struct FilterTypeOperator{T} <: InferableOperator end

on_call!(::Type{T}, ::Type{T}, operator::FilterTypeOperator{T}, source) where { T    }      = source
on_call!(::Type{L}, ::Type{T}, operator::FilterTypeOperator{T}, source) where { L, T <: L } = proxy(T, source, FilterTypeProxy{L, T}())
on_call!(::Type{L}, ::Type{T}, operator::FilterTypeOperator{T}, source) where { L, T }      = error("filter_type(T) operator usage error. Cannot filter $(T) elements from stream of $(L) elements.")

operator_right(::FilterTypeOperator{T}, ::Type{T}) where { T }         = T
operator_right(::FilterTypeOperator{T}, ::Type{L}) where { L, T <: L } = T
operator_right(::FilterTypeOperator{T}, ::Type{L}) where { L, T }      = error("filter_type(T) operator creation error. Cannot filter $(T) elements from stream of $(L) elements.")

struct FilterTypeProxy{L, T} <: ActorProxy end

actor_proxy!(::Type{T}, proxy::FilterTypeProxy{L, T}, actor::A) where { L, A, T } = FilterTypeActor{L, A, T}(actor)

struct FilterTypeActor{L, A, T} <: Actor{L}
    actor :: A
end

on_next!(actor::FilterTypeActor{L, A, T}, data::T) where { L, A, T } = next!(actor.actor, data)
on_next!(::FilterTypeActor, data)                                    = nothing

on_error!(actor::FilterTypeActor, err) = error!(actor.actor, err)
on_complete!(actor::FilterTypeActor)   = complete!(actor.actor)

Base.show(io::IO, ::FilterTypeOperator)         = print(io, "FilterTypeOperator()")
Base.show(io::IO, ::FilterTypeProxy)            = print(io, "FilterTypeProxy()")
Base.show(io::IO, ::FilterTypeActor{L}) where L = print(io, "FilterTypeActor($L)")


text:
using DataFrames
using Grids
using PrettyTables

struct RecurrentRVPTest <: TrainingTest
    key::String
    state::AbstractEnvState
    state_next::AbstractEnvState
    image::AbstractImage
    action::AbstractEnvAction
    action_image::AbstractImage
    # @todo - "actions" should be named better
    # @todo - enable rvp
    # rvp::PolicyValueReward
    # @todo - should be "actions::Vector{<:AbstractEnvAction}"
    actions::Vector{<:Integer}
    value::Float32
    reward::Float32
end

function Base.show(io::IO, m::MIME"text/plain", x::RecurrentRVPTest)
    println("\nRECURRENT RVP TEST #$(x.key)")
    println("\nTo Play:\n$(x.state.to_play.name)")
    println("\nAction:\n$(x.action.idx)")
    println("\nReward:\n$(x.reward)")
    println("\nNext To Play:\n$(x.state_next.to_play.name)")
    println("\nNext Optimal Action(s):\n$(join(map(y->y.idx, x.actions), ","))")
    println("\nNext Value:\n$(x.value)")
end

struct RecurrentRVPTestSVG
    test::RecurrentRVPTest
end
svg(x::RecurrentRVPTest) = RecurrentRVPTestSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::RecurrentRVPTestSVG)
    test = x.test
    s = test.state
    a = test.action
    r = test.reward
    s_next = test.state_next
    actions = test.actions
    state_grid = getgrid(
        s;
        header="to play: $(s.to_play.name)",
        footer="action: $(a.idx)",
        margin_left=0
    )
    state_next_grid = getgrid(
        s_next;
        header="to play: $(s_next.to_play.name)",
        margin_left=0
    )
    next_states = Rect[]
    from = 1
    len = length(actions)
    to = min(len, 3)
    for (idx,a) in enumerate(actions[from:to])
        push!(
            next_states,
            getgrid(s_next, a, r, getnextstate(s_next, a); header="action: $(a.idx)", margin_left=0)
        )
    end
    if length(actions) > 3
        push!(next_states, Grid(1,1, (r,c)->"..."; line_thickness=0))
    end
    hiddenstate_grid = Grid(
        [
            getgrid(
                s,
                a,
                0f0,
                s_next;
                header="to play: $(s_next.to_play.name)",
                margin_left=0
            )
        ],
        :,
        1,
        header="HIDDEN STATE",
        margin_left=0,
        line_thickness=0
    )
    policy_grid_header = "ACTION TEST"
    policy_grid = Grid(
        next_states,
        :,
        1,
        header=(text=policy_grid_header, height=40),
        margin_left=0,
        line_thickness=0,
        bg_color=RandColor()
    )
    g = Grid(
        [
            hiddenstate_grid,
            policy_grid
        ];
        cell_valign=:top,
        header="RECURRENT RVP TEST #$(test.key)",
        line_thickness=0
    )
    show(io, m, g)
end


text:
@testset "lessorqual.jl" begin

    @testset "isLessOrEqual()" begin
        trailer = SeaPearl.Trailer()
        b = SeaPearl.BoolVar("b", trailer)
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

        @test constraint in b.onDomainChange
        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint.active.value
    end
    @testset "propagate!(::isLessOrEqual)" begin
        @testset "boolean true" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.IntVar(2, 6, "x", trailer)
            y = SeaPearl.IntVar(2, 3, "y", trailer)

            constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, true)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

            @test length(x.domain) == 2
            @test !(5 in x.domain)
            @test 3 in x.domain
            @test prunedDomains == SeaPearl.CPModification("x" => [4, 5, 6])
            @test constraint.active.value

            z = SeaPearl.IntVar(1, 1, "z", trailer)
            constraint2 = SeaPearl.isLessOrEqual(b, x, z, trailer)

            @test !SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        end

        @testset "boolean false" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.IntVar(1, 4, "x", trailer)
            y = SeaPearl.IntVar(2, 6, "y", trailer)

            constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            SeaPearl.assign!(b, false)

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

            @test length(y.domain) == 2
            @test length(x.domain) == 2
            @test !(4 in y.domain)
            @test !(5 in y.domain)
            @test 3 in y.domain
            @test !(1 in x.domain)
            @test !(2 in x.domain)
            @test 3 in x.domain
            @test prunedDomains == SeaPearl.CPModification("y" => [4, 5, 6], "x" => [1, 2])
            @test constraint.active.value

            z = SeaPearl.IntVar(4, 4, "z", trailer)
            constraint2 = SeaPearl.isLessOrEqual(b, x, z, trailer)

            @test !SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        end

        @testset "boolean unassigned - maxx <= miny" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.IntVar(2, 4, "x", trailer)
            y = SeaPearl.IntVar(4, 6, "y", trailer)

            constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

            @test length(x.domain) == 3
            @test length(y.domain) == 3
            @test SeaPearl.assignedValue(b)
            @test prunedDomains == SeaPearl.CPModification("b" => [false])
            @test !constraint.active.value

        end

        @testset "boolean unassigned - minx > maxy" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.IntVar(5, 7, "x", trailer)
            y = SeaPearl.IntVar(2, 4, "y", trailer)

            constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

            @test length(x.domain) == 3
            @test length(y.domain) == 3
            @test !SeaPearl.assignedValue(b)
            @test prunedDomains == SeaPearl.CPModification("b" => [true])
            @test !constraint.active.value

        end

        @testset "boolean unassigned - casual x & y" begin

            trailer = SeaPearl.Trailer()
            b = SeaPearl.BoolVar("b", trailer)
            x = SeaPearl.IntVar(2, 6, "x", trailer)
            y = SeaPearl.IntVar(3, 7, "y", trailer)

            constraint = SeaPearl.isLessOrEqual(b, x, y, trailer)

            toPropagate = Set{SeaPearl.Constraint}()
            prunedDomains = SeaPearl.CPModification()

            @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

            @test length(x.domain) == 5
            @test length(y.domain) == 5
            @test length(b.domain) == 2
            @test prunedDomains == SeaPearl.CPModification()
            @test constraint.active.value

        end

    end
end


text:
# Autogenerated wrapper script for Jansson_jll for aarch64-linux-gnu
export libjansson

JLLWrappers.@generate_wrapper_header("Jansson")
JLLWrappers.@declare_library_product(libjansson, "libjansson.so.4")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libjansson,
        "lib/libjansson.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Base.Test, HDF5, JLD
using StateSpaceRoutines
using QuantEcon: solve_discrete_lyapunov

my_tests = [
            "kalman_filter"
            "chand_recursion"
            "tempered_particle_filter"
            "smoothers"
            ]

for test in my_tests
    @testset "$test.jl" begin
        include(string("$test.jl"))
    end
end


text:
using Weave

Weave.set_chunk_defaults!(:echo => false, :results => "raw")

include("../src/lotka-volterra-2d.jl")

weave("weave/lotka-volterra-2d-spark-glvprk.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-glspark.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-lobspark.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-pinternal.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-pmidpoint.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-psymmetric.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("weave/lotka-volterra-2d-spark-plobatto.jmd",
         out_path = "build/lotka-volterra-2d-spark",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)


text:
export measure,
       contract,
       norm,
       ent_entropy,
       enlarge_bond_dimension!,
       svd_truncate!,
       simplify!,
       save_mps,
       read_mps

"""
    expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int) where T<:Number

Compute expected value of the operator op_i in the state psi.
"""
function expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int) where T<:Number
    # Convert op_i to have the same type as psi and reshape so that it can pass
    # in prop_right3.
    c_op_i = convert.(T, reshape(op_i, 1, size(op_i, 1), size(op_i, 2), 1))

    L = ones(T, 1, 1)
    for s=1:psi.L
        if s != i
            L = prop_right2(L, psi.M[s], psi.M[s])
        else
            # Reshape L to a rank 3 tensor so that it can pass in prop_right3
            # and then reshape it back to rank 2.
            L3 = reshape(L, size(L, 1), 1, size(L, 2))
            L3 = prop_right3(L3, psi.M[s], c_op_i, psi.M[s])
            L = reshape(L3, size(L3, 1), size(L3, 3))
        end
    end
    return L[1, 1]
end

"""
    expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int) where T<:Number

Compute expected value of the operator (str_)op_i in the state psi.
"""
function expected(psi::Mps{T}, str_op_i::String, i::Int) where T<:Number
    return expected(psi, str_to_op(str_op_i), i)
end

"""
    expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int,
             op_j::AbstractMatrix{<:Number}, j::Int;
             ferm_op::AbstractMatrix{<:Number}=Matrix{T}(I, psi.d, psi.d)
             ) where T<:Number

Compute expected value of the operator op_i*op_j in the state psi. An optional
fermionic parity operator can be passed in ferm_op.
"""
function expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int,
                  op_j::AbstractMatrix{<:Number}, j::Int;
                  ferm_op::AbstractMatrix{<:Number}=Matrix{T}(I, psi.d, psi.d)
                  ) where T<:Number
    # Multiply operators acting on same site and call single operator method.
    if i == j
        return expected(psi, op_j*op_i, i)
    end

    # Convert op_i, op_j, and ferm_op to have the same type as psi and reshape
    # to a rank 4 tensor so that they can pass in prop_right3.
    c_op_i = convert.(T, reshape(op_i, 1, size(op_i, 1), size(op_i, 2), 1))
    c_op_j = convert.(T, reshape(op_j, 1, size(op_j, 1), size(op_j, 2), 1))
    c_ferm_op = convert.(T, reshape(ferm_op, 1, size(ferm_op, 1),
                                    size(ferm_op, 2), 1))

    # Sort operators.
    s1 = min(i, j)
    op_1 = i < j ? c_op_i : c_op_j
    s2 = max(i, j)
    op_2 = i > j ? c_op_i : c_op_j

    L = ones(T, 1, 1)
    for s=1:s1-1
        L = prop_right2(L, psi.M[s], psi.M[s])
    end
    # Reshape L to a rank 3 tensor so that it can pass in prop_right3.
    L3 = reshape(L, size(L, 1), 1, size(L, 2))
    L3 = prop_right3(L3, psi.M[s1], op_1, psi.M[s1])
    # Between min(i, j) and max(i, j) apply the fermionic parity operator.
    for s=s1+1:s2-1
        L3 = prop_right3(L3, psi.M[s], c_ferm_op, psi.M[s])
    end
    L3 = prop_right3(L3, psi.M[s2], op_2, psi.M[s2])
    # Reshape L3 back to a rank 2 tensor.
    L = reshape(L3, size(L3, 1), size(L3, 3))
    for s=s2+1:psi.L
        L = prop_right2(L, psi.M[s], psi.M[s])
    end

    return L[1, 1]
end

"""
    expected(psi::Mps{T}, op_i::AbstractMatrix{<:Number}, i::Int,
             op_j::AbstractMatrix{<:Number}, j::Int;
             ferm_op::AbstractMatrix{<:Number}=Matrix{T}(I, psi.d, psi.d)
             ) where T<:Number

Compute expected value of the operator op_i*op_j in the state psi. An optional
fermionic parity operator can be passed in ferm_op.
"""
function expected(psi::Mps{T}, op_i::String, i::Int, op_j::String, j::Int;
                  ferm_op::String="none") where T<:Number
    if ferm_op != "none"
        return expected(psi, str_to_op(op_i), i, str_to_op(op_j), j,
                        ferm_op=str_to_op(ferm_op))
   else
        return expected(psi, str_to_op(op_i), i, str_to_op(op_j), j)
   end
end

"""
    contract(psi::Mps{T}, phi::Mps{T}) where T<:Number

Contraction of two MPS: <psi|phi>.
"""
function contract(psi::Mps{T}, phi::Mps{T}) where T<:Number
    L = Matrix{T}(I, 1, 1)
    for i=1:psi.L
        L = prop_right2(L, phi.M[i], psi.M[i])
    end
    return L[1, 1]
end

"""
    norm(psi::Mps{<:Number})

Norm of a MPS: <psi|psi>. Extend the LinearAlgebra module function.
"""
function LinearAlgebra.norm(psi::Mps{<:Number})
    return real(contract(psi, psi))
end

"""
    schmidt_decomp(psi::Mps{T}, i::Int) where T<:Number

Compute the singular values of the Schmidt decomposition of state `psi` between
sites `i` and `i+1`.
"""
function schmidt_decomp(psi::Mps{T}, i::Int) where T<:Number
    # Write the tensors 1 -> i in left canonical form and the tensors i+1 -> L
    # in right canonical form.
    # Left canonical form.
    L = ones(T, 1, 1)
    for j=1:i
        L = absorb_fromleft(L, psi.M[j])
        _, L = factorize_qr(L)
    end
    R = ones(T, 1, 1)
    for j=reverse(i+1:psi.L)
        R = absorb_fromright(psi.M[j], R)
        R, _ = factorize_lq(R)
    end
    return svdvals!(L*R)
end

"""
    ent_entropy(psi::Mps{T}, i::Int) where T<:Number

Compute the entanglement entropy of state `psi` between sites `i` and `i+1`.
"""
function ent_entropy(psi::Mps{T}, i::Int) where T<:Number
    rho = schmidt_decomp(psi, i)
    S = 0.
    for j=1:length(rho)
        abs(rho[j]) < 1e-10 && continue
        S -= rho[j]*log2(rho[j])
    end
    return S
end

"""
    bond_dimension_with_m(L::Int, i::Int, m::Int, d::Int)

Compute the appropriate bond dimension at site `i` for an Mps of length `L`,
physical dimension `d`, and maximum bond dimension `m`.

There are L+1 possible bonds: the bond at i=1 has dimension 1 and so does the
bond at i=L+1. The bond at site i=2 has dimension d, like the bond at i=L.
"""
function bond_dimension_with_m(L::Int, i::Int, m::Int, d::Int)
    mi = m
    if i <= L>>1
        if (i-1)*log(d) < log(m)
            mi = d^(i-1)
        end
    else
        if (L+1-i)*log(d) < log(m)
            mi = d^(L+1-i)
        end
    end
    return mi
end

"""
    enlarge_bond_dimension!(psi::Mps{T}, max_D::Int) where T<:Number

Take a state and add 0's to the tensors until a maximum bond dimension is
reached.
"""
function enlarge_bond_dimension!(psi::Mps{T}, max_D::Int) where T<:Number
    L = psi.L

    for i=1:L
        m1 = size(psi.M[i], 1)
        m2 = size(psi.M[i], 3)
        # Check if m1 and m2 need to be resized.
        new_m1 = bond_dimension_with_m(psi.L, i, max_D, psi.d)
        new_m2 = bond_dimension_with_m(psi.L, i+1, max_D, psi.d)
        need_resize_m1 = m1 < new_m1
        need_resize_m2 = m2 < new_m2
        # Resize M with the appropriate dimensions, if needed.
        if need_resize_m1 || need_resize_m2
            new_M = zeros(T, new_m1, psi.d, new_m2)
            new_M[1:m1, :, 1:m2] = psi.M[i]
            psi.M[i] = new_M
        end
    end
    return psi
end

"""
    svd_truncate!(psi::Mps{T}, m::Int) where T<:Number

Truncate the bond dimension of psi to a maximum of m using SVD decomposition.
"""
function svd_truncate!(psi::Mps{T}, m::Int) where T<:Number
    US = ones(T, 1, 1)
    for i=reverse(2:psi.L)
        psi.M[i] = absorb_fromright(psi.M[i], US)
        US, psi.M[i] = factorize_svd_left(psi.M[i],
                                          cutoff=0.,
                                          dimcutoff=m,
                                          normalize_S=false)
    end
    # Contract and normalize last tensor.
    psi.M[1] = absorb_fromright(psi.M[1], US)
    psi.M[1] ./= norm(psi.M[1])
    return psi
end

"""
    simplify!(psi::Mps{T}, m::Int;
              max_sweeps::Int=500, tol::Float64=1e-6) where T<:Number

Variationally simplify psi so that its maximum bond dimension is m.
"""
function simplify!(psi::Mps{T}, m::Int;
                   max_sweeps::Int=500, tol::Float64=1e-6) where T<:Number
    # Normalize state if it is not normalized.
    norm_psi = norm(psi)
    psi.M[end] ./= sqrt(norm_psi)

    # Copy original state to make it the objective state to which we want to
    # approximate. Make it a new Mps.
    obj_M = deepcopy(psi.M)
    obj = Mps(obj_M, psi.L, psi.d)
    # Reduce the bond dimension of the state to the wanted one.
    svd_truncate!(psi, m)

    # Create left and right environments.
    Le = fill(ones(T, 1, 1), psi.L)
    Re = fill(ones(T, 1, 1), psi.L)
    # Initialize left environment.
    for i=2:psi.L
        Le[i] = prop_right2(Le[i-1], psi.M[i-1], obj.M[i-1])
    end

    # Compute distance of new state to the objective one.
    dist = [2. - 2*real(contract(psi, obj))]
    dist_converged = false
    nsweep = 2
    while dist[nsweep-1] > tol && nsweep <= max_sweeps+1 && !dist_converged
        # Do left and right sweeps.
        do_sweep_simplify!(psi, obj, Le, Re, -1)
        do_sweep_simplify!(psi, obj, Le, Re, +1)
        # Compute distance of `psi` to objective after sweeps. Use abs to
        # remove any phase that could enter in `psi`.
        push!(dist, 2. - 2*abs(real(contract(psi, obj))))
        # Update while loop control parameters.
        dist_converged = abs(dist[nsweep] - dist[nsweep-1]) < 1e-8
        nsweep += 1
    end

    return psi
end

"""
    do_sweep_simplify!(psi::Mps{T}, obj::Mps{T},
                       Le::Vector{Array{T, 2}}, Re::Vector{Array{T, 2}},
                       sense::Int) where T<:Number

Do a sweep to variationally simplify `psi`.
"""
function do_sweep_simplify!(psi::Mps{T}, obj::Mps{T},
                            Le::Vector{Array{T, 2}}, Re::Vector{Array{T, 2}},
                            sense::Int) where T<:Number

    sense == 1 || sense == -1 || throw("`Sense` must be either `+1` or `-1`.")
    # Order of sites to do the sweep.
    sweep_sites = sense == +1 ? (1:psi.L) : reverse(1:psi.L)
    for i in sweep_sites
        # Compute local tensor.
        @tensor Mi[l1, s, r1] := Le[i][l1, l2]*obj.M[i][l2, s, r2]*Re[i][r1, r2]
        Mi = conj(Mi)

        # Update left and right environments.
        if sense == +1
            Mi = reshape(Mi, (size(Le[i], 1)*psi.d, size(Re[i], 1)))
            Qa, Ra = qr(Mi)
            Qa = Matrix(Qa)
            Qa = reshape(Qa, (size(Le[i], 1), psi.d, size(Qa, 2)))

            # Update left and right environments at L[i+1] and R[i] and psi.
            psi.M[i] = Qa
            if i < psi.L
                Le[i+1] = prop_right2(Le[i], Qa, obj.M[i])
                Re[i] = Ra*Re[i]
            end
        else
            Mi = reshape(Mi, (size(Le[i], 1), psi.d*size(Re[i], 1)))
            La, Qa = lq(Mi)
            Qa = Matrix(Qa)
            Qa = reshape(Qa, (size(Qa, 1), psi.d, size(Re[i], 1)))

            # Update left and right environments at L[i] and R[i-1] and psi.
            psi.M[i] = Qa
            if i > 1
                Re[i-1] = prop_left2(Qa, obj.M[i], Re[i])
                Le[i] = transpose(La)*Le[i]
            end
        end
    end
    return psi
end

#
# INPUT/OUTPUT OF MPS.
#

"""
    save_mps(filename::String, psi::Mps{T}) where T<:Number

Save the state `psi` in file `filename` in HDF5 format:
- `L`: length of Mps.
- `d`: physical dimension.
- `T`: type of Mps tensor elements.
- `Ai`: left canonical tensors.

More info about the HDF5 format can be found here:
    https://github.com/JuliaIO/HDF5.jl/blob/master/doc/hdf5.md
"""
function save_mps(filename::String, psi::Mps{T}) where T<:Number
    h5write(filename, "L", psi.L)
    h5write(filename, "d", psi.d)
    h5write(filename, "T", "$T")
    is_complex = !(T <: Real)
    for i=1:psi.L
        if !is_complex
            h5write(filename, "M$i", psi.M[i])
        else
            h5write(filename, "realM_$i", real(psi.M[i]))
            h5write(filename, "imagM_$i", imag(psi.M[i]))
        end
    end
    return
end

"""
    read_mps(filename::String)

Read the state `psi` in stored in file `filename` inHDF5 format.
"""
function read_mps(filename::String)
    L = h5read(filename, "L")
    d = h5read(filename, "d")
    str_T = h5read(filename, "T")
    if str_T == "Float32"
        T = Float32
    elseif str_T == "Float64"
        T = Float64
    elseif str_T == "Complex{Float32}"
        T = ComplexF32
    elseif str_T == "Complex{Float64}"
        T = ComplexF64
    end
    is_complex = !(T <: Real)

    # Start Mps as full and then replace tensors.
    psi = Mps(T, L, "full")
    for i=1:psi.L
        if !is_complex
            psi.M[i] = h5read(filename, "M$i")
        else
            psi.M[i] = complex.(h5read(filename, "realM_$i"),
                                h5read(filename, "imagM_$i"))
        end
    end

    return psi
end


text:
# This example shows how to use different predictors for symbols prediction

# Make the package available
using DiscretePredictors

# Some psuedo data
data = [ 'a', 'a', 'a', 'b', 'a',
         'b', 'b', 'b', 'b', 'b',
         'a', 'a', 'b', 'c', 'c',
         'd', 'd', 'c', 'b', 'a',
         'a', 'a', 'a' ];

# Make a KOM predictor for 'Char' symbol
p   = KOM{Char}(3)

# Just to show some usages
println( "I'm a ", info_string(p), " predictor" )
println( "You can uniquely identify as ", unique_string(p) )

# Add first symbol
add!( p, data[1] )
corr_pred = 0   # Create a counter for correct predictions
for idx = 2:length(data)
    # Predict a symbol from model
    prediction = get_best_symbol( p )
    # Compare
    corr_pred += (prediction==data[idx])?1:0
    # Add next symbol
    add!( p, data[idx] )
end

# Print out accuracy
println( "Accuracy of ", info_string(p), " : ", corr_pred/(length(data)-1) )

# Get current prediction table
println( "Predictions: ", predict(p) )

# Print out size of model
println( "Size of ", info_string(p), ": ", size(p), " nodes" )

# Print out prediction model
println( "Model: \n", p.model )


text:
#=
1. Create set of N triangles T such that
	(a) they all meet in a common vertex v.
	(b) the edges E opposite to v are parallel to each other with uniform distance between
	pairs of triangles.
	(c) the length of each edge e in E equals the largest distance between to pairs of edges 		in E.
2. Create a copy T′ of T.
3. Rotate T′ around the axis through v and the center of the edges E by 90 degrees.
4. Measure time for T ∪ T′.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;