text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra

# testtype = "PathDependent1D"
testtype = "Plasticity1D"

# prop = Dict("name"=> testtype, "rho"=> 0.1, "E"=> 200, "B"=> 10.0,
#             "sigmaY"=>0.300, "K"=>1/9*200, "A0"=> 1.0, "eta"=> 10.0)

prop = Dict("name"=> testtype, "rho"=> 8000.0, "E"=> 200e9, "nu"=> 0.45,
            "sigmaY"=>0.3e9, "K"=>1/9*200e9, "B"=> 0.0, "A0"=> 0.005)

include("NNTrussPull_Domain.jl")

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)

# need to update state in domain from globdat
updateStates!(domain, globdat)

for i = 1:NT
    solver = NewmarkSolver(Δt, globdat, domain, -1.0, 0.0, 1e-5, 1e-5, 100) # ok
    # close("all")
    # scatter(nodes[:, 1], nodes[:,2], color="red")
    # u,v = domain.state[1:domain.nnodes], domain.state[domain.nnodes+1:end]
    # scatter(nodes[:, 1] + u, nodes[:,2] + v, color="blue")
    # savefig("$(@__DIR__)/Debug/$i.png")

end
# error()
# error()

if !isdir("$(@__DIR__)/Data")
    mkdir("$(@__DIR__)/Data")
end
if !isdir("$(@__DIR__)/Debug")
    mkdir("$(@__DIR__)/Debug")
end

write_data("$(@__DIR__)/Data/$tid.dat", domain)
@save "Data/domain$tid.jld2" domain
# close("all")
# scatter(nodes[:, 1], nodes[:,2], color="red")
# u,v = domain.state[1:domain.nnodes], domain.state[domain.nnodes+1:end]
# scatter(nodes[:, 1] + u, nodes[:,2] + v, color="blue")
# ylim(-0.2,1.5)
# savefig("$(@__DIR__)/Debug/terminal$tid.png")


text:
using Test
using DataFrames
using Random: randstring
using WeakRefStrings: StringVector

@testset "JDF.jl eachcol" begin
    df = DataFrame([collect(1:100) for i = 1:3000])
    df[!, :int_missing] =
        rand([rand(rand([UInt, Int, Float64, Float32, Bool])), missing], size(df, 1))

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i = 1:size(df, 1)]
    df[!, :stringarray] = StringVector([randstring(8) for i = 1:size(df, 1)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i = 1:size(df, 1)]
    df[!, :stringarray_missing] =
        StringVector([rand([missing, randstring(8)]) for i = 1:size(df, 1)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i = 1:size(df, 1)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:size(df, 1), 10), :char_missing] .= missing

    JDF.save("a.jdf", df)

    df2 = jdf"a.jdf"

    @test size(df2, 2) == 3009
    @test size(df2, 1) == 100

    df3 = [a for a in eachcol(df2)]

    df4 = DataFrame(df3)

    @test size(df4) == size(df)
    @test all([isequal(df4[!, n], df[!, n]) for n = 1:size(df4, 2)])

    # clean up
    rm("a.jdf", force = true, recursive = true)
end

ext: jl


text:
#A counter type

"""
    Accumulator{T, V<:Number}

A accumulator is a data structure that maintains an accumulated total for each key.
The particular case where those totals are integers is a counter.
"""
struct Accumulator{T, V<:Number} <: AbstractDict{T,V}
    map::Dict{T,V}
end

## constructors

Accumulator{T, V}() where {T,V<:Number} = Accumulator{T,V}(Dict{T,V}())
Accumulator(map::AbstractDict) = Accumulator(Dict(map))
Accumulator(ps::Pair...) = Accumulator(Dict(ps))

counter(T::Type) = Accumulator{T,Int}()
counter(dct::AbstractDict{T, V}) where {T, V<:Integer} = Accumulator{T, V}(Dict(dct))

"""
    counter(seq)

Returns an `Accumulator` object containing the elements from `seq`.
"""
function counter(seq)
    ct = counter(eltype_for_accumulator(seq))
    for x in seq
        inc!(ct, x)
    end
    return ct
end

eltype_for_accumulator(seq::T) where T = eltype(T)
function eltype_for_accumulator(seq::T) where {T<:Base.Generator}
    Base.@default_eltype(seq)
end

copy(ct::Accumulator) = Accumulator(copy(ct.map))

length(a::Accumulator) = length(a.map)

## retrieval

get(ct::Accumulator, x, default) = get(ct.map, x, default)
# need to allow user specified default in order to
# correctly implement "informal" AbstractDict interface

getindex(ct::Accumulator{T,V}, x) where {T,V} = get(ct.map, x, zero(V))

setindex!(ct::Accumulator, x, v) = setindex!(ct.map, x, v)

haskey(ct::Accumulator, x) = haskey(ct.map, x)

keys(ct::Accumulator) = keys(ct.map)

values(ct::Accumulator) = values(ct.map)

sum(ct::Accumulator) = sum(values(ct.map))

## iteration

iterate(ct::Accumulator, s...) = iterate(ct.map, s...)

# manipulation

"""
    inc!(ct::Accumulator, x, [v=1])

Increments the count for `x` by `v` (defaulting to one)
"""
inc!(ct::Accumulator, x, v::Number) = (ct[x] += v)
inc!(ct::Accumulator{T,V}, x) where {T,V} = inc!(ct, x, one(V))

# inc! is preferred over push!, but we need to provide push! for the Bag interpreation
# which is used by classified_collections.jl
push!(ct::Accumulator, x) = inc!(ct, x)
push!(ct::Accumulator, x, a::Number) = inc!(ct, x, a)

# To remove ambiguities related to Accumulator now being a subtype of AbstractDict
push!(ct::Accumulator, x::Pair)  = inc!(ct, x)

"""
    dec!(ct::Accumulator, x, [v=1])

Decrements the count for `x` by `v` (defaulting to one)
"""
dec!(ct::Accumulator, x, v::Number) = (ct[x] -= v)
dec!(ct::Accumulator{T,V}, x) where {T,V} = dec!(ct, x, one(V))

#TODO: once we are done deprecating `pop!` for `reset!` then add `pop!` as an alias for `dec!`

"""
    merge!(ct1::Accumulator, others...)

Merges the other counters into `ctl`,
summing the counts for all elements.
"""
function merge!(ct::Accumulator, other::Accumulator)
    for (x, v) in other
        inc!(ct, x, v)
    end
    ct
end

function merge!(ct1::Accumulator, others::Accumulator...)
    for ct in others
        merge!(ct1,ct)
    end
    return ct1
end

"""
     merge(counters...)

Creates a new counter with total counts equal to the sum of the counts in the counters given as arguments.

See also merge!
"""
function merge(ct1::Accumulator, others::Accumulator...)
    ct = copy(ct1)
    merge!(ct,others...)
end

"""
    reset!(ct::Accumulator, x)

Remove a key `x` from an accumulator, and return its current value
"""
reset!(ct::Accumulator{<:Any,V}, x) where V = haskey(ct.map, x) ? pop!(ct.map, x) : zero(V)

"""
     nlargest(acc::Accumulator, [n])

Returns a sorted vector of the `n` most common elements, with their counts.
If `n` is omitted, the full sorted collection is returned.

This corresponds to Python's `Counter.most_common` function.

Example
```
julia> nlargest(counter("abbbccddddda"))

4-element Array{Pair{Char,Int64},1}:
 'd'=>5
 'b'=>3
 'c'=>2
 'a'=>2

julia> nlargest(counter("abbbccddddda"),2)

2-element Array{Pair{Char,Int64},1}:
 'd'=>5
 'b'=>3

```
"""
nlargest(acc::Accumulator) = sort!(collect(acc), by=last, rev=true)
nlargest(acc::Accumulator, n) = partialsort!(collect(acc), 1:n, by=last, rev=true)

"""
     nsmallest(acc::Accumulator, [n])

Returns a sorted vector of the `n` least common elements, with their counts.
If `n` is omitted, the full sorted collection is returned.

This is the opposite of the `nlargest` function.
For obvious reasons this will not include zero counts for items not encountered.
(unless those elements are added to he accumulator directly, eg via `acc[foo]=0)
"""
nsmallest(acc::Accumulator) = sort!(collect(acc), by=last, rev=false)
nsmallest(acc::Accumulator, n) = partialsort!(collect(acc), 1:n, by=last, rev=false)

###########################################################
## Multiset operations

struct MultiplicityException{K,V} <: Exception
    k::K
    v::V
end

function Base.showerror(io::IO, err::MultiplicityException)
    print(io, "When using an `Accumulator` as a multiset, all elements must have positive multiplicity")
    print(io, " element `$(err.k)` has multiplicity $(err.v)")
end

drop_nonpositive!(a::Accumulator, k) = (a[k] > 0 || delete!(a.map, k))

function Base.setdiff(a::Accumulator, b::Accumulator)
    ret = copy(a)
    for (k, v) in b
        v > 0 || throw(MultiplicityException(k, v))
        dec!(ret, k, v)
        drop_nonpositive!(ret, k)
    end
    return ret
end

Base.issubset(a::Accumulator, b::Accumulator) = all(b[k] >= v for (k, v) in a)

Base.union(a::Accumulator, b::Accumulator, c::Accumulator...) = union(union(a,b), c...)
Base.union(a::Accumulator, b::Accumulator) = union!(copy(a), b)
function Base.union!(a::Accumulator, b::Accumulator)
    for (kb, vb) in b
        va = a[kb]
        vb >= 0 || throw(MultiplicityException(kb, vb))
        va >= 0 || throw(MultiplicityException(kb, va))
        a[kb] = max(va, vb)
    end
    return a
end

Base.intersect(a::Accumulator, b::Accumulator, c::Accumulator...) = insersect(intersect(a,b), c...)
Base.intersect(a::Accumulator, b::Accumulator) = intersect!(copy(a), b)
function Base.intersect!(a::Accumulator, b::Accumulator)
    for k in union(keys(a), keys(b)) # union not interection as we want to check both multiplicties
        va = a[k]
        vb = b[k]
        va >= 0 || throw(MultiplicityException(k, va))
        vb >= 0 || throw(MultiplicityException(k, vb))

        a[k] = min(va, vb)
        drop_nonpositive!(a, k) # Drop any that ended up zero
    end
    return a
end
function Base.show(io::IO, acc::Accumulator{T,V}) where {T,V}
    l = length(acc)
    if l>0
        print(io, "Accumulator(")
    else
        print(io,"Accumulator{$T,$V}(")
    end
    for (count, (k, v)) in enumerate(acc)
        print(io, k, " => ", v)
        if count < l
            print(io, ", ")
        end
    end
    print(io, ")")
end


text:
using Base.Test
using RNGTest

include("tests.jl")


text:
# Script to replicate figure 8 and 9

using DataFrames, XLSX, StatsPlots

# -------------- Import data ---------------------------------------------------------------------------------------------------------

t_fig_9_10 = DataFrame(XLSX.readtable("clean/t_tab2.xlsx", "Sheet1")...) # still a lot of NaN
transform!(t_fig_9_10, :country .=> ByRow(string) .=> :country, renamecols=false) # change type to String
transform!(t_fig_9_10, [:year, :sector] .=> ByRow(Int) .=> [:year, :sector], renamecols=false) # change type to Int64
transform!(t_fig_9_10, [:FD_GO, :VA_GO, :U, :D] .=> ByRow(Float64) .=> [:FD_GO, :VA_GO, :U, :D], renamecols=false) # change type to Float64

subset!(t_fig_9_10, :D => ByRow(x -> !isnan(x))) # take out NaN
transform!(t_fig_9_10, :sector => ByRow(x -> ifelse(x in goods_sectors, "goods", "services")) => :sector_dummy, renamecols=false) # sector dummy for grouping

# -------------- Figure 9 --------------------------------------------------------------------------------------------------------------------------------

function plot_fig_9(df::DataFrame, year::Int64, industry::String, a, b) # input vectors in so we can plot multiple countries for comparison

    subset!(df, :year => ByRow(x -> x == year))
    df = industry in ["goods", "services"] ? subset(df, :sector_dummy => ByRow(x -> x == industry)) : df

    p = @df df scatter(cols(a), cols(b), group=:sector_dummy, xlims=(-0.1,1.1), ylims=(-0.1,1.1),
        xlabel="$(String(a))", ylabel="$(String(b))", label=:none, alpha=0.65,
        title="$(industry)", smooth=:true, linecolor=:black, lw=3)

    return p
end

l = @layout [a b c] # 1×3 layout
p1 = plot_fig_9(t_fig_9_10, 2011, "goods", :VA_GO, :FD_GO)
p2 = plot_fig_9(t_fig_9_10, 2011, "services", :VA_GO, :FD_GO)
p3 = plot_fig_9(t_fig_9_10, 2011, "all industries", :VA_GO, :FD_GO)
title = "Figure 9: Correlation between FD and VA: Goods vs Services"
p_fig9 = plot(p1, p2, p3, layout=l, plot_title=title, plot_titlefontsize=10)
savefig(p_fig9, "images/figure9.png") # export image to folder

# -------------- Figure 10 --------------------------------------------------------------------------------------------------------------------------------

function plot_fig_10(df::DataFrame, year::Int64, industry::String, a, b) # input vectors in so we can plot multiple countries for comparison

    subset!(df, :year => ByRow(x -> x == year))
    df = industry in ["goods", "services"] ? subset(df, :sector_dummy => ByRow(x -> x == industry)) : df

    p = @df df scatter(cols(a), cols(b), group=:sector_dummy, xlims=(0.9,4.5), ylims=(0.9,4.5),
        xlabel="$(String(a))", ylabel="$(String(b))", label=:none, alpha=0.65,
        title="$(industry)", smooth=:true, linecolor=:black, lw=3)

    return p
end

l = @layout [a b c] # 1×3 layout
p1 = plot_fig_10(t_fig_9_10, 2011, "goods", :D, :U)
p2 = plot_fig_10(t_fig_9_10, 2011, "services", :D, :U)
p3 = plot_fig_10(t_fig_9_10, 2011, "all industries", :D, :U)
title = "Figure 9: Correlation between U and D: Goods vs Services"
p_fig10 = plot(p1, p2, p3, layout=l, plot_title=title, plot_titlefontsize=10)
savefig(p_fig10, "images/figure10.png") # export image to folder


text:
using HOODESolver
using LinearAlgebra
using DoubleFloats
using Random
using SparseArrays
using Test
using DiffEqBase

function test_operator(T::DataType)
#    mat = [0 1;-1 0]
    for mat in [
    T.([0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]),
    T.([0 1;-1 0]),
    T.([0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 0;
     0 0 0 0 1 0; 0 0 0 -1 0 0; 0 0 0 0 0 0])
]
        epsilon = T(1//1024)
        linop = LinearHOODEOperator((1/epsilon)*mat)
        @test isconstant(linop)
        @test isapprox(linop.epsilon, epsilon)
        @test isapprox(linop.A, mat)
        L = T(1//16) * linop
        @test isconstant(L)
        @test isapprox(L.epsilon, epsilon*16)
        @test isapprox(L.A, mat)
        linop2 = LinearHOODEOperator(linop)
        @test isconstant(linop2)
        @test isapprox(linop2.epsilon, epsilon)
        @test isapprox(linop2.A, mat)
        L = T(1//16) * linop2
        @test isconstant(L)
        @test isapprox(L.epsilon, epsilon*16)
        @test isapprox(L.A, mat)
        linop3 = LinearHOODEOperator(epsilon, mat)
        @test isconstant(linop3)
        @test isapprox(linop3.epsilon, epsilon)
        @test isapprox(linop3.A, mat)
        L = T(1//16) * linop3
        @test isconstant(L)
        @test isapprox(L.epsilon, epsilon*16)
        @test isapprox(L.A, mat)
    end
end

@testset "test operator" begin
    test_operator(Float64)
    test_operator(Double64)
    test_operator(BigFloat)
end

function fct4bis( du, u, p, t)
    B=[0.1 1.1 0.2 05 ; -0.9 -0.12 -0.7 0.4 ; 0.5 0.66 0.7 0.8 ; -0.34 0.8 0 0.3]
    du[:] = B*u
    missing
end

function testcommon_interface_fct(linop)
    linopok = LinearHOODEOperator(linop)
    B=[0.1 1.1 0.2 05 ; -0.9 -0.12 -0.7 0.4 ; 0.5 0.66 0.7 0.8 ; -0.34 0.8 0 0.3]
    u0 = 2rand(4)-ones(4)
    sol_ref = exp(1.0*(1/linopok.epsilon*linopok.A+B))*u0
    fct = (u,p,t) -> B*u
    prob = SplitODEProblem(linop, fct, u0, (0.0, 1.0))
    sol = solve(prob, HOODEAB(), nb_t=100)
    @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
    fct = (u,p) -> B*u
    prob = SplitODEProblem(linop, fct, u0, (0.0, 1.0))
    sol = solve(prob, HOODEAB(), dt=0.01)
    @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
    fct = (u) -> B*u
    prob = SplitODEProblem(linop, fct, u0, (0.0, 1.0))
    sol = solve(prob, HOODEAB())
    @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)
    # prob = SplitODEProblem(linop, fct4bis, u0, (0.0, 1.0), missing, A, epsilon)
    # sol = solve(prob)
    # @test isapprox(sol_ref, sol[end], atol=1e-7, rtol=1e-6)

end
function testcommon_interface_epsilon()
    @time @testset "test interface while epsilon varying" begin
        Random.seed!(199881)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        fct = (u,p,t) -> B*u
        order = 5
        for i=1:1
            epsilon = big"1.0"/big"10"^(i*3)
            eps_v = convert(Float32, epsilon)
            t_max = big"1.0"
            linop = LinearHOODEOperator(epsilon, A)
            sol_ref = exp(t_max*(1/epsilon*A+B))*u0
            prob = SplitODEProblem(linop, fct, u0, (big"0.0", t_max))
            println("epsilon=$eps_v sol_ref=$sol_ref")
            nb = 100
            res_err = zeros(BigFloat,5)
            ind = 1
            while nb <= 1000
                sol = solve(prob, HOODEAB(5), nb_t=nb, dense=false)
                res_err[ind] = norm(sol[end]-sol_ref,Inf)
 #               println("\nnb=$nb epsilon=$eps_v err=$(res_err[ind])")
                nb *= 10
                ind += 1
            end
            coef = Float64(log(res_err[1]) - log(res_err[2]))/log(10)
            @test isapprox( order, coef, atol=0.01)
            println("epsilon=$eps_v coef log = $coef")
        end
    end
end
function testcommon_interface_interpolate()
    @time @testset "test interface for interpolation" begin
        seed=123873
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        order = 6
        nb = 100
        t_max = big"1.0"
        for i=1:1
            seed = seed*1341+6576123
            Random.seed!(seed)
            u0=rand(BigFloat,4)
            println("u0=$u0")
            B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
            fct = (u,p,t) -> B*u
            epsilon = big"0.4"/2.0^i
            eps_v = convert(Float32, epsilon)
            t_max = big"1.0"
            linop = LinearHOODEOperator(epsilon, A)
            prob = SplitODEProblem(linop, fct, u0, (big"0.0",t_max))
   #         sol = solve(prob, getprecision=false)
            sol = solve(prob,HOODEAB())
            m = 1/epsilon*A+B
            reftol=norm(exp(t_max*m)*u0-sol[end], Inf)*10
            for j=1:10
                t=rand(BigFloat)
                res_ex=exp(t*m)*u0
                res_ap=sol(t)
                println("t=$t norm=$(norm(res_ex-res_ap, Inf))")
                @test isapprox(res_ex, res_ap, atol=reftol, rtol=reftol*10)
            end
        end
    end
end
function testcommon_interface_interpolate_float()
    @time @testset "test interface for interpolation" begin
        seed=124363
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        order = 6
        nb = 100
        t_max = 1.0
        for i=1:1
            seed = seed*1341+6576127
            Random.seed!(seed)
            u0=rand(4)
            println("u0=$u0")
            B = 2rand(4,4)-ones(4,4)
            fct = (u,p,t) -> B*u
            epsilon = 0.4/2.0^i
            eps_v = convert(Float32, epsilon)
            t_max = 1.0
            linop = LinearHOODEOperator(epsilon, A)
            prob = SplitODEProblem(linop, fct, u0, (0.0,t_max))
   #         sol = solve(prob, getprecision=false)
            sol = solve(prob, HOODEAB())
            m = 1/epsilon*A+B
            reftol=norm(exp(t_max*m)*u0-sol[end], Inf)*10
            @test reftol < 1e-6
            for j=1:10
                t=rand()
                res_ex=exp(t*m)*u0
                res_ap=sol(t)
#                println("type ex=$(typeof(res_ex)) type ap=$(typeof(res_ap))")
#                println("t=$t norm=$(norm(res_ex-res_ap, Inf))")
                @test isapprox(res_ex, res_ap, atol=reftol, rtol=reftol*10)
            end
        end
    end
end
function testcommon_interface_short()
    seed=981881
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 7
    nb = 5
     for i=1:10
        seed += 10
        Random.seed!(seed)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        fct = (u,p,t) -> B*u
        epsilon = big"0.00004"/2.0^i
        eps_v = convert(Float32, epsilon)
        t_max = big"0.01"
        linop = LinearHOODEOperator(epsilon, A)
        prob = SplitODEProblem(linop, fct, u0, (big"0.0",t_max))
        sol = solve(prob, HOODEAB(order), getprecision=false, nb_t=nb)
        resnorm=norm(exp(t_max*(1/epsilon*A+B))*u0-sol[end], Inf)
        println("sol=$(sol[end])")
        println("solexact=$(exp(t_max*(1/epsilon*A+B))*u0)")
        println("resnorm=$resnorm")
        @test resnorm < 1e-10
    end
end
function tts_time(t_begin, t_end)
    @time @testset "test interface time from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = rand(BigFloat, 4)
        fct = (u,p,t) -> B*u
        epsilon=big"0.000000345"
        linop = LinearHOODEOperator(epsilon, A)
        prob = SplitODEProblem(linop, fct, u0, (t_begin,t_end))
        sol = solve(prob, HOODEAB(4), getprecision=false, nb_t=100)
        m = 1/epsilon*A+B
        solref = exp((t_end-t_begin)*m)*u0
        println("sol=$(sol[end]) solref=$solref norm=$(norm(sol[end]-solref,Inf))")
        @test isapprox(sol[end], solref,atol=1e-8, rtol=1e-7)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=exp((t-t_begin)*m)*u0
            res_ap=sol(t)
            @test isapprox(res_ex, res_ap, atol=1e-8, rtol=1e-7)
        end
    end
end
function tts_time_time(t_begin, t_end)
    @time @testset "test interface time time from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = 2rand(BigFloat, 4)-ones(BigFloat,4)
        fct = (u,p,t) -> B*u +t*p[1] +p[2]
        tuple_p = (2rand(BigFloat,4)-ones(BigFloat,4),
    2rand(BigFloat,4)-ones(BigFloat,4))
        epsilon=big"0.000003467"
        nb = 100
        order = 4
        linop = LinearHOODEOperator(epsilon, A)
        prob = SplitODEProblem(linop, fct, u0, (t_begin,t_end), tuple_p; B=B)
        sol = solve(prob, HOODEAB(4), getprecision=false, nb_t=100)
#        solref = getexactsol(sol.par_u0.parphi, u0, t_end)
        solref = getexactsol(sol.destats.par_u0.parphi, u0, t_end)
        println("sol=$(sol[end]) solref=$solref norm=$(norm(sol[end]-solref,Inf))")
        @test isapprox(sol[end], solref,atol=1e-8, rtol=1e-7)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=getexactsol(sol.destats.par_u0.parphi, u0, t)
            res_ap=sol(t)
            @test isapprox(res_ex, res_ap, atol=1e-8, rtol=1e-7)
        end
    end
end
function testcommon_interface_time()
    seed=7887
    Random.seed!(seed)
    tts_time(big"0.0",big"1.0")
    tts_time(big"0.0",big"1.4565656565")
    tts_time(big"2.234566",big"3.45766790123")
    tts_time(-big"0.8457676",big"0.56716")
    tts_time(-big"1.8111457676",-big"0.345456716")
    tts_time(big"1.8457676",big"0.56716")
    tts_time(-big"0.8457676",-big"1.56716")
    tts_time(-big"0.845722676",-big"0.56716")
end
function testcommon_interface_time_time()
    seed=788227
    Random.seed!(seed)
    tts_time_time(big"0.0",big"1.0")
    tts_time_time(big"0.0",big"1.4565656565")
    tts_time_time(big"11.234566",big"12.45766790123")
    tts_time_time(-big"0.8457676",big"0.56716")
    tts_time_time(-big"1.8111457676",-big"0.345456716")
    tts_time_time(big"1.8457676",big"0.56716")
    tts_time_time(-big"0.8457676",-big"1.56716")
    tts_time_time(-big"0.845722676",-big"0.56716")
end
@time @testset "test interface type of function" begin
    A =  [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    epsilon = 0.0001
    linop = LinearHOODEOperator(epsilon, A)
    testcommon_interface_fct(linop)
    linop2 = LinearHOODEOperator(epsilon, sparse(A))
    testcommon_interface_fct(linop2)
    testcommon_interface_fct(DiffEqArrayOperator((1/epsilon)*A))
end

testcommon_interface_interpolate_float()
testcommon_interface_time()
testcommon_interface_time_time()
testcommon_interface_epsilon()
testcommon_interface_interpolate()
testcommon_interface_short()


text:
#####################################################
### Stability parameters and stability correction ### ---------------------------------------------------
#####################################################

#' Monin-Obukhov Length
#'
#' calculates the Monin-Obukhov length.
#'
#' - data      DataFrame or matrix containing all required variables
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - ustar     Friction velocity (m s-1)
#' - H         Sensible heat flux (W m-2)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  k - von Karman constant (-)
#'                  g - gravitational acceleration (m s-2)
#'
#' # Details
 The Monin-Obukhov length (L) is given by:
#'
#'              ``L = - (\\rho * cp * ustar^3 * Tair) / (k * g * H)``
#'
#'              where ``rho`` is air density (kg m-3).
#'
#' # Value
 - L -: Monin-Obukhov length (m)
#'
#' # Note
#' Note that L gets very small for very low ustar values with implications
#'       for subsequent functions using L as input. It is recommended to filter
#'       data and exclude low ustar values (ustar < ~0.2) beforehand.
#'
#' #References
#' Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
#'
#' #See also
#' [`stability_parameter`](@ref)
#'
#' ```@example; output = false
#' ```
#' Monin_Obukhov_length(Tair=25,pressure=100,ustar=seq(0.2,1,0.1),H=seq(40,200,20))
#'
"""
"""
function Monin_Obukhov_length(data,Tair="Tair",pressure="pressure",ustar="ustar",
                                 H="H",constants=BigleafConstants())

  check_input(data,list(Tair,pressure,ustar,H))

  rho  = air_density(Tair,pressure,constants)
  Tair = Tair + constants.Kelvin
  MOL  = (-rho*constants.cp]*ustar^3*Tair) / (constants[:k]*constants[:g*H)

  return(MOL)
end

#' Stability Parameter "zeta"
#'
#' calculates "zeta", a parameter characterizing stratification in
#'              the lower atmosphere.
#'
#' - data      DataFrame or matrix containing all required variables
#' - Tair      Air temperature (degC)
#' - pressure  Atmospheric pressure (kPa)
#' - ustar     Friction velocity (m s-1)
#' - H         Sensible heat flux (W m-2)
#' - zr        Instrument (reference) height (m)
#' - d         Zero-plane displacement height (m)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  k - von Karman constant (-)
#'                  g - gravitational acceleration (m s-2)
#'
#' # Details
 The stability parameter ``\\zeta`` is given by:
#'
#'            ``\\zeta = (zr - d) / L``
#'
#'          where L is the Monin-Obukhov length (m), calculated from the function
#'          [`Monin_Obukhov_length`](@ref). The displacement height d can
#'          be estimated from the function [`roughness_parameters`](@ref).
#'
#' # Value
 - ``\\zeta`` - : stability parameter (-)
#'
#' ```@example; output = false
#' ```
#' df = DataFrame(Tair=25,pressure=100,ustar=seq(0.2,1,0.1),H=seq(40,200,20))
#' stability_parameter(df,zr=40,d=15)
#'
#' @export
function stability_parameter(data,Tair="Tair",pressure="pressure",ustar="ustar",
                                H="H",zr,d,constants=BigleafConstants())

  check_input(data,list(Tair,pressure,ustar,H))

  MOL  = Monin_Obukhov_length(data,Tair,pressure,ustar,H,constants)
  zeta = (zr - d) / MOL

  return(zeta)

end

#' Integrated Stability Correction Functions for Heat and Momentum
#'
#' dimensionless stability functions needed to correct deviations
#'              from the exponential wind profile under non-neutral conditions.
#'
#' - zeta         Stability parameter zeta (-)
#' - formulation  Formulation for the stability function. Either `Val(:Dyer_1970)`,
#'                     or `Val(:Businger_1971)`
#'
#' # Details
 The functions give the integrated form of the universal functions. They
#'          depend on the value of the stability parameter ``\\zeta``,
#'          which can be calculated from the function [`stability_parameter`](@ref).
#'          The integration of the universal functions is:
#'
#'            ``\\psi = -x * zeta``
#'
#'          for stable atmospheric conditions (``\\zeta`` >= 0), and
#'
#'            ``\\psi = 2 * log( (1 + y) / 2) ``
#'
#'          for unstable atmospheric conditions (``\\zeta`` < 0).
#'
#'          The different formulations differ in their value of x and y.
#'
#' # Value
 a DataFrame with the following columns:
#'          - psi_h: the value of the stability function for heat and water vapor (-)
#'          - psi_m: the value of the stability function for momentum (-)
#'
#' #References
#' Dyer, A_J., 1974: A review of flux-profile relationships.
#'             Boundary-Layer Meteorology 7, 363-372.
#'
#'             Dyer, A. J., Hicks, B_B., 1970: Flux-Gradient relationships in the
#'             constant flux layer. Quart. J. R. Meteorol. Soc. 96, 715-721.
#'
#'             Businger, J_A., Wyngaard, J. C., Izumi, I., Bradley, E. F., 1971:
#'             Flux-Profile relationships in the atmospheric surface layer.
#'             J. Atmospheric Sci. 28, 181-189.
#'
#'             Paulson, C_A., 1970: The mathematical representation of wind speed
#'             and temperature profiles in the unstable atmospheric surface layer.
#'             Journal of Applied Meteorology 9, 857-861.
#'
#'             Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
#'
#' ```@example; output = false
#' ```
#' zeta = seq(-2,0.5,0.05)
#' stability_correction(zeta)
#' stability_correction(zeta,formulation=Val(:Businger_1971))
#'
#' @export
function stability_correction(zeta,formulation=c(Val(:Dyer_1970),Val(:Businger_1971)))

  formulation  = match_arg(formulation)

  check_input(nothing,zeta)

  psi_h = psi_m = numeric()

  # universal functions
  if (formulation == Val(:Businger_1971))
    x_h = -7.8
    x_m = -6
    y_h = 0.95 * ( 1 - 11.6 * zeta)^0.5
    y_m = (1 - 19.3*zeta)^0.25
elseif (formulation == Val(:Dyer_1970))
    x_h = x_m = -5
    y_h       = (1 - 16 * zeta)^0.5
    y_m       = (1 - 16 * zeta)^0.25
end

  # integration of universal functions (after Paulson_1970 and Foken 2008)
  # stable
  stable = zeta >= 0 | ismissing(zeta)
  psi_h[stable] = x_h * zeta[stable]
  psi_m[stable] = x_m * zeta[stable]
  # unstable
  unstable = zeta < 0 | ismissing(zeta)
  psi_h[unstable] = 2 * log( (1 + y_h[unstable] ) / 2)
  psi_m[unstable] = 2 * log( (1 + y_m[unstable] ) / 2) +
                     log( ( 1 + y_m[unstable]^2 ) / 2)
                     -2 * atan(y_m[unstable]) + pi/2

  return(DataFrame(psi_h,psi_m))

end


text:
import StatsBase

# Each particle represents an independently evolving trajectory.
struct JumpParticle{uType}
    u::uType
    weight::Float64
end

function JumpParticle(setup)
    u = sample_initial_condition(setup.ensemble)
    w = initial_log_likelihood(setup.ensemble, u, setup.configuration.x_traj)
    JumpParticle(u, w)
end

function JumpParticle(parent::JumpParticle, setup)
    JumpParticle(copy(parent.u), 0.0)
end

function propagate(p::JumpParticle, tspan::Tuple{T,T}, setup) where T
    u_end, weight = propagate(setup.configuration, setup.ensemble, p.u, tspan)
    JumpParticle(u_end, weight)
end

# While this implementation of a particle is slower than the one above
# it keeps track of the genealogy of a particle (i.e. its parent, grandparent, etc.).
# This is useful for debugging and for collecting statistics.
mutable struct JumpParticleSlow{uType}
    u::uType
    weight::Float64
    parent::Union{JumpParticleSlow{uType},Nothing}

    function JumpParticleSlow(setup)
        u = sample_initial_condition(setup.ensemble)
        w = initial_log_likelihood(setup.ensemble, u, setup.configuration.x_traj)
        new{typeof(u)}(u, w, nothing)
    end

    function JumpParticleSlow(parent, setup)
        new{typeof(parent.u)}(copy(parent.u), 0.0, parent)
    end
end

struct Setup{Configuration,Ensemble}
    configuration::Configuration
    ensemble::Ensemble
end

function propagate(p::JumpParticleSlow, tspan::Tuple{T,T}, setup) where T
    u_end, weight = propagate(setup.configuration, setup.ensemble, p.u, tspan)
    p.u = u_end
    p.weight = weight
    p
end

weight(p::JumpParticle) = p.weight
weight(p::JumpParticleSlow) = p.weight

# This is the main routine to compute the marginal probability in SMC-PWS.
function sample(nparticles, dtimes, setup; inspect=Base.identity, new_particle=JumpParticle)
    particle_bag = [new_particle(setup) for i = 1:nparticles]
    weights = zeros(nparticles)
    log_marginal_estimate = zeros(length(dtimes))

    # First, handle initial condition
    for j in eachindex(particle_bag)
        weights[j] += weight(particle_bag[j])
    end
    log_marginal_estimate[1] = logmeanexp(weights)

    # At each time we do the following steps
    # 1) propagate all particles forwards in time and compute their weights
    # 2) update the log_marginal_estimate
    # 3) check whether we should resample the particle bag, and resample if needed
    for (i, tspan) in enumerate(zip(dtimes[begin:end - 1], dtimes[begin + 1:end]))

        # PROPAGATE
        for j in eachindex(particle_bag)
            particle_bag[j] = propagate(particle_bag[j], tspan, setup)
            weights[j] += weight(particle_bag[j])
        end

        # UPDATE ESTIMATE
        log_marginal_estimate[i+1] += logmeanexp(weights)

        if (i + 1) == lastindex(dtimes)
            break
        end

        # RESAMPLE IF NEEDED
        prob_weights = StatsBase.weights(exp.(weights .- maximum(weights)))

        # We only resample if the effective sample size becomes smaller than 1/2 the number of particles
        effective_sample_size = 1/sum(p -> (p / sum(prob_weights)) ^ 2, prob_weights)
        if effective_sample_size < nparticles / 2
            # sample parent indices
            parent_indices = systematic_sample(prob_weights)

            particle_bag = map(parent_indices) do k
                new_particle(particle_bag[k], setup)
            end
            weights .= 0.0
            log_marginal_estimate[i+1:end] .= log_marginal_estimate[i+1]
        end
    end

    inspect(particle_bag) #< useful for collecting statistics
    log_marginal_estimate
end

"""
    systematic_sample(weights)::Vector{Int}

Take random samples from `1:len(weights)` using weights given by the corresponding values
in the `weights` array.

The number of samples returned is equal to the length of `weights`.
"""
function systematic_sample(weights)
    N = length(weights)
    inc = 1/N
    x = inc * rand()
    j=1
    y = weights[j] / sum(weights)
    result = zeros(Int, N)
    for i=1:N
        while y < x
            j += 1
            y += weights[j] / sum(weights)
        end
        result[i] = j
        x += inc
    end
    result
end

struct SMCEstimate
    num_particles::Int
end

name(x::SMCEstimate) = "SMC"

struct SMCResult <: SimulationResult
    log_marginal_estimate::Vector{Float64}
end

log_marginal(result::SMCResult) = result.log_marginal_estimate

function simulate(algorithm::SMCEstimate, initial, system; kwargs...)
    setup = Setup(initial, system)
    log_marginal_estimate = sample(algorithm.num_particles, system.dtimes, setup; kwargs...)
    SMCResult(log_marginal_estimate)
end


text:
# A `Samples` object represents all of the logged information about the
# execution of a benchmark. Each object contains five vectors of equal length.
#
# Fields:
#
#     evaluations::Vector{Float64}: The number of times the expression was
#         evaluated for each sample. NB: We use a `Float64` vector to avoid
#         type conversion when calling the `linreg()` function, although we
#         always evaluate functions an integral number of times.
#
#     elapsed_times::Vector{Float64}: The execution time in nanoseconds for
#         each sample. NB: We use a `Float64` vector to avoid type conversion
#         when calling the `linreg()` function, although the timing functions
#         use `UInt` values to represent nanoseconds.
#
#     gc_times::Vector{Float64}: The GC time in nanoseconds for each sample.
#         NB: We use a `Float64` vector for similarity with `elapsed_times`,
#         although the timing functions use `UInt` values to represent
#         nanoseconds.
#
#     bytes_allocated::Vector{Int}: The total number of bytes allocated during
#         each sample.
#
#     allocations::Vector{Int}: The total number of allocation operations
#         performed during each sample.

immutable Samples
    evaluations::Vector{Float64}
    elapsed_times::Vector{Float64}
    gc_times::Vector{Float64}
    bytes_allocated::Vector{Int}
    allocations::Vector{Int}

    function Samples()
        new(
            Vector{Float64}(0),
            Vector{Float64}(0),
            Vector{Float64}(0),
            Vector{Int}(0),
            Vector{Int}(0),
        )
    end
end

# Push values for all metrics to a `Samples` object.
#
# Arguments:
#
#     s::Samples: The `Samples` object whose contents should be emptied.
#
#     evaluations::Real: The number of evaluations of the core expression for
#         this core expression.
#
#     elapsed_time::Real: The time in nanoseconds to evaluate the core
#         expression.
#
#     gc_time::Real: The time spent in the GC while evaluating the core
#         expression.
#
#     bytes_allocated::Real: The number of bytes allocated while evaluating the
#         core expression.
#
#     allocations::Real: The number of allocation operations while evaluating
#         the core expression.

function Base.push!(
    s::Samples,
    evaluations::Real,
    elapsed_time::Real,
    gc_time::Real,
    bytes_allocated::Real,
    allocations::Real,
)
    push!(s.evaluations, evaluations)
    push!(s.elapsed_times, elapsed_time)
    push!(s.gc_times, gc_time)
    push!(s.bytes_allocated, bytes_allocated)
    push!(s.allocations, allocations)
    return
end

# Empty all of the five vectors from a `Samples` object.
#
# Arguments:
#
#     s::Samples: The `Samples` object whose contents should be emptied.

function Base.empty!(s::Samples)
    empty!(s.evaluations)
    empty!(s.elapsed_times)
    empty!(s.gc_times)
    empty!(s.bytes_allocated)
    empty!(s.allocations)
    return
end

# Pretty-print information about the `Samples` object.
#
# Arguments:
#
#     io::IO: An `IO` object to be written to.
#
#     s::Samples: The `Samples` object that we want to print to `io`.

function Base.show(io::IO, s::Samples)
    names = Compat.String["Number of samples"]
    values = Any[length(s.elapsed_times)]

    @printf(io, "================== Benchmark Samples ======================\n")
    max_length = maximum([length(n) for n in names]) +  1
    for (n, v) in zip(names, values)
        @printf(io, "%s: %s\n", lpad(n, max_length), v)
    end
end

# Log information about a set of samples to a CSV file.
#
# Arguments:
#
#     filename::AbstractString: The name of a file to which we'll write information
#         the `Samples` object, `s`.
#
#     s::Samples: The set of samples that we want to log to disk.
#
#     append::Bool: Should we overwrite the file or append to existing content?
#         Defaults to false.

function Base.writecsv(filename::AbstractString, s::Samples, append::Bool = false)
    if append
        io = open(filename, "a")
    else
        io = open(filename, "w")
    end
    println(
        io,
        join(
            [
                "evaluations",
                "elapsed_time",
                "gc_time",
                "bytes_allocated",
                "allocations",
            ],
            ','
        )
    )
    for i in 1:length(s.evaluations)
        println(
            io,
            join(
                [
                    string(s.evaluations[i]),
                    string(s.elapsed_times[i]),
                    string(s.gc_times[i]),
                    string(s.bytes_allocated[i]),
                    string(s.allocations[i]),
                ],
                ','
            )
        )
    end
    close(io)
end


text:
## http://docs.sympy.org/0.7.2/modules/integrals/integrals.html

integrals_sympy_methods = (:integrate, :Integral,
                           :mellin_transform,
                           :inverse_mellin_transform,
                           :laplace_transform,
                           :inverse_laplace_transform,
                           :fourier_transform,
                           :inverse_fourier_transform,
                           :sine_transform,
                           :inverse_sine_transform,
                           :cosine_transform,
                           :inverse_cosine_transform,
                           :hankel_transform,
                           :inverse_hankel_transform,
                           :line_integrate,
                           :deltaintegrate,
                           :ratint,
                           :heurisch,
                           :trigintegrate,
                           :line_integrate
                           )

integrals_instance_methods = (:as_sum,
                              :limits,
                              :transform,
                              :variables)

"""
Create a parameterized curve for line integrals

```
@vars t x y
C = Curve([exp(t)+1, exp(t)-1], (t, 0, log(Sym(2))))
line_integrate(x + y, C, [x,y])
```

"""
Curve{T<:Sym}(exs::Vector{T}, p) = sympy_meth(:Curve, exs, p)
export Curve

"""
Dirac delta for integration

[SymPy Documentation](http://docs.sympy.org/dev/modules/functions/special.html)
"""
DiracDelta(x::Number) = convert(Function, sympy[:DiracDelta])(x)
export DiracDelta

"""
Heaviside function for integration.

`H(x) = x < 0 ? 0 : (x > 0 ? 1 : 1/2)`

[SymPy Documentation](http://docs.sympy.org/dev/modules/functions/special.html)
"""
Heaviside(x::Number) = convert(Function, sympy[:Heaviside])(x)
export Heaviside

## Alternate interface for simple integral
"""

The `integrate` function has its limits specified with tuples of the type `(var, a, b)`. This
profides a simpler interface for one-dimensional integrals: `integrate(ex, var, a, b)`

"""
integrate(s::Sym, x::Sym, from::SymOrReal, to::SymOrReal) = integrate(s, (x, convert(Sym,from), convert(Sym,to)))

"Symbolically integrate a function"
function integrate(f::Function)
    x = Sym("x")
    integrate(f(x), x)
end

"Symbolically integrate a function over `[a,b]`"
function integrate(f::Function, from::SymOrReal, to::SymOrReal)
    x = Sym("x")
    integrate(f(x), x, from, to)
end

## summation
summations_sympy_methods = (
                            :summation, :Sum
                            )

summations_instance_methods  = (
                                :euler_maclaurin,
                                )
summations_object_properties  = (:is_zero, :is_number)


text:
"""
    $(SIGNATURES)

One dimesional trapezoidal integration over a regularly discretized axis.
"""
function trapz(x::StepRange, y::AbstractArray)
    @assert length(x) == length(y)
    trapz(x.step, y)
end
function trapz(step::Number, y::AbstractArray)
    step * (sum(y) - (y[1] + y[end]) / 2)
end
function integral!(f::Function, axis::AbstractArray,
                   y::AbstractArray, result::AbstractArray)
    @argcheck length(axis) == length(y)
    @argcheck length(y) == length(result)
    y0 = f(axis[1], y[1])
    @argcheck dimension(axis[1] * y0) == dimension(eltype(result))
    result[1] = 0 * unit(eltype(result))
    for i in 2:length(result)
        y1 = f(axis[i], y[i])
        result[i] = result[i - 1] + (axis[i] - axis[i - 1]) * (y1 + y0) / 2
        y0 = y1
    end
    result
end
function integral(f::Function, axis::AbstractArray, y::AbstractArray)
    x0, y0 = zero(eltype(axis)), zero(eltype(y))
    f0 = f(x0, y0)
    if dimension(f0) == dimension(x0)^-1
        U = typeof(ustrip(x0 * f0))
    else
        U = typeof(x0 * f0)
    end
    result = similar(y, U)
    integral!(f, axis, y, result)
end

function add_radial_hartree!(ρ::AbstractArray, potential::AbstractArray,
                             axis::AbstractArray)
    @argcheck axis[1] ≈ zero(eltype(axis))
    @argcheck length(ρ) == length(potential)
    @argcheck length(potential) == length(axis)
    @argcheck all(axis[2:end] .> (0 * unit(eltype(axis))))
    const N = length(axis)

    # e²/(4πϵ₀)4π/r \int_0^r r'^2 * ρ(r') dr'
    ∫r²ρdr = integral((x, y) -> x*x*y, axis, ρ)
    potential[2:end] .+= (view(∫r²ρdr, 2:N) ./ view(axis, 2:N))u"e₀^2*ϵ₀^-1"

    # e²/(4πϵ₀)4π \int_0^r r' * ρ(r') dr'
    ∫rρdr = integral((x, y) -> x*y, axis, ρ)
    potential[1:end] .+= (∫rρdr[end] - ∫rρdr)u"e₀^2*ϵ₀^-1"

    potential
end
@inline function radial_hartree!{T <: DFTUnits.Ρ, Q <: DFTUnits.Ε}(ρ::AxisArray{T},
                                                                   potential::AxisArray{Q})
    radial_hartree!(SimpleTraits.trait(AtomicDFT.HasSpinDim{typeof(ρ)}), ρ, potential)
end
@traitfn function radial_hartree!(ρ::::(!HasSpinDim), potential)
    @argcheck :radius ∈ axisnames(ρ)
    const i = findfirst(axisnames(ρ), :radius)
    add_radial_hartree!(ρ, potential, axisvalues(ρ)[i])
end
@traitfn function radial_hartree!(ρ::::(HasSpinDim), potential)
    @argcheck :radius ∈ axisnames(ρ)
    const i = findfirst(axisnames(ρ), :radius)
    add_radial_hartree!(view(ρ, Axis{:spin}(:↑)) .+ view(ρ, Axis{:spin}(:↓)),
                        potential, axisvalues(ρ)[i])
end
@lintpragma("Ignore use of undeclared variable similar_potential")
radial_hartree{Q <: DFTUnits.Ρ}(ρ::AxisArray{Q}) = radial_hartree!(ρ, similar_potential(ρ))

# immutable LDAPotential{T <: LibXC.CReal} <: AbstractPotential
#     functional::XCFunctional{T}
#     work::Array{T}
# end
#
#
# function RadialKohnSham(xc::Vararg{Symbol}; charge::Unitful.ChargeUnit=1u"e₀",
#                         polarized::Bool=false)
#     LDARadialKohnSham(map(x -> XCFunctional(x, polarized), xc), charge, ρ)
# end
#
# function add_potential!{Ρ <: DFTUnits.Ρ, 𝐕 <: DFTUnits.Ε}(func::XCFunctional,
#                                                           ρ::AxisArray{Q},
#                                                           potential::AxisArray{𝐕})
#    potential
# end
#
# function potential!{Ρ <: DFTUnits.Ρ, 𝐕 <: DFTUnits.Ε}(func::LDARadialKohnSham,
#                                                       ρ::AxisArray{Q},
#                                                       potential::AxisArray{𝐕})
#     @argcheck size(ρ) == size(potential)
#     fill!(potential, zero(eltype(potential)))
#
#     𝐯 = potential(func.ϵxcs[1], ρ)
#     potential += 𝐯
#
#     for i in eachindex(drop(func.ϵxcs, 1))
#         add_potential!(func.ϵxcs[1], ρ, 𝐯)
#     end
#
#     radial_hartree!(ρ, ∂ϵ_∂ρ)
# end


text:
using Yao, Yao.EasyBuild, Test

@vars α β γ
c = chain(put(3, 2 => Rx(α)), control(3, 2, 1 => Ry(β)), put(3, (1, 2) => rot(kron(X, X), γ)))
h = heisenberg(3);
es = expect(h, zero_state(Basic, 3) => c)
gs = expect'(h, zero_state(Basic, 3) => c).second
assign = Dict(α => 0.5, β => 0.7, γ => 0.8)

cn = subs(Float64, c, assign...)
en = expect(h, zero_state(3) => cn)
gn = expect'(h, zero_state(3) => cn).second
es_ = subs(es, assign...)
gs_ = map(x -> subs(x, assign...), gs)
@test ComplexF64.(gs_) ≈ gn
@test ComplexF64(es_) ≈ en


text:
using Flux

cd(@__DIR__)

include("pointcnn.jl")

struct Classifier
    cls::Int
    pcnn1::PointCNN
    pcnn2::PointCNN
    pcnn3::PointCNN
    pcnn4::PointCNN
    fc1::Dense
    fc2::Dense
end

Classifier(cls::Int) = Classifier(
    cls,
    PointCNN(1024, 0, 32, 16, 8),
    PointCNN(384, 32, 64, 8, 12;dilation=2),
    PointCNN(128, 64, 128, 16, 16;dilation=2),
    PointCNN(128, 128, 256, 32, 16;dilation=3),
    Dense(256, 128),
    Dense(128, cls)
    )

function (m::Classifier)(pts)
    # pts : (n_p, 1, 3, batch)
    p1, f1 = m.pcnn1(pts, nothing)
    p2, f2 = m.pcnn2(p1, f1)
    p3, f3 = m.pcnn3(p2, f2)
    p4, f4 = m.pcnn4(p3, f3)
    f = reshape(f4, (size(f4)[1], size(f4)[3], size(f4)[end]))
    fp = permutedims(f, [2, 1, 3])
    logpred = reshape(sum(cat(
        [m.fc2(m.fc1(fp[:,:,i])) for i in 1:size(fp)[3]]..., dims=3
    ), dims=2), (m.cls, :))

    return logpred
end

Flux.@functor Classifier


text:
using Griddly
using Test
using Makie

image_path = joinpath(@__DIR__,"..","..","resources","images")
shader_path = joinpath(@__DIR__,"..","..","resources","shaders")
gdy_path = joinpath(@__DIR__,"..","..","resources","games")

@testset "Griddly" begin

	gdy_reader = Griddly.GDYReader(image_path,shader_path)
	grid = Griddly.load!(gdy_reader,joinpath(gdy_path,"Single-Player/Mini-Grid/minigrid-doggo.yaml"))
	# Set of Test to see if we load from file correctly and if basic functions works
	@testset "LoadFromFile" begin
    	@test Griddly.get_tile_size(grid) == 24
    	@test Griddly.get_all_available_actions(grid) == ["move"]
    	@test Griddly.get_player_available_actions(grid) == ["move"]
    	@test Griddly.get_non_player_available_actions(grid) == []
    	@test Griddly.get_input_ids(grid,"move") == [1,2,3]
    	@test Griddly.get_avatar_object(grid) =="doggo"
    	Griddly.set_tile_size!(grid,16)
    	@test Griddly.get_tile_size(grid) == 16
        Griddly.set_tile_size!(grid,24)
    end

    @testset "Game Entity" begin
        game = Griddly.create_game(grid,Griddly.SPRITE_2D)
        @test Griddly.get_player_count(grid) == 1
        @test Griddly.get_num_players(game) == 0
        player1 = Griddly.register_player!(game,"Bob", Griddly.BLOCK_2D)
        @test Griddly.get_num_players(game) == 1
        Griddly.init!(game)
        Griddly.load_level!(grid,1)
        Griddly.reset!(game)
        @test Griddly.get_width(grid) == 7
        @test Griddly.get_height(grid) == 7

		# test a random game with different level
		function run()
			for l in 0:4
				Griddly.load_level!(grid,l)
				Griddly.reset!(game)
				for j in 1:1000
					dir = rand(0:5)
					reward, done = Griddly.step_player!(player1,"move", [dir])
				end
			end
		end
		println(@timed run())

		@testset "RenderTools" begin

	        Griddly.load_level!(grid,1)
	        Griddly.reset!(game)
	        # test if the render tool is working
			function run_window()
		        for l in 0:4
					render_window = RenderWindow(700,700)
					display(render_window.scene)
		            Griddly.load_level!(grid,l)
		            Griddly.reset!(game)
		            observation = Griddly.observe(game)
		            observation = Griddly.get_data(observation)
		            for j in 1:800
		    	        dir = rand(0:5)

		    	        reward, done = Griddly.step_player!(player1,"move", [dir])

		                observation = Griddly.observe(game)
		                observation = Griddly.get_data(observation)
		                render(render_window,observation)
		    	    end
		        end
			end

			function run_window_player()
		        for l in 0:4
					render_window_player = RenderWindow(700,700)
					display(render_window.scene)
		            Griddly.load_level!(grid,l)
		            Griddly.reset!(game)
		            for j in 1:1000
		                dir = rand(0:5)

		                reward, done = Griddly.step_player!(player1,"move", [dir])

		                player1_tiles = Griddly.observe(player1)
		                player1_tiles = Griddly.get_data(player1_tiles)
		                render(render_window_player, player1_tiles)
		            end
		        end
			end

	        # test if we can capture video
			function run_video()
				video = VideoRecorder((700,700),"test_video2";saving_path="julia/Griddly/test/")
				display(video.scene)
		        io = start_video(video)
		        Griddly.load_level!(grid,4)
		        Griddly.reset!(game)
		        observation = Griddly.observe(game)
		        observation = Griddly.get_data(observation)
		        for j in 1:1000
		            dir = rand(0:5)

		            reward, done = Griddly.step_player!(player1,"move", [dir])

		            observation = Griddly.observe(game)
		            observation = Griddly.get_data(observation)
		            add_frame!(video,io,observation)
		        end
		        save_video(video,io)
			end

			# test if we can render several observation
			function run_multiple_screen()
				Griddly.load_level!(grid,3)
				Griddly.reset!(game)
				screen = MultipleScreen(1000,1000;nb_scene=5)
				display(screen.scene)
				for j in 1:200
					dir = rand(0:5)

					reward, done = Griddly.step_player!(player1,"move", [dir])

					observation = Griddly.observe(game)
					observation = Griddly.get_data(observation)
					player1_tiles = Griddly.observe(player1)
					player1_tiles = Griddly.get_data(player1_tiles)
					render_multiple(screen, [player1_tiles,player1_tiles,player1_tiles,observation,observation])
				end
			end

			# test if we can save a frame
			function run_save_frame()
				Griddly.load_level!(grid,3)
				Griddly.reset!(game)
				observation = Griddly.observe(game)
				observation = Griddly.get_data(observation)
				for j in 1:10
					dir = rand(0:5)

					reward, done = Griddly.step_player!(player1,"move", [dir])

					observation = Griddly.observe(game)
					observation = Griddly.get_data(observation)
				end
				Griddly.save_frame(observation,(700,700),"test_capture_frame2";file_path="julia/Griddly/test/",format=".jpg")
			end

	        # println(@timed run_multiple_screen())
		end
    end
	@testset "LoadFromString & Multiplayer" begin
	    gdy_string = """Version: "0.1"
	Environment:
	  Name: GriddlyRTS
	  TileSize: 16
	  BackgroundTile: oryx/oryx_tiny_galaxy/tg_sliced/tg_world/tg_world_floor_panel_metal_a.png
	  Player:
	    Count: 2
	  Termination:
	    Lose:
	      - eq: [base:count, 0] # If the player has no bases
	    Win:
	      - eq: [_score, 10] # First player to 10 reward points
	  Levels:
	    - |
	      W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W
	      W  .  .  B1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  H1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  H1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  M  M  .  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  M  .  M  M  M  .  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  M  M  M  M  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  M  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  M  .  .  M  .  M  M  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  M  M  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  M  M  M  M  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  H2 .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  H2 .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  B2 .  .  W
	      W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W
	    - |
	      W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W
	      W  .  .  B1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  H1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  H1 .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  M  M  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  M  M  .  M  M  M  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  M  .  M  M  M  .  M  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W  W  W  W  W  w  W
	      W  W  W  W  W  W  W  W  W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  .  .  W  W  .  .  .  w  w  w  w  w  w  .  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  w  w  .  .  .  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  M  M  M  M  .  .  .  .  W  W  W  W  W  w  W
	      W  .  .  .  .  P1 w  .  .  .  .  .  .  .  .  M  M  M  M  M  M  .  W  W  W  W  W  W  W  W
	      W  .  .  .  .  P1 w  .  .  .  .  .  .  .  M  .  .  M  .  M  M  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  P2 w  .  .  .  .  .  .  .  M  M  M  M  M  M  M  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  P2 w  .  .  .  .  .  .  .  .  .  M  M  M  M  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  P1 w  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  H2 .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  H2 .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  W
	      W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  B2 .  .  W
	      W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W  W

	Actions:
	  - Name: gather
	    Behaviours:
	      - Src:
	          Object: harvester
	          Commands:
	            - incr: resources
	            - reward: 1
	        Dst:
	          Object: minerals
	          Commands:
	            - decr: resources
	      - Src:
	          Object: harvester
	          Commands:
	            - decr: resources
	            - reward: 1
	        Dst:
	          Object: base
	          Commands:
	            - incr: resources

	  - Name: move
	    Behaviours:
	      - Src:
	          Object: [harvester, puncher, pusher, movable_wall]
	          Commands:
	            - mov: _dest # mov will move the object, _dest is the destination location of the action
	        Dst:
	          Object: _empty

	      - Src:
	          Object: pusher
	          Commands:
	            - mov: _dest # mov will move the object, _dest is the destination location of the action
	        Dst:
	          Object: [movable_wall, harvester, puncher]
	          Commands:
	            - cascade: _dest # reapply the same action to the dest location of the action

	  - Name: punch
	    Behaviours:
	      - Src:
	          Object: puncher
	          Commands:
	            - reward: 1
	        Dst:
	          Object: [puncher, harvester, pusher, base]
	          Commands:
	            - decr: health
	            - eq:
	                Arguments: [0, health]
	                Commands:
	                  - remove: true

	Objects:
	  - Name: minerals
	    MapCharacter: M
	    Variables:
	      - Name: resources
	        InitialValue: 10
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_items/tg_items_crystal_green.png
	      Block2D:
	        - Shape: triangle
	          Color: [0.0, 1.0, 0.0]
	          Scale: 1.0

	  - Name: harvester
	    MapCharacter: H
	    Variables:
	      - Name: resources
	        InitialValue: 0
	      - Name: health
	        InitialValue: 10
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_monsters/tg_monsters_jelly_d1.png
	      Block2D:
	        - Shape: square
	          Color: [0.6, 0.2, 0.2]
	          Scale: 0.5

	  - Name: pusher
	    MapCharacter: P
	    Variables:
	      - Name: health
	        InitialValue: 10
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_monsters/tg_monsters_crawler_queen_d1.png
	      Block2D:
	        - Shape: square
	          Color: [0.2, 0.2, 0.6]
	          Scale: 1.0

	  - Name: puncher
	    MapCharacter: p
	    Variables:
	      - Name: health
	        InitialValue: 5
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_monsters/tg_monsters_beast_d1.png
	      Block2D:
	        - Color: [0.2, 0.6, 0.6]
	          Shape: square
	          Scale: 0.8

	  - Name: fixed_wall
	    MapCharacter: W
	    Observers:
	      Sprite2D:
	        - TilingMode: WALL_2 # Will tile walls with two images
	          Image:
	            - oryx/oryx_tiny_galaxy/tg_sliced/tg_world_fixed/img33.png
	            - oryx/oryx_tiny_galaxy/tg_sliced/tg_world_fixed/img40.png
	      Block2D:
	        - Color: [0.5, 0.5, 0.5]
	          Shape: square

	  - Name: movable_wall
	    MapCharacter: w
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_world_fixed/img282.png
	      Block2D:
	        - Color: [0.8, 0.8, 0.8]
	          Shape: square

	  - Name: base
	    MapCharacter: B
	    Variables:
	      - Name: resources
	        InitialValue: 0
	      - Name: health
	        InitialValue: 10
	    Observers:
	      Sprite2D:
	        - Image: oryx/oryx_tiny_galaxy/tg_sliced/tg_world_fixed/img324.png
	      Block2D:
	        - Color: [0.8, 0.8, 0.3]
	          Shape: triangle"""

	    grid = Griddly.load_string!(gdy_reader,gdy_string)
	    @test Griddly.get_tile_size(grid) == 16
		@test Griddly.get_all_available_actions(grid) == ["gather","move","punch"]
    	@test Griddly.get_player_available_actions(grid) == ["gather","move","punch"]
    	@test Griddly.get_non_player_available_actions(grid) == []
    	@test Griddly.get_input_ids(grid,"move") == [1,2,3,4]
    	@test Griddly.get_input_ids(grid,"punch") == [1,2,3,4]
		@test Griddly.get_input_ids(grid,"gather") == [1,2,3,4]
		actions_map = Griddly.get_player_actions_dict(grid)
		@test actions_map == Dict("move"=>[1,2,3,4],
								  "punch"=>[1,2,3,4],
								  "gather"=>[1,2,3,4])
		game = Griddly.create_game(grid,Griddly.SPRITE_2D)
        @test Griddly.get_player_count(grid) == 2
        @test Griddly.get_num_players(game) == 0
        player1 = Griddly.register_player!(game,"Bob", Griddly.BLOCK_2D)
        @test Griddly.get_num_players(game) == 1
		player2 = Griddly.register_player!(game,"Alice",Griddly.BLOCK_2D)
        @test Griddly.get_num_players(game) == 2
        Griddly.init!(game)
        Griddly.load_level!(grid,1)
        Griddly.reset!(game)
        @test Griddly.get_width(grid) == 30
        @test Griddly.get_height(grid) == 30

		# test a random game with different level and random action
		function run()
			for l in 0:1
				Griddly.load_level!(grid,l)
				Griddly.reset!(game)
				for j in 1:1000
					x = rand(1:Griddly.get_width(grid))
					y = rand(1:Griddly.get_height(grid))

					# Alternate between player_1 and player_2 actions
					action = collect(keys(actions_map))[rand(1:length(actions_map))]
					dir = rand(1:length(actions_map[action]))
					println("action:$action,dir:$dir")
			        if j % 2 == 0
			            player_1_step_result = Griddly.step_player!(player1,action, [x, y, dir])
			        else
			            player_2_step_result = Griddly.step_player!(player2,action, [x, y, dir])
					end
				end
			end
		end
		# println(@timed run())

		@testset "RenderTool" begin
			function run_multiple_screen()
				Griddly.load_level!(grid,1)
				Griddly.reset!(game)
				screen = MultipleScreen(1000,1000;nb_scene=3)
				display(screen.scene)
				for j in 1:1000
					x = rand(1:Griddly.get_width(grid))
					y = rand(1:Griddly.get_height(grid))

					# Alternate between player_1 and player_2 actions
					action = collect(keys(actions_map))[rand(1:length(actions_map))]
					dir = rand(1:length(actions_map[action]))
					println("action:$action,dir:$dir")
			        if j % 2 == 0
			            player_1_step_result = Griddly.step_player!(player1,action, [x, y, dir])
			        else
			            player_2_step_result = Griddly.step_player!(player2,action, [x, y, dir])
					end
					player1_tiles = Griddly.get_data(Griddly.observe(player1))
					player2_tiles = Griddly.get_data(Griddly.observe(player2))
					observation = Griddly.get_data(Griddly.observe(game))
					render_multiple(screen,[player1_tiles, player2_tiles,observation])
				end
			end
			# println(@timed run_multiple_screen())
		end
	end
end


text:
using Distributions
using FillArrays
using StatsPlots

using LinearAlgebra
using Random

# Set a random seed.
Random.seed!(3)

# Define Gaussian mixture model.
w = [0.5, 0.5]
μ = [-3.5, 0.5]
mixturemodel = MixtureModel([MvNormal(Fill(μₖ, 2), I) for μₖ in μ], w)

# We draw the data points.
N = 60
x = rand(mixturemodel, N);

scatter(x[1, :], x[2, :]; legend=false, title="Synthetic Dataset")

using Turing

@model function gaussian_mixture_model(x)
    # Draw the parameters for each of the K=2 clusters from a standard normal distribution.
    K = 2
    μ ~ MvNormal(Zeros(K), I)

    # Draw the weights for the K clusters from a Dirichlet distribution with parameters αₖ = 1.
    w ~ Dirichlet(K, 1.0)
    # Alternatively, one could use a fixed set of weights.
    # w = fill(1/K, K)

    # Construct categorical distribution of assignments.
    distribution_assignments = Categorical(w)

    # Construct multivariate normal distributions of each cluster.
    D, N = size(x)
    distribution_clusters = [MvNormal(Fill(μₖ, D), I) for μₖ in μ]

    # Draw assignments for each datum and generate it from the multivariate normal distribution.
    k = Vector{Int}(undef, N)
    for i in 1:N
        k[i] ~ distribution_assignments
        x[:, i] ~ distribution_clusters[k[i]]
    end

    return k
end

model = gaussian_mixture_model(x);

sampler = Gibbs(PG(100, :k), HMC(0.05, 10, :μ, :w))
nsamples = 100
nchains = 3
chains = sample(model, sampler, MCMCThreads(), nsamples, nchains);

let
    # Verify that the output of the chain is as expected.
    for i in MCMCChains.chains(chains)
        # μ[1] and μ[2] can switch places, so we sort the values first.
        chain = Array(chains[:, ["μ[1]", "μ[2]"], i])
        μ_mean = vec(mean(chain; dims=1))
        @assert isapprox(sort(μ_mean), μ; rtol=0.1)
    end
end

plot(chains[["μ[1]", "μ[2]"]]; colordim=:parameter, legend=true)

plot(chains[["w[1]", "w[2]"]]; colordim=:parameter, legend=true)

chain = chains[:, :, 1];

# Model with mean of samples as parameters.
μ_mean = [mean(chain, "μ[$i]") for i in 1:2]
w_mean = [mean(chain, "w[$i]") for i in 1:2]
mixturemodel_mean = MixtureModel([MvNormal(Fill(μₖ, 2), I) for μₖ in μ_mean], w_mean)

contour(
    range(-7.5, 3; length=1_000),
    range(-6.5, 3; length=1_000),
    (x, y) -> logpdf(mixturemodel_mean, [x, y]);
    widen=false,
)
scatter!(x[1, :], x[2, :]; legend=false, title="Synthetic Dataset")

assignments = [mean(chain, "k[$i]") for i in 1:N]
scatter(
    x[1, :],
    x[2, :];
    legend=false,
    title="Assignments on Synthetic Dataset",
    zcolor=assignments,
)

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
# AbstractString interface for JSON.String

"""
Processes escape sequences and return a `Base.Substring`
(without copying where possible).
"""
Base.convert(::Type{SubString}, s::JSON.String) = parse_string(s.s, s.i)[1]

function parse_string(s, i)
    last_i, has_escape = scan_string(s, i)
    if !has_escape
        return SubString(s, next_i(s, i), prevind(s, last_i)), last_i
    else
        return unescape_json_string(s, i, last_i), last_i
    end
end

Base.convert(::Type{Base.String}, s::JSON.String) =
    convert(Base.String, convert(Base.SubString, s))

Base.String(s::JSON.String) = convert(Base.String, s)
Base.SubString(s::JSON.String) = convert(Base.SubString, s)

# Abstract String Interface

Base.IteratorSize(::Type{JSON.String{T}}) where T = Base.SizeUnknown()

Base.ncodeunits(s::JSON.String) = scan_string(s.s, s.i)[1] - s.i - 1

@propagate_inbounds(
Base.isvalid(s::JSON.String, i::Integer) = string_index_isvalid(s.s, s.i + i))

Base.codeunit(s::JSON.String) = codeunit(s.s)

@propagate_inbounds(
Base.codeunit(s::JSON.String, i::Integer) = codeunit(s.s, s.i + i))

@propagate_inbounds(
Base.iterate(s::JSON.String, i::Integer=1) = _iterate(s, i))
@propagate_inbounds(
function _iterate(s::JSON.String, i)
    i, c = json_char(s.s, s.i + i)
    if c == nothing
        return nothing
    end
    return c, i - s.i
end)

# Unescaping JSON Strings

"""
https://tools.ietf.org/html/rfc7159#section-7

    char = unescaped /
        escape (
            %x22 /          ; "    quotation mark  U+0022
            %x5C /          ; \\    reverse solidus U+005C
            %x2F /          ; /    solidus         U+002F
            %x62 /          ; b    backspace       U+0008
            %x66 /          ; f    form feed       U+000C
            %x6E /          ; n    line feed       U+000A
            %x72 /          ; r    carriage return U+000D
            %x74 /          ; t    tab             U+0009
            %x75 4HEXDIG )  ; uXXXX                U+XXXX
"""
unescape_c(c) = c == '"'  ? c :
               c == '\\' ? c :
               c == '/'  ? c :
               c == 'b'  ? UInt8('\b') :
               c == 'f'  ? UInt8('\f') :
               c == 'n'  ? UInt8('\n') :
               c == 'r'  ? UInt8('\r') :
               c == 't'  ? UInt8('\t') :
               c == 'u'  ? c : 0x00

"""
Lookup tables for hex digits.
"""
hexvalue(c)      = isxdigit(c) ? tryparse(UInt16, "0x$c") : UInt16(0)
const unhex_1    = [hexvalue(c) for c in Char(1):Char(255)]
const unhex_10   = unhex_1 * 0x10
const unhex_100  = unhex_1 * 0x100
const unhex_1000 = unhex_1 * 0x1000

"""
Read 4HEXDIG from a String `s` at byte index `i` + 1.
FIXME more instructions but less memory access might be faster. Need to test.
"""
function unescape_hex4(s, i)
    i, c = next_ic(s, i) ; u  = @inbounds(unhex_1000[c])
    i, c = next_ic(s, i) ; u += @inbounds(unhex_100[c])
    i, c = next_ic(s, i) ; u += @inbounds(unhex_10[c])
    i, c = next_ic(s, i) ; u += @inbounds(unhex_1[c])
    return i, u
end

"""
Unescape bytes of a String `s` up to byte index `l`, starting at byte index `i`.
Return a new String.
"""
function unescape_json_string(s, i, l)

    i, c = next_ic(s, i)

    out = Base.StringVector(l - i)      # Larger than needed because escape
    j = 1                               # characters will be removed.

    while i < l
        if c != '\\'
            @inbounds out[j] = c
            j += 1
        else
            last_i, cp = json_unescape_char(s, i, c, l)
            if i == last_i
                @inbounds out[j] = cp
                j += 1
            else
                i = last_i
                j = setc_utf8(out, j, cp)
            end
        end
        i, c = next_ic(s, i)
    end

    s = Base.String(out)
    return SubString(s, 1, prevind(s, j))
end

"""
Return a `Char` from a JSON string `s` at index `i`.
Unescpe character if needed.
"""
@propagate_inbounds(
function json_char(s, i)::Tuple{Int, Union{Char,Nothing}}

    c = getc(s, i)
    l = ncodeunits(s)

    if c != '\\' || next_i(s, i) > l
        c = s[i]
        if c == '"'
            c = nothing
        end
        return nextind(s, i), c
    end

    i, c = json_unescape_char(s, i, c, l)
    return next_i(s, i), Char(c)
end)

"""
Return the character code for an escaped character in a JSON string `s`
starting at index `i` (the `\\` character).
"""
function json_unescape_char(s, i, c, l)::Tuple{Int, Union{UInt8,UInt16,UInt32}}

    i, c = next_ic(s, i)
    uc = unescape_c(c)
    if uc == 0x00 ||                    # If the character after '\' was not
       uc == 'u' && next_i(s, i, 4) > l # escapable, or if there aren't enough
        return i, UInt8('\\')           # bytes for \uXXXX, leave the \ as is.
    end

    if uc != 'u'                        # Simple single ecaped character.
        return i, uc                    # Return 7-bit character code.
    end

    i, c16 = unescape_hex4(s, i)        # \uXXXX escaped UTF16 codepoint.
    if c16 in 0xd800:0xdbff          && # Check for surrogate pair...
       next_i(s, i, 6) <= l          &&
       getc(s, next_i(s, i)) == '\\' &&
       getc(s, next_i(s, i, 2)) == 'u'

        j = next_i(s, i, 2)
        j, tail = unescape_hex4(s, j)
        if tail in 0xdc00:0xdfff
            c32 = UInt32(c16 - 0xd7f7)  # Recombine surrogate pair to produce
            c32 <<= 10                  # 32-bit character code.
            c32 += tail
            return j, c32
        end
    end
    return i, c16                       # Return 16-bit character code.
end

"""
Is `i` the index of the start of an escape sequence in a JSON string `s`?
Yes if s[i] == '\\', unless the escape is escaped itself.
"""
function isescape(s, i)
    if getc(s, i) != '\\'               # Check for the escape character.
        return false
    end
    j = i - 1
    while j > 0 && getc(s, j) == '\\'   # Count the number of preceding escapes.
        j -= 1                          # "foo \\\\\uXXXX"
    end                                 #     ^    ^
    return (i - j) % 2 == 1             #     j    i  -- (i - j)=5 % 2 = 1
end

"""
Is `i` the index of the start of a character sequence in a JSON string `s`?
A valid `i` can point to:
 - a single-byte (7-bit) character;
 - a '\\X' escaped 7-bit character;
 - the start of a multi-byte UTF8 character;
 - the start of a `\\uXXXX` escaped UTF16 character;
"""
function string_index_isvalid(s, i)

# FIXME assumes dense indexes

    c = getc(s, i)

    if (c == '"'  ||
        c == '\\' ||
        c == '/'  ||
        c == 'n'  ||
        c == 'r'  ||
        c == 't'  ||
        c == 'b'  ||
        c == 'f'  ||
        c == 'u') &&                    # Escapable character preceded by
        isescape(s, i - 1)              # escape are not a valid index.

        return false
    end

    if c == '\\' &&                     # \uXXXX sequence is invalid if
       getc(s, i + 1) == 'u' &&         # it is preceded by another \uXXXX
       i > 6 &&                         # sequence ...
       getc(s, i - 5) == 'u' &&
       isescape(s, i - 6)

        _, a = unescape_hex4(s, i - 5)
        _, b = unescape_hex4(s, i + 1)

        if a in 0xd800:0xdbff &&        # ... and if the two form a valid UTF16
           b in 0xdc00:0xdfff           # surrogate pair.

            return false
        end
    end

    if ishexdigit(c)                    # A index pointing to a hex digit is
        j = i                           # invalid if part of a \uXXXX sequence.
        m = max(3, i - 3)
        while (c = getc(s, j); ishexdigit(c)) && j >= m
            j -= 1
        end
        if c == 'u' && isescape(s, j - 1)
           return false
        end
    end

    return isvalid(s, i)
end

ishexdigit(c::UInt8) =
    c in UInt8('0'):UInt8('9') ||
    c in UInt8('A'):UInt8('F') ||
    c in UInt8('a'):UInt8('f')

"""
Write a Unicode chatacter `c` into a String `s` as UTF8 at byte index `i`.
"""
function setc_utf8(s::Vector{UInt8}, i, c)
    bytes = bswap(reinterpret(UInt32, Char(c)))
    while true
        @inbounds s[i] = bytes % UInt8
        i += 1
        bytes >>= 8
        if bytes == 0
            return i
        end
    end
end

# IOString Wrappers

#FIXME
#https://github.com/JuliaLang/julia/commit/1f0c6fa35ab64ad66a5bb413fad474e2c722c290#r27686527
@propagate_inbounds function Base.iterate(s::IOString, i::Int = 1)
    b = codeunit(s, i)
    i > ncodeunits(s) && return nothing
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return next_continued(s, i, u)
end

function next_continued(s::IOString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

#FIXME
#https://github.com/JuliaLang/julia/commit/8de25f5ac6c8a8ef9a8872f2d9aaaee9ddbf6bf7#r27684956
Base.isvalid(s::IOString, i::Integer) = Base._thisind_str(s, i) == i

Base.convert(::Type{SubString}, j::JSON.String{IOString{T}}) where T =
    pump(() -> parse_string(j.s, j.i)[1], j.s)

Base.ncodeunits(j::JSON.String{IOString{T}}) where T =
    pump(() -> scan_string(j.s, j.i)[1] - j.i - 1, j.s)

Base.iterate(j::JSON.String{IOString{T}}, i::Integer=1) where T =
    pump(() -> _iterate(j, i), j.s)


text:
@testset "Login" begin
    @test user isa IBMQJulia.IBMQUser
    @test user.id isa String
end

@testset "Backends" begin
    @test reg isa IBMQJulia.IBMQReg
    @test reg.device == "ibmq_qasm_simulator"
    @test reg.id == user.id
end

@testset "job" begin
    @test job isa IBMQJulia.Job
    @test job.reg == reg
    @test job.data isa IBMQJulia.Qobj
    @test job.data.data["experiments"][1]["instructions"] == IBMQJulia.generate_inst(qc)
    @test job.jobid isa String
end

@testset "status" begin
    @test stat == "COMPLETED" || stat == "RUNNING" || stat == "QUEUED"
end

if stat == "COMPLETED"
    @testset "result" begin
        @test res isa Array{Any, 1}
        @test res[1] isa Dict{String, Any}
        @test res[1]["status"] == "DONE"
        @test res[1]["success"] == true
    end
end


text:
using ..HydrostaticBoussinesqModel

function init(model::HydrostaticBoussinesqModel,
              mpicomm::MPI.Comm,
              dg::DGModel,
              Q::MPIStateArray,
              starttime::String)
    @warn "No diagnostics implemented for HydrostaticBoussinesqModel"
end

"""
    collect(bl, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function collect(model::HydrostaticBoussinesqModel, currtime)
    return nothing
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
function iterative_greedy_construction(tuning_run::Run,
                                       reference::RemoteChannel;
                                       cutoff::Integer = 10_000)
    key_set            = collect(keys(tuning_run.starting_point.parameters))
    name               = "Iterative Greedy Construction"
    iteration          = 1
    cost_calls         = tuning_run.cost_evaluations
    stopping_criterion = @task tuning_run.stopping_criterion(tuning_run.duration)
    stop               = consume(stopping_criterion)

    while !stop
        iteration += 1
        for key in key_set
            result                     = greedy_construction(tuning_run,
                                                             key,
                                                             cutoff = cutoff)
            cost_calls                += result.cost_calls
            result.cost_calls          = cost_calls
            result.start               = tuning_run.starting_point
            result.technique           = name
            result.iterations          = iteration
            result.current_iteration   = iteration
            tuning_run.starting_point  = result.minimum
            tuning_run.starting_cost   = result.cost_minimum
            put!(reference, result)
        end
        stop = consume(stopping_criterion)
        shuffle!(key_set)
    end
end


text:
#=
All units in cgs.
=#

module Constants

export G, Msun, kpc, km_per_s

G = 6.67384e-08
Msun = 1.9891e+33
kpc = 3.08568e+21
km_per_s = 1e5
Gyr = 3.15576e+16

end


text:
using QuBase
using QuDynamics
using BenchmarkLite

type SampleProc{Alg} <: Proc end

# Jaynes-Cummings type which takes the configuration parameters.
# These parameters are used for benchmarking, creating various
# JaynesCummings types with various parameters and using
# the solvers to get the timing results.
type JaynesCummings
    N::Int
    wc::Float64
    wa::Float64
    g::Float64
    kappa::Float64
    gamma::Float64
    use_rwa::Bool
    tlist
end

# Name of the solver, used in the creating the table of timing analysis
Base.string{Alg}(::SampleProc{Alg}) = lowercase(string(Alg))

# Constructing a valid configuration, this is set to `true` as we have
# set the type of the required fields in the construct `JaynesCummings` itself.
Base.isvalid{Alg}(p::SampleProc{Alg}, cfg::JaynesCummings) = true

# This returns the `hamitlonian`, `collapse operators`, `initial density matrix`,
# `observable whose expectation value is to be calculated` taking in the `cfg`
# configuration of the system.
function operator_set(p::SampleProc, cfg::JaynesCummings)
    idc = QuArray(eye(cfg.N))
    ida = QuArray(eye(2))
    a  = tensor(lowerop(cfg.N), ida)
    sm = tensor(idc, lowerop(2))
    if cfg.use_rwa
        H = cfg.wc * a' * a + cfg.wa * sm' * sm + cfg.g * (a' * sm + a * sm')
    else
        H = cfg.wc * a' * a + cfg.wa * sm' * sm + cfg.g * (a' + a) * (sm + sm')
    end
    c_op_list = Array(QuBase.AbstractQuMatrix, 0)
    n_th_a = 0.0
    rate_1 = cfg.kappa * (1 + n_th_a)
    if rate_1 > 0.0
        push!(c_op_list, full(sqrt(rate_1) * a))
    end
    rate = cfg.kappa * n_th_a
    if rate > 0.0
        push!(c_op_list, full(sqrt(rate) * a'))
    end
    rate = cfg.gamma
    if rate > 0.0
        push!(c_op_list, full(sqrt(rate) * sm))
    end
    psi = complex(tensor(statevec(1, FiniteBasis(cfg.N)), statevec(2, FiniteBasis(2))))
    return full(H), c_op_list, psi*psi', a*a'
end

# Start states are pre-allocated, that is they are constructed before the evaluation
# here we return the `hamiltonian`, `collapse operator list`, `initial density matrix`,
# `observable whose expectation is to be calculated` using the above function
# i.e., operator_set. The first three remain fixed, the fourth can be varied as required.
Base.start(p::SampleProc, cfg::JaynesCummings) = (operator_set(p, cfg))

Base.length(p::SampleProc, cfg::JaynesCummings) = cfg.N

# This function gets the timing done for calculation of trace, using the pre-allocated
# constructions from the start method i.e., the hamiltonain, collapse operator list,
# initial density matrix and observable whose expectation is to be calculated.
function Base.run{Alg}(p::SampleProc{Alg}, cfg::JaynesCummings, s::(QuBase.QuArray, Array, QuBase.QuArray, QuBase.QuArray))
    qprop = QuPropagator(s[1], s[2], s[3], cfg.tlist, Alg())
    for (t, psi) in qprop
        trace(coeffs(psi*(s[4])))
    end
end

Base.done(p::SampleProc, cfg, s) = nothing

# Creating the array of methods on which benchmarks have to be performed
procs = Proc[ SampleProc{QuExpmV}(),
              SampleProc{QuExpokit}(),
              SampleProc{QuODE45}(),
              SampleProc{QuODE78}(),
              SampleProc{QuODE23s}()]

# Setting up various configurations.
cfgs = JaynesCummings[JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.:0.1:2*pi),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.:0.1:2*pi),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.:0.1:2*pi),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.:0.1:2*pi),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.:0.1:2*pi),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.:0.1:2*pi),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.:0.1:2*pi),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.:0.1:2*pi),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, true, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(2, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(3, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(6, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.0:0.20134228187919462:30.0),
                      JaynesCummings(10, 1.0 * 2 * pi, 1.0 * 2 * pi, 0.05 * 2 * pi, 0.05, 0.15, false, 0.0:0.20134228187919462:30.0)]

# Generating the timing analysis.
rtable = run(procs, cfgs)

# Table display from the above method.
show(rtable; unit=:usec)


text:
module ContinuousDPs

# stdlib
using LinearAlgebra, Random

using QuantEcon
import QuantEcon:
    bellman_operator, bellman_operator!, compute_greedy!, compute_greedy,
    evaluate_policy, DDPAlgorithm, solve, simulate, simulate!

const DPAlgorithm = DDPAlgorithm

include("cdp.jl")
include("lq_approx.jl")

export
    ContinuousDP, evaluate_policy!, set_eval_nodes!, simulate, approx_lq, LQA

end # module


text:
using HoloProcessing
using Test
using FFTW
using BenchmarkTools
using LinearAlgebra: norm

@testset "HoloProcessing.jl" begin
    img = [
        16.0 2 3 13
        5 11 10 8
        9 7 6 12
        4 14 15 1
    ]

    @test reconst(img, plan_rfft(img), 1) |> norm ==
          reconst(img, plan_rfft(img), 1; nthreads = true) |> norm ≈
          1.1215103508786834
    @test reconst(img, plan_fft(img), 1) |> norm ==
          reconst(img, plan_fft(img), 1; nthreads = true) |> norm ≈
          1.1375929179890423
    @test reconst(img, plan_fft(img), 1; shift = true) |> norm ==
          reconst(img, plan_fft(img), 1; shift = true, nthreads = true) |> norm ≈
          1.137592917989042
    @test normalize(img) == normalize!(img)
    @test img |> color .|> Float64 == img
    # path = "/home/qling/Documents/Denoising/hologram"
    # holo = load_holo(path, "head.bmp")
    # FFTW.set_num_threads(Sys.CPU_THREADS)
    # P1 = plan_rfft(holo)
    # @benchmark sdm($holo, $(size(holo) .÷ 2), ($50, $100), $P1, $600)

    # @benchmark rse($holo, $(size(holo) .÷ 2), ($50, $100), $P1, $600)

    # P2 = plan_rfft(similar(holo, size(holo) .÷ 2))

    # @benchmark ldr($holo, $2, ($50, $100), $P2, $600)
end


text:
function optimize(obj::MultivariateOptimizer, arg0::Vector{OptimizationData})
    return jcall(obj, "optimize", PointValuePair, (Vector{OptimizationData},), arg0)
end