text:
using CSV, DataFrames, Distributed, Dates, LinearAlgebra

# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Get number of workers
if length(ARGS) == 1
    addprocs(parse(Int64, ARGS[1]))
elseif length(ARGS) > 1
    throw(ArgumentError("Only one command line argument (cores)."))
end

# Import packages needed for all workers
@everywhere  begin
    using Jedi
    using Distributions
    using DelimitedFiles
end

# Parameters
reps = 200
steps = 5 * 10^8
rho = [0, 0.1, 0.5, 1., 2]
l_0 = 15
N = 1000
nu = 1/N
emat = 2 * (ones(4, 4) - Matrix{Float64}(I, 4, 4))
f0 = 50/2N
fl = 0.3/2N

# Function to run one simulation
@everywhere function run(N, f0, fl, rho, nu, l_0, emat, steps)
    # Initiate population
    pop = mono_pop(N=1000, l=15)
    initiate!(pop, opt=true)
    rand_rho = rand(steps)
    rand_nu = rand(steps)
    f = fermi_fitness(f0=f0, fl=fl)
    Gamma = zeros(Float64, 1000)
    l_arr = zeros(Int64, 1000)
    for i in 1:steps
        bp_substitution!(pop, emat, f)
        if rand_rho[i] < rho/N
            driver_mutation!(pop)
        end
        if rand_nu[i] < 0.001
            l_substitution!(pop, emat, f)
        end
        # Recover lost sites
        if length(pop.seqs) < 7
            initiate!(pop, opt=true)
        end
        if i % 200000 == 0
            Gamma[i ÷ 200000] = get_energy(pop, emat)
            l_arr[i ÷ 200000] = length(pop.seqs)
        end
    end
    return Gamma, l_arr
end
#=
# Run simulations on all available workers
@sync @distributed for j in 1:reps
    for r in 1:length(rho)
        E[r, j, :], L[r, j, :] = run(N, f0, fl, rho[r], nu, l_0, emat, steps)
        RHO[r, j, :] = rho[r] .* ones(Float64, 1000)
    end
    println("Run $j done.")
end

# Save results
df = DataFrame(gamma=[(E...)...], l=[(L...)...], rho=[(RHO...)...])
CSV.write(date*"_script4_results.csv", df)
=#
# Write Metadata
open(date*"_script4_results_METADATA.txt", "a") do io
    write(io, "N=$N\n")
    write(io, "f0=$f0\n")
    write(io, "fl=$fl\n")
    write(io, "repetitions=$reps\n")
    write(io, "steps=$steps\n")
    write(io, "rho=$rho\n")
    write(io, "nu=$nu\n")
    write(io, "l_0=$l_0\n")
end

open(date*"_script4_results.csv", "w") do io
   write(io, "rho\tl\tGamma\n")
end

# Run simulations on all available workers
@sync @distributed for j in 1:reps
    for r in 1:length(rho)
        E, L = run(N, f0, fl, rho[r], nu, l_0, emat, steps)
        RHO = rho[r] .* ones(Float64, 1000)
        open(date*"_script4_results.csv", "a") do io
            writedlm(io, [RHO L E])
        end
    end
    println("Run $j done.")
end


text:
module Findpeaks

export findpeaks

"""
`findpeaks(y::Array{T},
x::Array{S}=collect(1:length(y))
;min_height::T=minimum(y), min_prom::T=minimum(y),
min_dist::S=0, threshold::T=0 ) where {T<:Real,S}`\n
Returns indices of local maxima (sorted from highest peaks to lowest)
in 1D array of real numbers. Similar to MATLAB's findpeaks().\n
*Arguments*:\n
`y` -- data\n
*Optional*:\n
`x` -- x-data\n
*Keyword*:\n
`min_height` -- minimal peak height\n
`min_prom` -- minimal peak prominence\n
`min_dist` -- minimal peak distance (keeping highest peaks)\n
`threshold` -- minimal difference (absolute value) between
 peak and neighboring points\n
"""
function findpeaks(
                   y :: AbstractVector{T},
                   x :: AbstractVector{S} = collect(1:length(y))
                   ;
                   min_height :: T = minimum(y),
                   min_prom :: T = zero(y[1]),
                   min_dist :: S = zero(x[1]),
                   threshold :: T = zero(y[1]),
                  ) where {T <: Real, S}

    dy = diff(y)

    peaks = in_threshold(dy, threshold)

    yP = y[peaks]
    peaks = with_prominence(y, peaks, min_prom)

    #minimal height refinement
    peaks = peaks[y[peaks] .> min_height]
    yP = y[peaks]

    peaks = with_distance(peaks, x, y, min_dist)

    peaks
end

"""
Select peaks that are inside threshold.
"""
function in_threshold(
                      dy :: AbstractVector{T},
                      threshold :: T,
                     ) where {T <: Real}

    peaks = 1:length(dy) |> collect

    k = 0
    for i = 2:length(dy)
        if dy[i] <= -threshold && dy[i-1] >= threshold
            k += 1
            peaks[k] = i
        end
    end
    peaks[1:k]
end

"""
Select peaks that have a given prominence
"""
function with_prominence(
                         y :: AbstractVector{T},
                         peaks :: AbstractVector{Int},
                         min_prom::T,
                        ) where {T <: Real}

    #minimal prominence refinement
    peaks[prominence(y, peaks) .> min_prom]
end

"""
Calculate peaks' prominences
"""
function prominence(y::AbstractVector{T}, peaks::AbstractVector{Int}) where {T <: Real}
    yP = y[peaks]
    proms = zero(yP)

    for (i, p) in enumerate(peaks)
        lP, rP = 1, length(y)
        for j = (i-1):-1:1
            if yP[j] > yP[i]
                lP = peaks[j]
                break
            end
        end
        ml = minimum(y[lP:p])
        for j = (i+1):length(yP)
            if yP[j] > yP[i]
                rP = peaks[j]
                break
            end
        end
        mr = minimum(y[p:rP])
        ref = max(mr,ml)
        proms[i] = yP[i] - ref
    end

    proms
end

"""
Select only peaks that are further apart than `min_dist`
"""
function with_distance(
                       peaks :: AbstractVector{Int},
                       x :: AbstractVector{S},
                       y :: AbstractVector{T},
                       min_dist::S,
                      ) where {T <: Real, S}

    peaks2del = zeros(Bool, length(peaks))
    inds = sortperm(y[peaks], rev=true)
    permute!(peaks, inds)
    for i = 1:length(peaks)
        for j = 1:(i-1)
            if abs(x[peaks[i]] - x[peaks[j]]) <= min_dist
                if !peaks2del[j]
                    peaks2del[i] = true
                end
            end
        end
    end

    peaks[.!peaks2del]
end

end # module

ext: jl


text:
using ParallelAccelerator

@acc function covariance(m::Int, n::Int, float_n::Float64, data::Matrix{Float64})
    cov = Matrix{Float64}(m, m)
    mean = Vector{Float64}(m)
    data_m_mean = Matrix{Float64}(n, m)
    for j in 1:m
        mean[j] = 0.0
        for i = 1:n
            mean[j] += data[i, j]
        end
        mean[j] /= float_n
    end
    for i = 1:n
        for j = 1:m
            data_m_mean[i, j] = data[i, j] - mean[j]
        end
    end
    for i = 1:m
        for j = i:m
            cov[i, j] = 0.0
            for k = 1:n
                cov[i, j] += data_m_mean[k, i] * data_m_mean[k, j]
            end
            cov[i, j] /= (float_n - 1.0)
            cov[j, i] = cov[i, j]
        end
    end
    return cov
end

@inbounds function covariance_o(m::Int, n::Int, float_n::Float64, data::Matrix{Float64})
    cov = Matrix{Float64}(m, m)
    mean = Vector{Float64}(m)
    data_m_mean = Matrix{Float64}(n, m)
    tic()
    for j in 1:m
        mean[j] = 0.0
        for i = 1:n
            mean[j] += data[i, j]
        end
        mean[j] /= float_n
    end
    for i = 1:n
        for j = 1:m
            data_m_mean[i, j] = data[i, j] - mean[j]
        end
    end
    for i = 1:m
        for j = i:m
            cov[i, j] = 0.0
            for k = 1:n
                cov[i, j] += data_m_mean[k, i] * data_m_mean[k, j]
            end
            cov[i, j] /= (float_n - 1.0)
            cov[j, i] = cov[i, j]
        end
    end
    println(" Julia $(toq())")
    return cov
end

m = 2300
n = 2300
float_n = Float64(n)
data = Matrix{Float64}(n, m)
data = rand(n, m)

res_tiramisu = covariance(m, n, float_n, data)
res_julia  = covariance_o(m, n, float_n, data)

function compare(A, B)
    for i in 1:m
        for j in 1:m
            if A[i, j] != B[i, j]
                println("different $((i, j)) : $(A[i, j]) != $(B[i, j])")
                return false
            end
        end
    end
    return true
end

println(compare(res_tiramisu, res_julia))

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
abstract type AbstractBarycentricMatrix{T} <: AbstractLowRankMatrix{T} end

size(B::AbstractBarycentricMatrix) = (B.b-B.a+1, B.d-B.c+1)

struct EvenBarycentricMatrix{T} <: AbstractBarycentricMatrix{T}
    a::Int
    b::Int
    c::Int
    d::Int
    x::Vector{T}
    λ::Vector{T}
    w::Vector{T}
    β::Vector{T}
    W::Matrix{T}
    F::Matrix{T}
end

function EvenBarycentricMatrix(::Type{T}, f::Function, a::Int, b::Int, c::Int, d::Int) where T
    n = BLOCKRANK(T)
    x = chebyshevpoints(T, n)
    λ = chebyshevbarycentricweights(T, n)

    w = zeros(T, b-a+1)
    β = zeros(T, n)
    @inbounds for i = a:b
        for k = 1:n
            w[i+1-a] += λ[k]*inv(2i-a-b-(b-a)*x[k])
        end
    end

    W = zeros(T, n, b-a+1)
    @inbounds for i = a:b
        for k = 1:n
            W[k, i-a+1] = λ[k]*inv((2i-a-b-(b-a)*x[k])*w[i+1-a])
        end
    end

    F = zeros(T, d-c+1, n)
    @inbounds for k = 1:n
        for j = c:d
            F[j-c+1, k] = f(T,(a+b)/2+(b-a)*x[k]/2,j)
        end
    end

    EvenBarycentricMatrix(a, b, c, d, x, λ, w, β, W, F)
end

function getindex(B::EvenBarycentricMatrix{T}, i::Int, j::Int) where T
    ret = zero(T)

    if iseven(size(B, 1)+size(B, 2)+i+j)
        @inbounds for k = 1:length(B.x)
            ret += B.F[j, k]*B.W[k, i]
        end
    end

    ret
end

function barycentricmatrix(::Type{T}, f::Function, a::Int, b::Int, c::Int, d::Int) where T
    n = BLOCKRANK(T)
    x = chebyshevpoints(T, n)
    λ = chebyshevbarycentricweights(T, n)

    w = zeros(T, b-a+1)
    β = zeros(T, n)
    @inbounds for i = a:b
        for k = 1:n
            w[i+1-a] += λ[k]*inv(2i-a-b-(b-a)*x[k])
        end
    end

    W = zeros(T, b-a+1, n)
    @inbounds for k = 1:n
        for i = a:b
            W[i-a+1, k] = λ[k]*inv((2i-a-b-(b-a)*x[k])*w[i+1-a])
        end
    end

    F = zeros(T, d-c+1, n)
    @inbounds for k = 1:n
        for j = c:d
            F[j-c+1, k] = f(T,(a+b)/2+(b-a)*x[k]/2,j)
        end
    end

    LowRankMatrix(W, Diagonal(ones(T, n)), F)
end

function chebyshevpoints(::Type{T}, n::Int; kind::Int = 1) where T
    x = zeros(T, n)
    nd2 = n÷2
    if kind == 1
        @inbounds for k = 1:nd2
            x[k] = sinpi((n-2k+one(T))/2n)
        end
        @inbounds for k=1:nd2
            x[n+1-k] = -x[k]
        end
    else
        @inbounds for k = 1:nd2
            x[k] = sinpi((n-2k+one(T))/(2*(n-1)))
        end
        @inbounds for k=1:nd2
            x[n+1-k] = -x[k]
        end
    end
    x
end

function chebyshevbarycentricweights(::Type{T}, n::Int; kind::Int = 1) where T
    λ = zeros(T, n)
    nd2 = n÷2
    if kind == 1
        @inbounds for k = 1:nd2+1
            λ[k] = sinpi((2k-one(T))/2n)
        end
        @inbounds for k=1:nd2
            λ[n+1-k] = λ[k]
        end
        @inbounds for k=2:2:n
            λ[k] *= -1
        end
    else
        fill!(λ, one(T))
        @inbounds for k=2:2:n
            λ[k] *= -1
        end
        λ[1] *= half(T)
        λ[n] *= half(T)
    end
    λ
end

struct BarycentricPoly2D{T}
    x::Vector{T}
    y::Vector{T}
    λx::Vector{T}
    λy::Vector{T}
    F::Matrix{T}
end

function BarycentricPoly2D(::Type{T}, f::Function, a::T, b::T, c::T, d::T) where T
    r = BLOCKRANK(T)
    rx = r
    ry = r
    x = chebyshevpoints(T, rx)
    y = chebyshevpoints(T, ry)
    λx = chebyshevbarycentricweights(T, rx)
    λy = chebyshevbarycentricweights(T, ry)

    ab2 = half(T)*(a+b)
    ba2 = half(T)*(b-a)
    @inbounds for p in 1:rx
        x[p] = ab2+ba2*x[p]
    end

    cd2 = half(T)*(c+d)
    dc2 = half(T)*(d-c)
    @inbounds for q in 1:ry
        y[q] = cd2+dc2*y[q]
    end

    F = zeros(T, rx, ry)

    @inbounds for n in 1:ry
        yn = y[n]
        for m in 1:rx
            F[m,n] = f(x[m], yn)
        end
    end

    BarycentricPoly2D(x, y, λx, λy, F)
end

function evaluate(B::BarycentricPoly2D{T}, x::T, y::T) where T
    xr, yr, λx, λy, F = B.x, B.y, B.λx, B.λy, B.F
    r = length(xr)
    temp1, temp2, temp3 = zero(T), zero(T), zero(T)
    ret = zero(T)
    for p in 1:r
        temp1 += λx[p]*inv(x-xr[p])
    end
    for q in 1:r
        temp2 += λy[q]*inv(y-yr[q])
    end

    for n in 1:r
        temp3 = zero(T)
        for m in 1:r
            temp3 += λx[m]*F[m,n]*inv(x-xr[m])
        end
        temp3 *= inv(temp1)
        ret += λy[n]*temp3*inv(y-yr[n])
    end
    ret *= inv(temp2)

    ret
end

(B::BarycentricPoly2D)(x, y) = evaluate(B, x, y)

struct BarycentricMatrix2D{T} <: AbstractBarycentricMatrix{T}
    B::BarycentricPoly2D{T}
    x::Vector{T}
    y::Vector{T}
    ir::UnitRange{Int}
    jr::UnitRange{Int}
    U::Matrix{T}
    V::Matrix{T}
    temp1::Vector{T}
    temp2::Vector{T}
end

size(B::BarycentricMatrix2D) = (length(B.ir), length(B.jr))

function getindex(B::BarycentricMatrix2D{T}, i::Int, j::Int) where T
    r = size(B.B.F, 1)
    ret = zero(T)
    for k in 1:r
        temp = zero(T)
        for l in 1:r
            temp += B.B.F[k, l]*B.V[j, l]
        end
        ret += B.U[i,k]*temp
    end

    ret
end

function BarycentricMatrix2D(::Type{T}, f::Function, a::T, b::T, c::T, d::T, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}) where T
    BarycentricMatrix2D(BarycentricPoly2D(T, f, a, b, c, d), x, y, ir, jr)
end

function BarycentricMatrix2D(B::BarycentricPoly2D{T}, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}) where T
    r = length(B.x)
    U = zeros(T, length(ir), r)
    V = zeros(T, length(jr), r)

    update!(BarycentricMatrix2D(B, x, y, ir, jr, U, V, zeros(T, r), zeros(T, r)), x, y, ir, jr)
end

function update!(B::BarycentricMatrix2D{T}, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}) where T
    @assert length(ir) == length(B.ir)
    @assert length(jr) == length(B.jr)
    r = length(B.B.x)
    U, V = B.U, B.V

    ishift = 1-first(ir)
    jshift = 1-first(jr)
    for m in 1:r
        λxm = B.B.λx[m]
        xm = B.B.x[m]
        for i in ir
            B.U[i+ishift, m] = λxm*inv(x[i]-xm)
        end
    end
    for i in ir
        tempU = zero(T)
        for m in 1:r
            tempU += U[i+ishift, m]
        end
        for m in 1:r
            U[i+ishift, m] /= tempU
        end
    end
    for n in 1:r
        λyn = B.B.λy[n]
        yn = B.B.y[n]
        for j in jr
            V[j+jshift, n] = λyn*inv(y[j]-yn)
        end
    end
    for j in jr
        tempV = zero(T)
        for n in 1:r
            tempV += V[j+jshift, n]
        end
        for n in 1:r
            V[j+jshift, n] /= tempV
        end
    end

    for i in ir
        B.x[i] = x[i]
    end
    for j in jr
        B.y[j] = y[j]
    end

    B
end

function indsplit(x::Vector{T}, ir::UnitRange{Int}, a::T, b::T) where T
    i = first(ir)
    ab2 = half(T)*(a+b)
    while x[i] ≥ ab2
        i += 1
        i > last(ir) && break
    end
    first(ir):i-1, i:last(ir)
end


text:
module Ranges

include("positional.jl")
include("pair.jl")
include("partial.jl")
include("keyword.jl")
include("length_step.jl")

# Alternate module for alternate positional syntax
# include("alternate.jl")

end


text:
#######################################################################################
#                                 Load Packages
#######################################################################################
# set the working directory to the directory in which this file is contained
cd(@__DIR__)
# load the package manager
using Pkg
# activate the project environment
Pkg.activate("../../")
using Revise, AdaptiveDesignOptimization, Random, UtilityModels, Distributions
includet("TAX_Model.jl")
Random.seed!(25974)
#######################################################################################
#                               Define Experiment Design
#######################################################################################
# model with default uniform prior
model = Model(;loglike)

parm_list = (
    δ = range(-2, 2, length=10),
    β = range(.5, 1.5, length=10),
    γ = range(.5, 1.2, length=10),
    θ = range(.5, 3, length=10)
)

dist = Normal(0,10)
n_vals = 3
n_choices = 2
design_vals = map(x->random_design(dist, n_vals, n_choices), 1:1000)
filter!(x->abs_zscore(x) ≤ .4, design_vals)
design_names = (:p1,:v1,:p2,:v2)
design_list = (design_names,design_vals[1:100])

data_list = (choice=[true, false],)
#######################################################################################
#                              Optimize Experiment
#######################################################################################
using DataFrames
true_parms = (δ=-1.0, β=1.0, γ=.7, θ=1.5)
n_trials = 100
optimizer = Optimizer(;design_list, parm_list, data_list, model)
design = optimizer.best_design
df = DataFrame(design=Symbol[], trial=Int[], mean_δ=Float64[], mean_β=Float64[],
    mean_γ=Float64[], mean_θ=Float64[], std_δ=Float64[], std_β=Float64[],
    std_γ=Float64[], std_θ=Float64[])
new_data = [:optimal, 0, mean_post(optimizer)..., std_post(optimizer)...]
push!(df, new_data)

for trial in 1:n_trials
    data = simulate(true_parms..., design...)
    design = update!(optimizer, data)
    new_data = [:optimal, trial, mean_post(optimizer)..., std_post(optimizer)...]
    push!(df, new_data)
end
#######################################################################################
#                              Random Experiment
#######################################################################################
randomizer = Optimizer(;design_list, parm_list, data_list, model, design_type=Randomize);
design = randomizer.best_design
new_data = [:random, 0, mean_post(randomizer)..., std_post(randomizer)...]
push!(df, new_data)

for trial in 1:n_trials
    data = simulate(true_parms..., design...)
    design = update!(randomizer, data)
    new_data = [:random, trial, mean_post(randomizer)..., std_post(randomizer)...]
    push!(df, new_data)
end
#######################################################################################
#                                 Plot Results
#######################################################################################
using StatsPlots
@df df plot(:trial, :mean_δ, xlabel="trial", ylabel="mean δ", group=:design, grid=false)
hline!([true_parms.δ], label="true")

@df df plot(:trial, :mean_β, xlabel="trial", ylabel="mean β", group=:design, grid=false)
hline!([true_parms.β], label="true")

@df df plot(:trial, :mean_γ, xlabel="trial", ylabel="mean γ", group=:design, grid=false)
hline!([true_parms.γ], label="true")

@df df plot(:trial, :mean_θ, xlabel="trial", ylabel="mean θ", group=:design, grid=false)
hline!([true_parms.θ], label="true")

@df df plot(:trial, :std_δ, xlabel="trial", ylabel="σ of δ", grid=false, group=:design, ylims=(0,1.5))

@df df plot(:trial, :std_β, xlabel="trial", ylabel="σ of β", grid=false, group=:design, ylims=(0,.5))

@df df plot(:trial, :std_γ, xlabel="trial", ylabel="σ of γ", grid=false, group=:design, ylims=(0,.5))

@df df plot(:trial, :std_θ, xlabel="trial", ylabel="σ of θ", grid=false, group=:design, ylims=(0,.8))


text:
using PyPlot, SparseArrays, LinearAlgebra # Packages used

rows = [1,3,4,2,1,3,1,4,1,5]
cols = [1,1,1,2,3,3,4,4,5,5]
vals = [5,-2,-4,5,-3,-1,-2,-10,7,9]

A = sparse(rows, cols, vals, 5, 5)

Array(A)

spy(A, marker=".", markersize=24);   ## Note - 0-based row and columns

B = A  - 4.3A   # Will automatically convert datatype of values to Float64
B[:,4] .= -1.1  # OK since B now has Float64 values (otherwise use Float64.(A) to convert)
C = A * A'      # Matrix multiplication (note: typically increases nnz)
Matrix([B C])   # Concatenation, again automatic conversion (of C)

"""
Incremental matrix construction using the sparse-format
Not recommended: Insertion into existing matrix very slow
"""
function incremental_test_1(N)
    A = spzeros(N,N)
    for k = 1:10N
        i,j = rand(1:N, 2)
        A[i,j] = rand()
    end
    return A
end

incremental_test_1(10);   # Force compile before timing
for N in [100,1000,10000]
  @time incremental_test_1(N);
end

"""
Incremental matrix construction using COO and a single call to sparse
Fast approach, avoids incremental insertion into existing array
"""
function incremental_test_2(N)
    rows = Int64[]
    cols = Int64[]
    vals = Float64[]
    for i = 1:10N
        push!(rows, rand(1:N))
        push!(cols, rand(1:N))
        push!(vals, rand())
    end
    return sparse(rows, cols, vals, N, N)
end

incremental_test_2(10);   # Force compile before timing
for N in [100,1000,10000,100000,1000000]
  @time incremental_test_2(N);
end


text:
if 1==1
    using EconPDEs, Plots
    stategrid = OrderedDict(:k => range(-7.0, 7.0, length = 1000))
    solend = OrderedDict(:V => ones(1000))
    function f(state::NamedTuple, sol::NamedTuple)
        #ρ = 0.05; δ = 0.10; z = 0.20;
        ρ = 0.15; δ = 0.10; z = 0.27;
        k = state.k
        V, Vk_up, Vk_down, Vkk = sol.V, sol.Vk_up, sol.Vk_down, sol.Vkk
        #
        Vk = Vk_up
        iter = 0
        @label start
        i = Vk - 1.0
        μk = i - δ*k
        if (iter == 0) & (μk <= 0)
            iter += 1
            Vk = Vk_down
            @goto start
        end
        Vt = - (z*k - i -0.5*(i^2.0) + μk*Vk - ρ*V)
        (Vt = Vt,)
    end
    y, residual_norm =  pdesolve(f, stategrid, solend)
    plot(y[:V])
    #To Do: simulate
end


text:
using DifferentialEquations
using Plots

μ = 0.15 # Drift
σ = 1.0 # Vol of Vol
θ = 0.04 # long-term variance
κ = 1.0 # Mean reversion speed
ρ = -0.75 # spot-vol correlation
u0 = [100.0, 0.0225]
tspan = (0.0, 1.0)
dt = 1.0/365
prob = HestonProblem(μ, κ, θ, σ, ρ, u0, (0.0, 1.0))

sol = solve(prob, EM(), dt=dt)

print(sol.u[1][1])

path_matrix = hcat(sol.u)

print(typeof(path_matrix))
# print(path_matrix)

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
#=
"""
    DiagnosticsGroupParams

Base type for any extra parameters needed by a diagnostics group.
"""
abstract type DiagnosticsGroupParams end

"""
    DiagnosticsGroup

A diagnostics group contains a set of diagnostic variables that are
collected at the same interval and written to the same file. One or more
diagnostics groups are placed in a [`DiagnosticsConfiguration`](@ref)
which is used by [`ClimateMachine.invoke!`](@ref).

A `DiagnosticsGroup` can be easily created with
[`@diagnostics_group`](@ref).
"""
mutable struct DiagnosticsGroup{
    DGP <: Union{Nothing, DiagnosticsGroupParams},
    DGI <: Union{Nothing, InterpolationTopology},
}
    name::String
    init::Function
    collect::Function
    fini::Function
    interval::String
    out_prefix::String
    writer::AbstractWriter
    interpol::DGI
    onetime::Bool
    params::DGP

    DiagnosticsGroup(
        name,
        init,
        fini,
        collect,
        interval,
        out_prefix,
        writer,
        interpol,
        onetime,
        params,
    ) = new{typeof(params), typeof(interpol)}(
        name,
        init,
        collect,
        fini,
        interval,
        out_prefix,
        writer,
        interpol,
        onetime,
        params,
    )
end

# `GenericCallbacks` implementations for `DiagnosticsGroup`s
function GenericCallbacks.init!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @info @sprintf(
        """
    Diagnostics: %s
        initializing at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.init(dgngrp, t)
    dgngrp.collect(dgngrp, t)
    return nothing
end
function GenericCallbacks.call!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @tic diagnostics
    @info @sprintf(
        """
    Diagnostics: %s
        collecting at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.collect(dgngrp, t)
    @toc diagnostics
    return nothing
end
function GenericCallbacks.fini!(dgngrp::DiagnosticsGroup, solver, Q, param, t)
    @info @sprintf(
        """
    Diagnostics: %s
        finishing at %8.2f""",
        dgngrp.name,
        t,
    )
    dgngrp.collect(dgngrp, t)
    dgngrp.fini(dgngrp, t)
    return nothing
end
=#
abstract type InterpolationType end
struct NoInterpolation <: InterpolationType end
struct InterpolateAfterCollection <: InterpolationType end
struct InterpolateBeforeCollection <: InterpolationType end

"""
    @diagnostics_group

Generate the functions needed to establish and
use a [`DiagnosticsGroup`](@ref) containing the named
[`DiagnosticVar`](@ref)s. In particular, this creates `$(name)()`,
which creates the diagnostics group.

# Arguments
- `name`: a string that uniquely identifies the group.
- `config_type`: a `ClimateMachineConfigType`.
- `params_type`: a subtype of `DiagnosticsGroupParams` or `Nothing`.
  An instance of this type is created in `dgngrp.params` when the
  group is set up.
- `init_fun`: a function that is called when the group is initialized
  (called with `(dgngrp, curr_time); `may be `(_...) -> nothing`). Use
  this to initialize any required state, such as in `dgngrp.params`.
- `interpolate`: one of `InterpolateBeforeCollection`,
  `InterpolateAfterCollection` or `NoInterpolation`.
- `dvarnames`: one or more diagnostic variables. Together with the
  `config_type`, identify the [`DiagnosticVar`](@ref)s to be included
  in the group.
"""
macro diagnostics_group(
    name,
    config_type,
    params_type,
    init_fun,
    interpolate,
    dvarnames..., # TODO: need to separately specify which ones to output
)
    CT = getfield(ConfigTypes, config_type)
    dvars = DiagnosticVar[AllDiagnosticVars[CT][String(dv)] for dv in dvarnames]

    # Partition the diagnostic variables by type.
    dvtype_dvars_map = Dict{DataType, Array{DiagnosticVar, 1}}()
    for dvar in dvars
        push!(
            get!(dvtype_dvars_map, supertype(typeof(dvar)), DiagnosticVar[]),
            dvar,
        )
    end

    gen_params = (
        name,
        config_type,
        params_type,
        init_fun,
        interpolate,
        dvtype_dvars_map,
    )

    using_exprs = quote
        using KernelAbstractions
        using OrderedCollections
        using ClimateMachine.BalanceLaws
        using ClimateMachine.DGMethods
        using ClimateMachine.Mesh.Interpolation
        using ClimateMachine.Mesh.Topologies
        using ClimateMachine.MPIStateArrays
        using ClimateMachine.VariableTemplates
        using ClimateMachine.Writers
    end
    init_fun = esc(prewalk(unblock, generate_init_fun(gen_params...)))
    fini_fun = esc(prewalk(unblock, generate_fini_fun(gen_params...)))
    collect_fun = esc(prewalk(unblock, generate_collect_fun(gen_params...)))
    setup_fun = esc(prewalk(unblock, generate_setup_fun(gen_params...)))

    return Expr(:block, using_exprs, init_fun, fini_fun, collect_fun, setup_fun)
end

include("group_gen.jl")


text:
@testset "Function counter" begin
    prob = OptimTestProblems.UnconstrainedProblems.examples["Rosenbrock"]

    let
        global fcount = 0
        global fcounter
        function fcounter(reset::Bool = false)
            if reset
                fcount = 0
            else
                fcount += 1
            end
            fcount
        end
        global gcount = 0
        global gcounter
        function gcounter(reset::Bool = false)
            if reset
                gcount = 0
            else
                gcount += 1
            end
            gcount
        end
        global hcount = 0
        global hcounter
        function hcounter(reset::Bool = false)
            if reset
                hcount = 0
            else
                hcount += 1
            end
            hcount
        end
    end

    f(x) = begin
        fcounter()
        MVP.objective(prob)(x)
    end
    g!(out, x) = begin
        gcounter()
        MVP.gradient(prob)(out, x)
    end
    h!(out, x) = begin
        hcounter()
        MVP.hessian(prob)(out, x)
    end

    options = Optim.Options(; Optim.default_options(IPNewton())...)
    # TODO: Run this on backtrack_constrained as well (when we figure out what it does)
    for ls in [Optim.backtrack_constrained_grad,]
               #Optim.backtrack_constrained]

        fcounter(true); gcounter(true); hcounter(true)

        df = TwiceDifferentiable(f, g!, h!, prob.initial_x)
        infvec = fill(Inf, size(prob.initial_x))
        # TODO: fix `optimize`so that empty constraints work
        #constraints = TwiceDifferentiableConstraints([], [])
        constraints = TwiceDifferentiableConstraints(-infvec, infvec)

        res = optimize(df, constraints, prob.initial_x,
                       IPNewton(linesearch = ls),
                       options)
        @test fcount == Optim.f_calls(res)
        @test gcount == Optim.g_calls(res)
        @test hcount == Optim.h_calls(res)
    end

    # TODO: Test a constrained problem

    # TODO: Test constraints call counter (if we need to make one)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Extract data from a 2D/3D DG solution and prepare it for visualization as a heatmap/contour plot.
#
# Returns a tuple with
# - x coordinates
# - y coordinates
# - nodal 2D data
# - x vertices for mesh visualization
# - y vertices for mesh visualization
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, levels, ndims,
                     unstructured_data, n_nodes,
                     grid_lines=false, max_supported_level=11, nvisnodes=nothing,
                     slice=:xy, point=(0.0, 0.0, 0.0))
  # Determine resolution for data interpolation
  max_level = maximum(levels)
  if max_level > max_supported_level
    error("Maximum refinement level $max_level is higher than " *
          "maximum supported level $max_supported_level")
  end
  max_available_nodes_per_finest_element = 2^(max_supported_level - max_level)
  if nvisnodes === nothing
    max_nvisnodes = 2 * n_nodes
  elseif nvisnodes == 0
    max_nvisnodes = n_nodes
  else
    max_nvisnodes = nvisnodes
  end
  nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes)
  resolution = nvisnodes_at_max_level * 2^max_level
  nvisnodes_per_level = [2^(max_level - level)*nvisnodes_at_max_level for level in 0:max_level]
  # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate
  # level-0-cell) that contains the number of visualization nodes for any
  # refinement level to visualize on an equidistant grid

  if ndims == 3
    (unstructured_data, coordinates, levels,
        center_level_0) = unstructured_3d_to_2d(unstructured_data,
        coordinates, levels, length_level_0, center_level_0, slice,
        point)
  end

  # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
  n_elements = length(levels)
  normalized_coordinates = similar(coordinates)
  for element_id in 1:n_elements
    @views normalized_coordinates[:, element_id] .= (
          (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 ))
  end

  # Interpolate unstructured DG data to structured data
  (structured_data =
      unstructured2structured(unstructured_data, normalized_coordinates,
                              levels, resolution, nvisnodes_per_level))

  # Interpolate cell-centered values to node-centered values
  node_centered_data = cell2node(structured_data)

  # Determine axis coordinates for contour plot
  xs = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[1]
  ys = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[2]

  # Determine element vertices to plot grid lines
  if grid_lines
    mesh_vertices_x, mesh_vertices_y = calc_vertices(coordinates, levels, length_level_0)
  else
    mesh_vertices_x = Matrix{Float64}(undef, 0, 0)
    mesh_vertices_y = Matrix{Float64}(undef, 0, 0)
  end

  return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y
end

# Extract data from a 1D DG solution and prepare it for visualization as a line plot.
# This returns a tuple with
# - x coordinates
# - nodal 1D data
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_data_1d(original_nodes, unstructured_data, nvisnodes)
  # Get the dimensions of u; where n_vars is the number of variables, n_nodes the number of nodal values per element and n_elements the total number of elements.
  n_nodes, n_elements, n_vars = size(unstructured_data)

  # Set the amount of nodes visualized according to nvisnodes.
  if nvisnodes === nothing
    max_nvisnodes = 2 * n_nodes
  elseif nvisnodes == 0
    max_nvisnodes = n_nodes
  else
    @assert nvisnodes >= 2 "nvisnodes must be zero or >= 2"
    max_nvisnodes = nvisnodes
  end

  interpolated_nodes = Array{eltype(original_nodes),    2}(undef, max_nvisnodes, n_elements)
  interpolated_data  = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes, n_elements, n_vars)

  for j in 1:n_elements
    # Interpolate on an equidistant grid.
    interpolated_nodes[:, j] .= range(original_nodes[1,1,j], original_nodes[1,end,j], length = max_nvisnodes)
  end

  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)
  nodes_out = collect(range(-1, 1, length = max_nvisnodes))

  # Calculate vandermonde matrix for interpolation.
  vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)

  # Iterate over all variables.
  for v in 1:n_vars
    # Interpolate data for each element.
    for element in 1:n_elements
      multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]),
        vandermonde, @view(unstructured_data[:, element, v]))
    end
  end
  # Return results after data is reshaped
  return vec(interpolated_nodes), reshape(interpolated_data, :, n_vars), vcat(original_nodes[1, 1, :], original_nodes[1, end, end])
end

# Change order of dimensions (variables are now last) and convert data to `solution_variables`
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_unstructured_data(u, solution_variables, mesh, equations, solver, cache)

  if solution_variables === cons2cons
    raw_data = u
    n_vars = size(raw_data, 1)
  else
    # FIXME: Remove this comment once the implementation following it has been verified
    # Reinterpret the solution array as an array of conservative variables,
    # compute the solution variables via broadcasting, and reinterpret the
    # result as a plain array of floating point numbers
    # raw_data = Array(reinterpret(eltype(u),
    #        solution_variables.(reinterpret(SVector{nvariables(equations),eltype(u)}, u),
    #                   Ref(equations))))
    # n_vars = size(raw_data, 1)
    n_vars_in = nvariables(equations)
    n_vars = length(solution_variables(get_node_vars(u, equations, solver), equations))
    raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...)
    reshaped_u = reshape(u, n_vars_in, :)
    reshaped_r = reshape(raw_data, n_vars, :)
    for idx in axes(reshaped_u, 2)
      reshaped_r[:, idx] = solution_variables(get_node_vars(reshaped_u, equations, solver, idx), equations)
    end
  end

  unstructured_data = Array{eltype(raw_data)}(undef,
                                              ntuple((d) -> nnodes(solver), ndims(equations))...,
                                              nelements(solver, cache), n_vars)
  for variable in 1:n_vars
    @views unstructured_data[.., :, variable] .= raw_data[variable, .., :]
  end

  return unstructured_data
end

end # @muladd


text:
function testDCT1d(N=32)
  Random.seed!(1235)
  x = zeros(ComplexF64, N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2)) .+ 1im*rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = DCTOp(Float64,(N^2,),2)
  D2 = sqrt(2/N^2)*[cos(pi/(N^2)*j*(k+0.5)) for j=0:N^2-1,k=0:N^2-1]
  D2[1,:] .*= 1/sqrt(2)
  D3 = DCTOp(Float64,(N^2,),4)
  D4 = sqrt(2/N^2)*[cos(pi/(N^2)*(j+0.5)*(k+0.5)) for j=0:N^2-1,k=0:N^2-1]

  y1 = D1*x
  y2 = D2*x

  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  y3 = D3*x
  y4 = D4*x
  @test norm(y3 - y4) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1)*y1
  @test norm(x1 - x) / norm(x) ≈ 0 atol=0.01

  x2 = adjoint(D3)*y3
  @test norm(x2 - x) / norm(x) ≈ 0 atol=0.01
end

function testFFT1d(N=32,shift=true)
  Random.seed!(1234)
  x = zeros(N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = FFTOp(Float64,(N^2,),shift)
  D2 =  1.0/N*[exp(-2*pi*im*j*k/N^2) for j=0:N^2-1,k=0:N^2-1]

  y1 = D1*x
  if shift
    y2 = fftshift(D2*fftshift(x))
  else
    y2 = D2*x
  end
  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1) * y1
  if shift
    x2 = ifftshift(adjoint(D2)*ifftshift(y2))
  else
    x2 = adjoint(D2)*y2
  end
  @test norm(x1 - x2) / norm(x1) ≈ 0 atol=0.01
end

function testFFT2d(N=32,shift=true)
  Random.seed!(1234)
  x = zeros(N^2)
  for i=1:5
    x .+= rand()*cos.(rand(1:N^2)*collect(1:N^2))
  end
  D1 = FFTOp(Float64,(N,N),shift)

  idx = CartesianIndices((N,N))[collect(1:N^2)]
  D2 =  1.0/N*[ exp(-2*pi*im*((idx[j][1]-1)*(idx[k][1]-1)+(idx[j][2]-1)*(idx[k][2]-1))/N) for j=1:N^2, k=1:N^2 ]

  y1 = D1*x
  if shift
    y2 = D2*vec(fftshift(reshape(x,N,N)))
    y2 = vec(fftshift(reshape(y2,N,N)))
  else
    y2 = D2*x
  end
  @test norm(y1 - y2) / norm(y1) ≈ 0 atol=0.01

  x1 = adjoint(D1) * y1
  if shift
    x2 = adjoint(D2)*vec(ifftshift(reshape(y2,N,N)))
    x2 = vec(ifftshift(reshape(x2,N,N)))
  else
    x2 = adjoint(D2)*y2
  end
  @test norm(x1 - x2) / norm(x1) ≈ 0 atol=0.01
end

function testWeighting(N=512)
  Random.seed!(1234)
  x1 = rand(N)
  weights = rand(N)
  W = WeightingOp(weights)
  y1 = W*x1
  y = weights .* x1

  @test norm(y1 - y) / norm(y) ≈ 0 atol=0.01

  x2 = rand(2*N)
  W2 = WeightingOp(weights,2)
  y2 = W2*x2
  y = repeat(weights,2) .* x2

  @test norm(y2 - y) / norm(y) ≈ 0 atol=0.01
end

@testset "Linear Operators" begin
  @info "test DCT-II and DCT-IV"
  for N in [2,8,16,32]
    testDCT1d(N)
  end
  @info "test FFT"
  for N in [8,16,32]
    testFFT1d(N,false)
    testFFT1d(N,true)
    testFFT2d(N,false)
    testFFT2d(N,true)
  end
  @info "test Weighting"
  testWeighting(512)
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# Node.jl: helper function to access interface to Camera nodes

# Utility functions
function available(nodeName)
  pbAvailable = Ref(bool8_t(false))
  spinNodeIsAvailable(nodeName[], pbAvailable)
  return (pbAvailable[] == 1)
end

function readable(nodeName)
  if available(nodeName)
   pbReadable = Ref(bool8_t(false))
   spinNodeIsReadable(nodeName[], pbReadable)
   return (pbReadable[] == 1)
 else
   return false
 end
end

function writable(nodeName)
  if available(nodeName)
   pbWriteable = Ref(bool8_t(false))
   spinNodeIsReadable(nodeName[], pbWriteable)
   return (pbWriteable[] == 1)
 else
   return false
 end
end

function _getnode(cam, name::String, nodemap)
  hNodeMap = Ref(spinNodeMapHandle(C_NULL))
  _nodemap!(cam, hNodeMap, nodemap)
  hNode = Ref(spinNodeHandle(C_NULL))
  spinNodeMapGetNode(hNodeMap[], name, hNode);
  return hNode
end

abstract type AbstractSpinNode end

#
# String nodes
#

struct SpinStringNode <: AbstractSpinNode
  name::String
  hNode::Ref{spinNodeHandle}
  SpinStringNode(cam, name::String, nodemap=CameraNodeMap()) = new(name, _getnode(cam, name, nodemap))
end

function get(node::SpinStringNode)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  strbuf = Vector{UInt8}(undef, MAX_BUFFER_LEN)
  strlen = Ref(Csize_t(MAX_BUFFER_LEN))
  spinStringGetValue(node.hNode[], strbuf, strlen)
  return unsafe_string(pointer(strbuf))
end

#
# Integer nodes
#

struct SpinIntegerNode <: AbstractSpinNode
  name::String
  hNode::Ref{spinNodeHandle}
  SpinIntegerNode(cam, name::String, nodemap=CameraNodeMap()) = new(name, _getnode(cam, name, nodemap))
end

function range(node::SpinIntegerNode)
  hMin = Ref(Int64(0.0))
  hMax = Ref(Int64(0.0))
  spinIntegerGetMin(node.hNode[], hMin)
  spinIntegerGetMax(node.hNode[], hMax)
  return (hMin[], hMax[])
end

function set!(node::SpinIntegerNode, value::Number; clampwarn::Bool = true)
  if !writable(node.hNode)
    throw(ErrorException("Node $(node.name) is not writable"))
  end
  noderange = range(node)
  if clampwarn
    value < noderange[1] && @warn "Requested value ($value) is smaller than minimum ($(noderange[1])), value will be clamped."
    value > noderange[2] && @warn "Requested value ($value) is greater than maximum ($(noderange[2])), value will be clamped."
  end
  spinIntegerSetValue(node.hNode[], Int64(clamp(value, noderange[1], noderange[2])))
  get(node)
end

function get(node::SpinIntegerNode)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  hval = Ref(Int64(0))
  spinIntegerGetValue(node.hNode[], hval)
  return hval[]
end

#
# Floating point nodes
#

struct SpinFloatNode <: AbstractSpinNode
  name::String
  hNode::Ref{spinNodeHandle}
  SpinFloatNode(cam, name::String, nodemap=CameraNodeMap()) = new(name, _getnode(cam, name, nodemap))
end

function range(node::SpinFloatNode)
  hMin = Ref(Float64(0.0))
  hMax = Ref(Float64(0.0))
  spinFloatGetMin(node.hNode[], hMin)
  spinFloatGetMax(node.hNode[], hMax)
  return (hMin[], hMax[])
end

function set!(node::SpinFloatNode, value::Number; clampwarn::Bool = true)
  if !writable(node.hNode)
    throw(ErrorException("Node $(node.name) is not writable"))
  end
  noderange = range(node)
  if clampwarn
    value < noderange[1] && @warn "Requested value ($value) is smaller than minimum ($(noderange[1])), value will be clamped."
    value > noderange[2] && @warn "Requested value ($value) is greater than maximum ($(noderange[2])), value will be clamped."
  end
  spinFloatSetValue(node.hNode[], Float64(clamp(value, noderange[1], noderange[2])))
  get(node)
end

function get(node::SpinFloatNode)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  hval = Ref(Float64(0))
  spinFloatGetValue(node.hNode[], hval)
  return hval[]
end

#
# Integer enumeration nodes
#

struct SpinEnumNode <: AbstractSpinNode
  name::String
  hNode::Ref{spinNodeHandle}
  SpinEnumNode(cam, name::String, nodemap=CameraNodeMap()) = new(name, _getnode(cam, name, nodemap))
end

function get(node::SpinEnumNode)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  hNodeEntry = Ref(spinNodeHandle(C_NULL))
  strbuf = Vector{UInt8}(undef, MAX_BUFFER_LEN)
  strlen = Ref(Csize_t(MAX_BUFFER_LEN))
  spinEnumerationGetCurrentEntry(node.hNode[], hNodeEntry)
  spinEnumerationEntryGetSymbolic(hNodeEntry[], strbuf, strlen)
  return unsafe_string(pointer(strbuf))
end

function set!(node::SpinEnumNode, value)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  hNodeEntry = Ref(spinNodeHandle(C_NULL))
  hNodeVal = Ref(Int64(0))
  spinEnumerationGetEntryByName(node.hNode[], value, hNodeEntry)

  # Get integer value from string
  if !readable(hNodeEntry)
    throw(ErrorException("Node $(node.name) entry is not readable"))
  end
  spinEnumerationEntryGetIntValue(hNodeEntry[], hNodeVal)

  # Set value
  if !writable(node.hNode)
    throw(ErrorException("Node $(node.name) is not writable"))
  end
  spinEnumerationSetIntValue(node.hNode[], hNodeVal[])

  # Readback
  get(node)
end

#
# Boolean mnodes
#

struct SpinBooleanNode <: AbstractSpinNode
  name::String
  hNode::Ref{spinNodeHandle}
  SpinBooleanNode(cam, name::String, nodemap=CameraNodeMap()) = new(name, _getnode(cam, name, nodemap))
end

function get(node::SpinBooleanNode)
  if !readable(node.hNode)
    throw(ErrorException("Node $(node.name) is not readable"))
  end
  hval = Ref(bool8_t(0))
  spinBooleanGetValue(node.hNode[], hval)
  return hval[] == 1 ? true : false
end

function set!(node::SpinBooleanNode, value::Bool)
  if !writable(node.hNode)
    throw(ErrorException("Node $(node.name) is not writable"))
  end
  spinBooleanSetValue(node.hNode[], value)
  get(node)
end


text:
__precompile__(true)

module GeneralQP

using LinearAlgebra
using Polynomials

include("linear_algebra.jl")
include("change_constraints.jl")
include("qp.jl")
include("printing.jl")
export solve
export UpdatableQR, NullspaceHessianLDL, NullspaceHessian
export add_column!, remove_column!
export add_constraint!, remove_constraint!

end


text:
""" Lay out and draw wiring diagrams using Graphviz.

This module requires Graphviz v2.42 or higher.
"""
module GraphvizWiringDiagrams
export to_graphviz, graphviz_layout

using Compat: isnothing
import JSON
using LinearAlgebra: normalize
using MLStyle: @match
using StaticArrays

import ...Theories: HomExpr
using ...WiringDiagrams, ...WiringDiagrams.WiringDiagramSerialization
using ...CategoricalAlgebra.CSets, ...Graphs, ..GraphvizGraphs
import ..Graphviz
import ..GraphvizGraphs: to_graphviz
using ..WiringDiagramLayouts: BoxLayout, PortLayout, WirePoint,
  LayoutOrientation, LeftToRight, RightToLeft, TopToBottom, BottomToTop,
  is_horizontal, is_vertical, box_label, wire_label, port_sign, svector

# Default Graphviz font. Reference: http://www.graphviz.org/doc/fontfaq.txt
const default_graphviz_font = "Serif"

# Directed drawing
##################

# Default graph, node, edge, and cell attributes for directed wiring diagrams.
const default_directed_graphviz_attrs = (
  graph = Graphviz.Attributes(
    :fontname => default_graphviz_font,
  ),
  node = Graphviz.Attributes(
    :fontname => default_graphviz_font,
    :shape => "none",
    :width => "0",
    :height => "0",
    :margin => "0",
  ),
  edge = Graphviz.Attributes(
    :arrowsize => "0.5",
    :fontname => default_graphviz_font,
  ),
  cell = Graphviz.Attributes(
    :border => "1",
    :cellpadding => "4",
  )
)

struct GraphvizBox
  stmts::Vector{Graphviz.Statement} # Usually Graphviz.Node
  input_ports::Vector{Graphviz.NodeID}
  output_ports::Vector{Graphviz.NodeID}
end

""" Draw a wiring diagram using Graphviz.

The input can also be a morphism expression, in which case it is first converted
into a wiring diagram. This function requires Graphviz v2.42 or higher.

# Arguments
- `graph_name="G"`: name of Graphviz digraph
- `orientation=TopToBottom`: orientation of layout.
  One of `LeftToRight`, `RightToLeft`, `TopToBottom`, or `BottomToTop`.
- `node_labels=true`: whether to label the nodes
- `labels=false`: whether to label the edges
- `label_attr=:label`: what kind of edge label to use (if `labels` is true).
  One of `:label`, `:xlabel`, `:headlabel`, or `:taillabel`.
- `port_size="24"`: minimum size of ports on box, in points
- `junction_size="0.05"`: size of junction nodes, in inches
- `outer_ports=true`: whether to display the outer box's input and output ports.
  If disabled, no incoming or outgoing wires will be shown either!
- `anchor_outer_ports=true`: whether to enforce ordering of the outer box's
  input and output, i.e., ordering of the incoming and outgoing wires
- `graph_attrs=Dict()`: top-level graph attributes
- `node_attrs=Dict()`: top-level node attributes
- `edge_attrs=Dict()`: top-level edge attributes
- `cell_attrs=Dict()`: main cell attributes in node HTML-like label
"""
function to_graphviz(f::WiringDiagram;
    graph_name::String="G", orientation::LayoutOrientation=TopToBottom,
    node_labels::Bool=true, labels::Bool=false, label_attr::Symbol=:label,
    port_size::String="24", junction_size::String="0.05",
    outer_ports::Bool=true, anchor_outer_ports::Bool=true,
    graph_attrs::AbstractDict=Dict(), node_attrs::AbstractDict=Dict(),
    edge_attrs::AbstractDict=Dict(), cell_attrs::AbstractDict=Dict())::Graphviz.Graph
  @assert label_attr in (:label, :xlabel, :headlabel, :taillabel)

  # State variables.
  stmts = Graphviz.Statement[]
  port_map = Dict{Port,Graphviz.NodeID}()
  update_port_map! = (v::Int, kind::PortKind, node_ids) -> begin
    for (i, node_id) in enumerate(node_ids)
      port_map[Port(v,kind,i)] = node_id
    end
  end

  # Invisible nodes for incoming and outgoing wires.
  if outer_ports
    gv_box = graphviz_outer_box(f;
      anchor=anchor_outer_ports, orientation=orientation)
    append!(stmts, gv_box.stmts)
    update_port_map!(input_id(f), OutputPort, gv_box.input_ports)
    update_port_map!(output_id(f), InputPort, gv_box.output_ports)
  end

  # Visible nodes for boxes.
  default_attrs = default_directed_graphviz_attrs
  cell_attrs = merge(default_attrs.cell, Graphviz.as_attributes(cell_attrs))
  for v in box_ids(f)
    gv_box = graphviz_box(box(f,v), box_id([v]),
      orientation=orientation, labels=node_labels, port_size=port_size,
      junction_size=junction_size, cell_attrs=cell_attrs)
    append!(stmts, gv_box.stmts)
    update_port_map!(v, InputPort, gv_box.input_ports)
    update_port_map!(v, OutputPort, gv_box.output_ports)
  end

  # Edges.
  for (i, wire) in enumerate(wires(f))
    source, target = wire.source, wire.target
    if !(haskey(port_map, source) && haskey(port_map, target))
      continue
    end
    # Use the port value to label the wire. We take the source port.
    # In most wiring diagrams, the source and target ports should yield the
    # same label, but that is not guaranteed. An advantage of choosing the
    # source port over the target port is that it will carry the
    # "more specific" type when implicit conversions are allowed.
    port = port_value(f, source)
    attrs = Graphviz.Attributes(
      :id => wire_id(Int[], i),
      :comment => edge_label(port),
    )
    if labels
      attrs[label_attr] = edge_label(port)
    end
    edge = Graphviz.Edge(port_map[source], port_map[target]; attrs...)
    push!(stmts, edge)
  end

  # Graph.
  Graphviz.Digraph(graph_name, stmts; prog="dot",
    graph_attrs=merge(default_attrs.graph, Graphviz.as_attributes(graph_attrs),
      Graphviz.Attributes(
        :rankdir => rank_dir(orientation)
      )),
    node_attrs=merge(default_attrs.node, Graphviz.as_attributes(node_attrs)),
    edge_attrs=merge(default_attrs.edge, Graphviz.as_attributes(edge_attrs)))
end

function to_graphviz(f::HomExpr; kw...)::Graphviz.Graph
  to_graphviz(to_wiring_diagram(f); kw...)
end

""" Graphviz box for a generic box.
"""
function graphviz_box(box::AbstractBox, node_id;
    orientation::LayoutOrientation=TopToBottom,
    labels::Bool=true, port_size::String="0",
    cell_attrs::AbstractDict=Dict(), kw...)
  # Main node.
  nin, nout = length(input_ports(box)), length(output_ports(box))
  text_label = labels ? node_label(box.value) : ""
  html_label = box_html_label(nin, nout, text_label;
    orientation=orientation, port_size=port_size, attrs=cell_attrs)
  # Note: The `id` attribute is included in the Graphviz output but is not used
  # internally by Graphviz. It is for use by downstream applications.
  # Reference: http://www.graphviz.org/doc/info/attrs.html#d:id
  node = Graphviz.Node(node_id,
    id = node_id,
    comment = node_label(box.value),
    label = html_label,
  )

  # Input and output ports.
  graphviz_port = (kind::PortKind, port::Int) -> begin
    Graphviz.NodeID(node_id, port_name(kind, port),
                    port_anchor(kind, orientation))
  end
  inputs = [ graphviz_port(InputPort, i) for i in 1:nin ]
  outputs = [ graphviz_port(OutputPort, i) for i in 1:nout ]

  GraphvizBox([node], inputs, outputs)
end

""" Graphviz box for a junction.
"""
function graphviz_box(junction::Junction, node_id;
    junction_size::String="0", kw...)
  graphviz_junction(junction, node_id;
    comment = "junction",
    label = "",
    shape = "circle",
    style = "filled",
    fillcolor = "black",
    width = junction_size,
    height = junction_size,
  )
end

function graphviz_junction(junction::Junction, node_id; kw...)
  node = Graphviz.Node(node_id; id=node_id, kw...)
  nin, nout = length(input_ports(junction)), length(output_ports(junction))
  GraphvizBox([node],
    repeat([Graphviz.NodeID(node_id)], nin),
    repeat([Graphviz.NodeID(node_id)], nout))
end

""" Create "HTML-like" node label for a box.
"""
function box_html_label(nin::Int, nout::Int, text_label::String;
    orientation::LayoutOrientation=TopToBottom, port_size::String="0",
    attrs::AbstractDict=Dict())::Graphviz.Html
  if is_vertical(orientation)
    input_label = ports_horizontal_html_label(InputPort, nin, port_size)
    output_label = ports_horizontal_html_label(OutputPort, nout, port_size)
    Graphviz.Html("""

$(orientation == TopToBottom ? input_label : output_label)

$(escape_html(text_label))

$(orientation == BottomToTop ? input_label : output_label)

""")
  else
    input_label = ports_vertical_html_label(InputPort, nin, port_size)
    output_label = ports_vertical_html_label(OutputPort, nout, port_size)
    Graphviz.Html("""

      $(orientation == LeftToRight ? input_label : output_label)
      $(escape_html(text_label))
      $(orientation == RightToLeft ? input_label : output_label)

""")
  end
end

""" Create horizontal "HTML-like" label for the input or output ports of a box.
"""
function ports_horizontal_html_label(kind::PortKind, nports::Int,
    port_size::String="0")::Graphviz.Html
  cols = if nports > 0
    join("""	"""
         for i in 1:nports)
  else
    """	"""
  end
  Graphviz.Html("""

$cols

""")
end

""" Create vertical "HTML-like" label for the input or output ports of a box.
"""
function ports_vertical_html_label(kind::PortKind, nports::Int,
    port_size::String="0")::Graphviz.Html
  rows = if nports > 0
    join("""

"""
         for i in 1:nports)
  else
    """

"""
  end
  Graphviz.Html("""

$rows
""")
end

""" Graphviz box for the outer box of a wiring diagram.
"""
function graphviz_outer_box(f::WiringDiagram;
    anchor::Bool=true, orientation::LayoutOrientation=TopToBottom)
  # Subgraphs containing invisible nodes.
  stmts = Graphviz.Statement[]
  ninputs, noutputs = length(input_ports(f)), length(output_ports(f))
  if ninputs > 0
    push!(stmts, graphviz_outer_ports(InputPort, ninputs;
      anchor=anchor, orientation=orientation))
  end
  if noutputs > 0
    push!(stmts, graphviz_outer_ports(OutputPort, noutputs;
      anchor=anchor, orientation=orientation))
  end

  # Input and output ports.
  inputs = map(1:ninputs) do i
    Graphviz.NodeID(port_node_name(InputPort, i),
                    port_anchor(OutputPort, orientation))
  end
  outputs = map(1:noutputs) do i
    Graphviz.NodeID(port_node_name(OutputPort, i),
                    port_anchor(InputPort, orientation))
  end

  GraphvizBox(stmts, inputs, outputs)
end

""" Create invisible nodes for the input or output ports of an outer box.
"""
function graphviz_outer_ports(kind::PortKind, nports::Int;
    anchor::Bool=true, orientation::LayoutOrientation=TopToBottom)::Graphviz.Subgraph
  @assert nports > 0
  port_width = "$(round(24/72,digits=3))" # port width in inches
  nodes = [ port_node_name(kind, i) for i in 1:nports ]
  stmts = Graphviz.Statement[
    Graphviz.Node(nodes[i], id=port_name(kind, i)) for i in 1:nports
  ]
  if anchor && length(nodes) >= 2
    push!(stmts, Graphviz.Edge(nodes))
  end
  Graphviz.Subgraph(
    stmts,
    graph_attrs=Graphviz.Attributes(
      :rank => kind == InputPort ? "source" : "sink",
      :rankdir => is_vertical(orientation) ? "LR" : "TB"
    ),
    node_attrs=Graphviz.Attributes(
      :style => "invis",
      :shape => "none",
      :label => "",
      :width => is_vertical(orientation) ? port_width : "0",
      :height => is_horizontal(orientation) ? port_width : "0",
    ),
    edge_attrs=Graphviz.Attributes(
      :style => "invis",
    ),
  )
end
port_node_name(kind::PortKind, port::Int) =
  string(kind == InputPort ? "n0in" : "n0out", port)

""" Graphviz rank direction (`rankdir`) for layout orientation.

Reference: https://www.graphviz.org/doc/info/attrs.html#k:rankdir
"""
rank_dir(orientation::LayoutOrientation) = @match orientation begin
  &TopToBottom => "TB"
  &BottomToTop => "BT"
  &LeftToRight => "LR"
  &RightToLeft => "RL"
end

""" Graphviz anchor for port.
"""
function port_anchor(kind::PortKind, orientation::LayoutOrientation)
  @match (kind, orientation) begin
    (&InputPort, &TopToBottom) => "n"
    (&OutputPort, &TopToBottom) => "s"
    (&InputPort, &BottomToTop) => "s"
    (&OutputPort, &BottomToTop) => "n"
    (&InputPort, &LeftToRight) => "w"
    (&OutputPort, &LeftToRight) => "e"
    (&InputPort, &RightToLeft) => "e"
    (&OutputPort, &RightToLeft) => "w"
  end
end

# Undirected drawing
####################

const default_undirected_graphviz_attrs = (
  graph = Graphviz.Attributes(
    :fontname => default_graphviz_font,
  ),
  node = Graphviz.Attributes(
    :fontname => default_graphviz_font,
    :shape => "ellipse",
    :margin => "0.05,0.025",
    :width => "0.5",
    :height => "0.5",
  ),
  edge = Graphviz.Attributes(
    :fontname => default_graphviz_font,
    :len => "0.5",
  )
)

""" Draw an undirected wiring diagram using Graphviz.

Creates an undirected, bipartite Graphviz graph, with the boxes and outer ports
of the diagram becoming nodes of one kind and the junctions of the diagram
becoming nodes of the second kind.

# Arguments
- `graph_name="G"`: name of Graphviz graph
- `prog="neato"`: Graphviz program, usually "neato" or "fdp"
- `box_labels=false`: if boolean, whether to label boxes with their number;
   if a symbol, name of data attribute for box label
- `port_labels=false`: whether to label ports with their number
- `junction_labels=false`: if boolean, whether to label junctions with their
  number; if a symbol, name of data attribute for junction label
- `junction_size="0.075"`: size of junction nodes, in inches
- `implicit_junctions=false`: whether to represent a junction implicity as a
  wire when it has exactly two incident ports
- `graph_attrs=Dict()`: top-level graph attributes
- `node_attrs=Dict()`: top-level node attributes
- `edge_attrs=Dict()`: top-level edge attributes
"""
function to_graphviz(d::UndirectedWiringDiagram; kw...)::Graphviz.Graph
  to_graphviz(to_graphviz_property_graph(d; kw...))
end

function to_graphviz_property_graph(d::UndirectedWiringDiagram;
    graph_name::String="G", prog::String="neato",
    box_labels::Union{Bool,Symbol}=false, port_labels::Bool=false,
    junction_labels::Union{Bool,Symbol}=false,
    junction_size::String="0.075", implicit_junctions::Bool=false,
    graph_attrs::AbstractDict=Dict(), node_attrs::AbstractDict=Dict(),
    edge_attrs::AbstractDict=Dict())::SymmetricPropertyGraph
  default_attrs = default_undirected_graphviz_attrs
  graph = SymmetricPropertyGraph{Any}(
    name = graph_name, prog = prog,
    graph = merge(default_attrs.graph, Graphviz.as_attributes(graph_attrs)),
    node = merge(default_attrs.node, Graphviz.as_attributes(node_attrs)),
    edge = merge(default_attrs.edge, Graphviz.as_attributes(edge_attrs)),
  )

  # Create nodes for boxes.
  box_vs = add_vertices!(graph, nboxes(d))
  set_vprop!(graph, box_vs, :id, [ "box$b" for b in boxes(d) ])
  labels = if box_labels isa Symbol
    node_label.(subpart(d, box_vs, box_labels))
  else
    box_labels ? string.(boxes(d)) : ""
  end
  set_vprop!(graph, box_vs, :label, labels)

  # Create nodes for outer ports.
  outer_vs = add_vertices!(graph, length(ports(d, outer=true)))
  set_vprop!(graph, outer_vs, :id, [ "outer$p" for p in ports(d, outer=true) ])
  for v in outer_vs
    set_vprops!(graph, v, style="invis", shape="none", label="",
                width="0", height="0", margin="0")
  end

  # Create nodes and edge for junctions.
  for j in junctions(d)
    # Get all ports, including outer ports, incident to junction.
    incident_ports = [
      map(ports_with_junction(d, j)) do port
        b = box(d, port)
        (b, findfirst(p -> p == port, ports(d, b)))
      end;
      map(ports_with_junction(d, j, outer=true)) do port
        (nboxes(d) + port, port)
      end
    ]

    # Case 1: Create edge for junction with exactly two incident ports.
    if implicit_junctions && length(incident_ports) == 2
      (src, src_port), (tgt, tgt_port) = incident_ports
      e, e_inv = add_edge!(graph, src, tgt)
      if port_labels
        set_eprops!(graph, e, taillabel=string(src_port),
                    headlabel=string(tgt_port))
      end
    # Case 2: Create junction node and add edge for every incident port.
    else
      jv = add_vertex!(graph, id="junction$j", comment="junction",
                       shape="circle", label="",
                       style="filled", fillcolor="black",
                       width=junction_size, height=junction_size)
      if junction_labels isa Symbol
        set_vprop!(graph, jv, :xlabel,
                   node_label(subpart(d, j, junction_labels)))
      elseif junction_labels
        set_vprop!(graph, jv, :xlabel, string(j))
      end
      for (v, port) in incident_ports
        e, e_inv = add_edge!(graph, v, jv)
        if port_labels
          set_eprop!(graph, e, :taillabel, string(port))
        end
      end
    end
  end

  return graph
end

# Graphviz labels
#################

""" Create a label for the main content of a box.
"""
node_label(box_value) = box_label(MIME("text/plain"), box_value)

""" Create a label for an edge.
"""
edge_label(port_value) = wire_label(MIME("text/plain"), port_value)

""" Encode attributes for Graphviz HTML-like labels.
"""
function html_attributes(attrs::AbstractDict)::String
  join([ "$(uppercase(string(k)))=\"$v\""
         for (k,v) in Graphviz.as_attributes(attrs) ], " ")
end

""" Escape special HTML characters: &, <, >, ", '

Adapted from [HttpCommon.jl](https://github.com/JuliaWeb/HttpCommon.jl).
"""
function escape_html(s::AbstractString)
  reduce(replace, [
    "&" => "&",
    "\"" => """,
    "'" => "'",
    "<" => "<",
    ">" => ">",
  ], init=s)
end

# Directed layout
#################

""" Lay out directed wiring diagram using Graphviz.

Note: At this time, only the positions and sizes of the boxes, and the positions
of the outer ports, are used. The positions of the box ports and the splines for
the wires are ignored.
"""
function graphviz_layout(diagram::WiringDiagram; kw...)
  graph = to_graphviz(diagram; kw...)
  doc = JSON.parse(Graphviz.run_graphviz(graph, format="json0"))
  graphviz_layout(diagram, parse_graphviz(doc))
end

function graphviz_layout(diagram::WiringDiagram, graph::PropertyGraph)
  # Graphviz uses the standard Cartesian coordinate system (with origin in
  # bottom left corner), while our layout system uses a coordinate system with
  # centered origin and positive y-axis pointing downwards.
  orientation = inverse_rank_dir(get_gprop(graph, :rankdir))
  bounds, pad = get_gprop(graph, :bounds), get_gprop(graph, :pad)
  diagram_size = bounds + 2*pad
  transform_point(p) = SVector(1,-1) .* (p - diagram_size/2)

  # Assumes vertices are in same order as created by `to_graphviz`:
  # 1. Input ports of outer box
  nin, nout = length(input_ports(diagram)), length(output_ports(diagram))
  main_dir = svector(orientation, 1.0, 0.0)
  inputs = map(enumerate(input_ports(diagram))) do (i, value)
    attrs = vprops(graph, i)
    pos = transform_point(attrs[:position])
    PortLayout(; value=value, position=pos, normal=-main_dir)
  end

  # 2. Output ports of outer box
  outputs = map(enumerate(output_ports(diagram))) do (i, value)
    attrs = vprops(graph, nin + i)
    pos = transform_point(attrs[:position])
    PortLayout(; value=value, position=pos, normal=main_dir)
  end

  # 3. Inner boxes
  # TODO: Use port positions from Graphviz layout. Obtain these from either
  # the HTML label or, more likely, an incident edge.
  layout = WiringDiagram(BoxLayout(size=diagram_size), inputs, outputs)
  node_map = Dict{Int,Int}()
  for (i, box) in enumerate(boxes(diagram))
    node = nin + nout + i
    attrs = vprops(graph, node)
    shape = Symbol(get(attrs, :shape, :ellipse))
    if shape == :none; shape = :rectangle end # XXX: Assume HTML label.
    @assert shape == :rectangle "Only the rectangle shape is implemented so far"
    pos, size = transform_point(attrs[:position]), attrs[:size]
    node_map[node] = add_box!(layout, Box(
      BoxLayout(; value=box.value, shape=shape, position=pos, size=size),
      layout_linear_ports(InputPort, input_ports(box), size, orientation),
      layout_linear_ports(OutputPort, output_ports(box), size, orientation)
    ))
  end

  # Add wires using spline points from Graphviz edge layout.
  function map_port(node::Int, portname::Union{String,Nothing})
    if node <= nin
      Port(input_id(layout), OutputPort, node)
    elseif node <= nin + nout
      Port(output_id(layout), InputPort, node - nin)
    else
      # XXX: Parsing the port names is hacky. The alternative is parsing the
      # ordered ports in the HTML node label and matching the port names.
      kind = startswith(portname, "out") ? OutputPort : InputPort
      port = parse(Int, portname[findfirst(r"[0-9]+", portname)])
      Port(node_map[node], kind, port)
    end
  end
  for edge in edges(graph)
    attrs = eprops(graph, edge)
    # Chop off start and end points to obtain intermediate points.
    points = attrs[:spline][3:end-2]
    wire_layout = map(Iterators.partition(points, 3)) do (c1, p, c2)
      WirePoint(transform_point(p),
                normalize(transform_point(c2) - transform_point(c1)))
    end
    add_wire!(layout, Wire(wire_layout,
      map_port(src(graph, edge), get(attrs, :tailport, nothing)),
      map_port(tgt(graph, edge), get(attrs, :headport, nothing))
    ))
  end
  layout
end

""" Lay out ports linearly, equispaced along a rectangular box.

FIXME: Should this be in `WiringDiagramLayouts` as an alternative layout method?
"""
function layout_linear_ports(port_kind::PortKind, port_values::Vector,
                             box_size::StaticVector{2}, orientation::LayoutOrientation)
  n = length(port_values)
  sign = port_sign(port_kind, orientation)
  normal = svector(orientation, sign, 0.0)
  map(port_values, range(-1,1,length=n+2)[2:end-1]) do value, coeff
    pos = svector(orientation, sign, coeff) .* box_size/2
    PortLayout(; value=value, position=pos, normal=normal)
  end
end

""" Layout orientation for Graphviz rank direction (`rankdir`).

Reference: https://www.graphviz.org/doc/info/attrs.html#k:rankdir
"""
inverse_rank_dir(rank_dir::String) = @match rank_dir begin
  "TB" => TopToBottom
  "BT" => BottomToTop
  "LR" => LeftToRight
  "RL" => RightToLeft
end

end


text:
using BinDeps

@BinDeps.setup
group = library_group("SOFA")

libsofa_c = library_dependency("libsofa_c", aliases = ["libsofa_c", "libsofa"], runtime = true, group = group)

prefix = joinpath(BinDeps.depsdir(libsofa_c), "usr")
libdir = joinpath(BinDeps.depsdir(libsofa_c), "usr", "lib")
srcdir = joinpath(BinDeps.depsdir(libsofa_c), "src", "sofa")

provides(SimpleBuild,
    (@build_steps begin
    CreateDirectory(prefix)
    CreateDirectory(libdir)
        @build_steps begin
            ChangeDirectory(srcdir)
            `make clean`
            `make`
            `cp libsofa_c.so libsofa_c.dylib`
            `mv libsofa_c.so $prefix/lib`
            `mv libsofa_c.dylib $prefix/lib`
            `make clean`
        end
    end), libsofa_c, os = :Unix)

@BinDeps.install Dict(:libsofa_c => :libsofa_c)


text:
"""
    struct LambdaComponents
        function_name::String
        aws_config::AWSConfig
        local_image::Union{Nothing, LocalImage}
        remote_image::Union{Nothing, RemoteImage}
        lambda_function::Union{Nothing, LambdaFunction}
    end
"""
@with_kw mutable struct LambdaComponents
  function_name::String
  aws_config::AWSConfig
  local_image::Union{Nothing, LocalImage}
  remote_image::Union{Nothing, RemoteImage}
  lambda_function::Union{Nothing, LambdaFunction}
end
Base.show(l::LambdaComponents) = "$(l.local_image)\t$(l.remote_image)\t$(l.lambda_function)"

function get_from_any_component(
    get_func::Function,
    l::LambdaComponents,
  )
  for (f_name, f_type) in zip(fieldnames(LambdaComponents), fieldtypes(LambdaComponents))
    if f_type <: LambdaComponent
      try
        out = get_func(getfield(l, f_name))
        return out
      catch e
        continue
      end
    end
  end
end

"""
    create_lambda_components(
        res::AbstractResponder;
        image_suffix::Union{Nothing, String} = nothing,
        aws_config::Union{Nothing, AWSConfig} = nothing,
        image_tag::String = "latest",
        no_cache::Bool = false,
        julia_base_version::String = "1.6.1",
        julia_cpu_target::String = "x86-64",
        package_compile::Bool = false,
        user_defined_labels::AbstractDict{String, String} = OrderedDict{String, String}(),
      )::LambdaComponents

Creates a `LocalImage` from the given responder, creates a `LambdaComponents` object to store the
local image, and then returns the `LambdaComponents` object.

Acts as an alternative to `create_local_image`, but returns a `LambdaComponents` rather than just
the local image. This can be more convenient for keeping the components of a lambda function
organized - for example:
`create_lambda_components(responder) |> with_remote_image! |> with_lambda_function!` will run through
the entire process of creating a local image, pushing that image to ECR, and then creating a Lambda
function.
"""
function create_lambda_components(
    res::AbstractResponder;
    image_suffix::Union{Nothing, String} = nothing,
    aws_config::Union{Nothing, AWSConfig} = nothing,
    image_tag::String = "latest",
    no_cache::Bool = false,
    julia_base_version::String = "1.6.3",
    julia_cpu_target::String = "x86-64",
    package_compile::Bool = false,
    user_defined_labels::AbstractDict{String, String} = OrderedDict{String, String}(),
  )::LambdaComponents
  local_image = create_local_image(res;
                                   image_suffix = image_suffix,
                                   aws_config = aws_config,
                                   image_tag = image_tag,
                                   no_cache = no_cache,
                                   julia_base_version = julia_base_version,
                                   julia_cpu_target = julia_cpu_target,
                                   package_compile = package_compile,
                                   user_defined_labels = user_defined_labels)
  LambdaComponents(get_response_function_name(local_image),
                   get_aws_config(),
                   local_image,
                   nothing,
                   nothing,
  )
end

"""
    with_remote_image!(l::LambdaComponents)::LambdaComponents

Adds a 'RemoteImage` object to the passed `LambdaComponents` instance. Will error if the instance
has neither an existing remote image or local image.
"""
function with_remote_image!(l::LambdaComponents)::LambdaComponents
  if isnothing(l.local_image) && isnothing(l.remote_image)
    error("Unable to add remote image to LambdaComponents as it does not have a local image")
  end
  if isnothing(l.remote_image)
    l.remote_image = push_to_ecr!(l.local_image)
    l
  else
    l
  end
end

"""
    with_lambda_function!(l::LambdaComponents)::LambdaComponents

Adds a 'LambdaFunction` instance to the passed `LambdaComponents` instance. Will error if the instance
has neither an existing remote image, local image or lambda function.
"""
function with_lambda_function!(l::LambdaComponents)::LambdaComponents
  with_remote = if isnothing(l.remote_image) && isnothing(l.lambda_function)
    try
      with_remote_image!(l)
    catch e
      error("Unable to create lambda function from LambdaComponents; it has neither a local image nor a remote image")
    end
  else
    l
  end
  if isnothing(with_remote.lambda_function)
    l.lambda_function = create_lambda_function(with_remote.remote_image)
    l
  else
    l
  end
end

"""
    delete!(l::LambdaComponents)

Deletes the local docker image, remote image, and lambda function associated with the
`LambdaComponents` instance.
"""
function delete!(l::LambdaComponents)
  !isnothing(l.lambda_function) && delete!(l.lambda_function)
  !isnothing(l.remote_image) && delete!(l.remote_image)
  !isnothing(l.local_image) && delete!(l.local_image)
end

"""
    run_test(
        l::LambdaComponents;
        function_argument::Any = "",
        expected_response::Any = nothing;
      )::Tuple{Bool, Float64}

Tests the passed `LambdaComponents` instance.

The test runs on the most downstream object. So if the instance has a `LambdaFunction`, this will
be tested. Otherwise, the attached `LocalImage` will be tested. If the `LambdaComponents` object
has neither a local image or a lambda function, then it has nothing that can be tested and the
function will throw an error.

Returns a tuple of {Test pass/fail, Test time taken in seconds}.
"""
function run_test(
    l::LambdaComponents,
    function_argument::Any = "",
    expected_response::Any = nothing,
  )::Tuple{Bool, Union{Missing, Float64}}
  if !isnothing(l.lambda_function)
    run_test(l.lambda_function, function_argument, expected_response)
  elseif !isnothing(l.local_image)
    run_test(l.local_image, function_argument, expected_response)
  else
    error("Unable to test LambdaComponents object; it has neither a local image or a lambda function")
  end
end

function matches(res::AbstractResponder, local_image::LocalImage)::Bool
  tree_hash = get_tree_hash(local_image)
  isnothing(tree_hash) ? false : (get_tree_hash(res) == tree_hash)
end
matches(local_image::LocalImage, res::AbstractResponder) = matches(res, local_image)

function matches(local_image::LocalImage, ecr_repo::ECRRepo)::Bool
  local_image.Repository == ecr_repo.repositoryUri
end
matches(ecr_repo::ECRRepo, local_image::LocalImage) = matches(local_image, ecr_repo)

function matches(local_image::LocalImage, remote_image::RemoteImage)::Bool
  @debug local_image.Digest
  @debug remote_image.imageDigest
  local_image.Digest == remote_image.imageDigest
end
matches(remote_image::RemoteImage, local_image::LocalImage) = matches(local_image, remote_image)

function matches(res::AbstractResponder, remote_image::RemoteImage)::Bool
  tree_hash = get_tree_hash(remote_image)
  isnothing(tree_hash) ? false : get_tree_hash(res) == tree_hash
end
matches(remote_image::RemoteImage, res::AbstractResponder) = matches(res, remote_image)

function matches(remote_image::RemoteImage, lambda_function::LambdaFunction)::Bool
  hash_only = split(remote_image.imageDigest, ':')[2]
  hash_only == lambda_function.CodeSha256
end
matches(lambda_function::LambdaFunction, remote_image::RemoteImage) = matches(remote_image, lambda_function)

function combine_if_matches(l1::LambdaComponents, l2::LambdaComponents)::Union{Nothing, LambdaComponents}
  get_non_nothing_type(x::Type{IT}) where {IT} = typeof(x) == Union ? x.b : x

  lambda_types = map(get_non_nothing_type, fieldtypes(LambdaComponents))
  lambda_names = fieldnames(LambdaComponents)
  l1_fields = [getfield(l1, name) for name in lambda_names]
  l2_fields = [getfield(l2, name) for name in lambda_names]

  function match_across_fields(
      l1_flds::Vector,
      l2_flds::Vector,
    )::Bool
    for (fieldtype_1, val_1) in zip(lambda_types, l1_flds)
      for (fieldtype_2, val_2) in zip(lambda_types, l2_flds)
        if !isnothing(val_1) && !isnothing(val_2) && hasmethod(matches, Tuple{fieldtype_1, fieldtype_2})
          if matches(val_1, val_2)
            return true
          end
        end
      end
    end
    return false
  end

  function cmb(c1::Union{Nothing, T}, c2::Union{Nothing, T})::Union{Nothing, T} where {T}
    if isnothing(c1)
      c2
    elseif isnothing(c2)
      c1
    else
      c1 != c2 && error("Found non-matching element in matched lambda")
      c1
    end
  end

  function cmb(c1::T, c2::T)::T where {T}
    c1 != c2 && error("Found non-matching element in matched lambda")
    c1
  end

  if match_across_fields(l1_fields, l2_fields)
    cmb_fields = Dict(sym => cmb(f_1, f_2) for (sym, f_1, f_2) in zip(lambda_names, l1_fields, l2_fields))
    LambdaComponents(; cmb_fields...)
  else
    nothing
  end
end

struct TableComponent
  name::String
  value_function::Function
  highlighter_f::Union{Function, Nothing}
end

const not_present = "-"

function_name_f(l::LambdaComponents)::String = l.function_name
const function_name_component = TableComponent("Function Name", function_name_f, nothing)

function account_id_f(l::LambdaComponents)::String
  l.aws_config.account_id
end
const account_id_component = TableComponent("Account ID", account_id_f, nothing)

function responder_source_f(l::LambdaComponents)::String
  src = get_labels(l).RESPONDER_PKG_SOURCE
  isnothing(src) ? not_present : src
end
function responder_source_h_f(
    headers::OrderedDict{String, TableComponent}, lambdas::Vector{LambdaComponents},
  )::Vector{Highlighter}
  h1 = Highlighter(bold = true, foreground = :blue) do table_data, i, j
    comps = values(headers) |> collect
    if comps[j] == responder_source_component
      responder_path = table_data[i, j]
      lc_tree_hash = get_tree_hash(lambdas[i])
      if ispath(responder_path)
        get_tree_hash(dirname(responder_path)) == lc_tree_hash
      else
        false
      end
    else
      false
    end
  end

  h2 = Highlighter(foreground = :dark_gray) do table_data, i, j
    comps = values(headers) |> collect
    if comps[j] == responder_source_component
      responder_path = table_data[i, j]
      ispath(responder_path) ? false : true
    else
      false
    end
  end
  [h1, h2]
end
const responder_source_component = TableComponent("Responder Source", responder_source_f, responder_source_h_f)

function tree_hash_f(l::LambdaComponents)::String
  isnothing(l.local_image) && return not_present
  hsh = get_tree_hash(l.local_image)
  isnothing(hsh) ? not_present : hsh[1:docker_hash_limit]
end
const tree_hash_component = TableComponent("Tree Hash", tree_hash_f, nothing)

local_image_name_f(l::LambdaComponents)::String = isnothing(l.local_image) ? not_present : get_lambda_name(l.local_image)
const local_image_name_component = TableComponent("Image Name", local_image_name_f, nothing)

local_image_id_f(l::LambdaComponents)::String = isnothing(l.local_image) ? not_present : l.local_image.ID
const local_image_id_component = TableComponent("Image ID", local_image_id_f, nothing)

local_image_tag_f(l::LambdaComponents)::String = isnothing(l.local_image) ? not_present : l.local_image.Tag
const local_image_tag_component = TableComponent("Image Tag", local_image_tag_f, nothing)

function remote_image_tag_f(l::LambdaComponents)::String
  isnothing(l.remote_image) && return not_present
  itag = l.remote_image.imageTag
  ismissing(itag) ? not_present : itag
end
const remote_image_tag_component = TableComponent("Image Tag", remote_image_tag_f, nothing)

function remote_image_digest_f(l::LambdaComponents)::String
  isnothing(l.remote_image) && return "-"
  digest = l.remote_image.imageDigest
  if ismissing(digest)
    "-"
  else
    hash_only = split(digest, ':') |> last
    hash_only[begin:docker_hash_limit]
  end
end
const remote_image_digest_component = TableComponent("Image Digest", remote_image_digest_f, nothing)

function lambda_function_name_f(l::LambdaComponents)::String
  isnothing(l.lambda_function) && return not_present
  fn = l.lambda_function.FunctionName
  ismissing(fn) ? not_present : fn
end
const lambda_function_name_component = TableComponent("Function Name", lambda_function_name_f, nothing)

function lambda_function_last_modified_f(l::LambdaComponents)::String
  isnothing(l.lambda_function) && return not_present
  lm = l.lambda_function.LastModified
  ismissing(lm) ? not_present : lm
end
const lambda_function_last_modified_component = TableComponent(
  "Last Modified", lambda_function_last_modified_f, nothing
)

function get_table_data(
    headers::OrderedDict{String, TableComponent},
    lambdas::Vector{LambdaComponents},
  )::Matrix{String}

  # TODO refactor common interface for all getters - check component then check sub
  # TODO highlight responder path if not current - probably put it in grey
  all_funcs = values(headers) |> collect
  data_rows = [map(tc -> tc.value_function(l), all_funcs) for l in lambdas]
  data_rows = map(row-> reshape(row, (1, :)), data_rows)
  data = vcat(data_rows...)
end

"""
    show_lambdas()::Nothing

Displays a table of all objects generated using Jot.jl.

Each row of the table shows at least one of a Responder, a local docker image, a remote (hosted
on AWS ECR) docker image, and an AWS-hosted Lambda function. The local docker image, remote docker
image and lambda function on a given row of the table are guaranteed to share the same underlying
function code.

The Responder column is colour-coded:
- Grey indicates that this path no longer exists.
- White indicates that this path still exists, but the code has changed since the objects shown
in the row were created.
- Blue indicates that the underlying code for this row (eg the code present in the local image,
remote image etc) is the same as is currently present at this path.
"""
function show_lambdas()
  @info "Collecting lambda components; this may take a few seconds..."
  lambdas = get_all_lambdas()

  table_headers = OrderedDict(
    "Function Name" => function_name_component,
    "Responder" => responder_source_component,
    "Local Image" => local_image_id_component,
    "Remote Image" => remote_image_digest_component,
    "Lambda Function" => lambda_function_name_component,
  )

  table_data = get_table_data(table_headers, lambdas)
  headers_matrix = (keys(table_headers) |> collect, [tc.name for tc in values(table_headers)])

  table_comps = filter(tc -> !isnothing(tc.highlighter_f), values(table_headers) |> collect)
  highlighter_funcs = map(tc -> tc.highlighter_f, table_comps)
  highlighters = [h for hf in highlighter_funcs for h in hf(table_headers, lambdas)]
  highlighters = tuple(highlighters...)

  pretty_table(
    table_data;
    header=headers_matrix,
    show_row_number=true,
    crop=:none,
    maximum_columns_width=30,
    highlighters=highlighters,
  )
end

function get_all_lambdas()::Vector{LambdaComponents}
  all_local = get_all_local_images()
  all_remote = get_all_remote_images()
  all_functions = get_all_lambda_functions()
  aws_config = get_aws_config()
  local_lambdas = [
    LambdaComponents(get_labels(l) |> get_responder_full_function_name, aws_config, l, nothing, nothing)
    for l in all_local if (is_lambda(l) && is_jot_generated(l))
  ]
  remote_lambdas = [
    LambdaComponents(get_labels(r) |> get_responder_full_function_name, aws_config, nothing, r, nothing)
    for r in all_remote if is_jot_generated(r)
  ]
  func_lambdas = [
    LambdaComponents(get_labels(f) |> get_responder_full_function_name, aws_config, nothing, nothing, f)
    for f in all_functions if is_jot_generated(f)
  ]
  all_lambdas = [ local_lambdas ; remote_lambdas ; func_lambdas ]

  function match_off_lambdas(
      to_match::Vector{LambdaComponents},
      matched::Vector{LambdaComponents}
  )::Vector{LambdaComponents}
    if length(to_match) == 0
      matched
    else
      match_head = to_match[1]; match_tail = to_match[2:end]
      add_to_matched = match_head
      for (i, m) in enumerate(matched)
        cmb = combine_if_matches(match_head, m)
        if !isnothing(cmb)
          add_to_matched = cmb
          deleteat!(matched, i)
          break
        end
      end
      match_off_lambdas(match_tail, [matched; [add_to_matched]])
    end
  end
  match_off_lambdas(all_lambdas, Vector{LambdaComponents}())
end

function group_by_function_name(lambdas::Vector{LambdaComponents})::Dict{String, Vector{LambdaComponents}}
  has_local_image = filter(l -> !isnothing(l.local_image), lambdas)
  func_names = map(l -> get_response_function_name(l.local_image), has_local_image)
  lambdas_by_function = Dict()
  for (func_name, lambda) in zip(func_names, has_local_image)
    if !isnothing(lambda.local_image)
      if !isnothing(func_name)
        lambdas_for_name = get(lambdas_by_function, mod_func_name, Vector{LambdaComponents}())
        lambdas_by_function[mod_func_name] = [lambdas_for_name ; [lambda]]
      end
    end
  end
  lambdas_by_function
end


text:
function test_fba(model, obj_ider, cost_ider = nothing; summary = false)
    hascost = cost_ider in model.rxns
    fbaout = !hascost ? ChLP.fba(model, obj_ider) : ChLP.fba(model, obj_ider, cost_ider)

    # info
    obj_val = ChU.av(model, fbaout, obj_ider)
    exp_obj_val = maximum(NanchenData.val("D"))
    cost_ider = !hascost ? "Not included" : cost_ider
    cost_val = !hascost ? "Not included" : ChU.av(model, fbaout, cost_ider)
    @info("FBA SOLUTION", obj_val, exp_obj_val, cost_ider, cost_val)

    summary && ChU.summary(model, fbaout)

    nothing
end

function test_fba(exp::Int, model, obj_ider, cost_ider = nothing; summary = false)
    hascost = cost_ider in model.rxns
    fbaout = !hascost ? ChLP.fba(model, obj_ider) : ChLP.fba(model, obj_ider, cost_ider)

    # info
    obj_val = ChU.av(model, fbaout, obj_ider)
    exp_obj_val = NanchenData.val("D", exp)
    cost_ider = !hascost ? "Not included" : cost_ider
    cost_val = !hascost ? "Not included" : ChU.av(model, fbaout, cost_ider)
    @info("FBA SOLUTION", obj_val, exp_obj_val, cost_ider, cost_val)

    summary && ChU.summary(model, fbaout)
    nothing
end


text:
# Paulis
@test X * Id == X
@test Id * X == X
@test Y * Id == Y
@test Z * Id == Z

@test X * X == Id
@test Y * Y == Id
@test Z * Z == Id

@test X * Y == im∘Z
@test Y * Z == im∘X
@test Z * X == im∘Y
@test Y * X == -im∘Z

@test Id < X
@test X < Y
@test Y < Z

@test II < IX
@test IX < XX
@test XI < XX
@test IX < IY
@test XI < YI
@test XX < YY
@test XX < XY
@test YX > XX

# generators
@test generators(X) == [X]
@test generators(Z) == [Z]
@test generators(Y) == [im∘X, Z]

@test generators(XX) == [XI, IX]
@test generators(ZZ) == [ZI, IZ]
@test generators(XY) == [im∘XI, IX, IZ]
@test generators(YX) == [im∘XI, ZI, IX]
@test generators(YY) == [-XI, ZI, IX, IZ]

XII = kron(X,Id,Id)
IXI = kron(Id,X,Id)
IIX = kron(Id,Id,X)
ZII = kron(Z,Id,Id)
IZI = kron(Id,Z,Id)
IIZ = kron(Id,Id,Z)

@test generators(kron(X,X,X)) == [XII, IXI, IIX]
@test generators(kron(Z,Z,Z)) == [ZII, IZI, IIZ]
@test generators(kron(X,Y,Z)) == [im∘XII, IXI, IZI, IIZ]

@test Pauli([1 0; 0 1]) == Id
@test Pauli([0 1; 1 0]) == X
@test Pauli([0 -im; im 0]) == Y
@test Pauli([1 0; 0 -1]) == Z

@test complex(paulieye(2)) == eye(4)

# Factors
@test factor(XII)  == [ X,Id,Id]
@test factor(XII*IZI)  == [ X,Z,Id]
@test factor(-XII)  == [-X,Id,Id]
@test factor(-XII*IZI)  == [-X,Z,Id]
@test factor(-IXI)  == [-Id,X,Id]
@test factor(-IXI*IIZ)  == [-Id,X,Z]
@test factor(im∘XII)  == [im∘X,Id,Id]
@test factor(im∘XII*IZI)  == [im∘X,Z,Id]
@test factor(im∘IXI)  == [im∘Id,X,Id]
@test factor(im∘IXI*IIZ)  == [im∘Id,X,Z]

# Weight
for i=1:20
  v = rand(10)
  vX = (v .> .25) .* (v .< .5)
  vY = 2*(v .> .5) .* (v .< .75)
  vZ = 3*(v .> .75)
  r = rand()
  rS = (r > .25) + (r > .5) + (r > .75)
  @test weight(Pauli(vX + vY + vZ,rS)) == sum(v .> .25)
end

# Matrix conversions
@test complex(X) == [0 1; 1 0]
@test complex(Y) == [0 -1im; 1im 0]
@test complex(Z) == [1 0; 0 -1]

# Matrix arthimetic
@test iszero(X - [0 1; 1 0])
@test iszero(Y - [0 -1im;1im 0])
@test iszero(Z - [1 0; 0 -1])

# weight returns Int
@test isa(weight(Pauli(1)), Int)

@testset "Construct from matrix" begin
    xmat = Complex{Float64}[4.702883518586134e-7 + 0.0im 0.9999991609182552 + 0.0im; 0.9999990516115099 + 0.0im 5.793263216649176e-7 + 0.0im]
    @test Pauli(xmat) == X
    ymat = -im * Complex{Float64}[7.6121754928998e-7 + 0.0im 5.577787991615943e-7 - 0.9999984033382637im; 2.066866029117562e-7 + 0.9999984033382637im 8.745277847991395e-8 + 0.0im]
    @test Pauli(ymat) == -im ∘ Y
end


text:
using Genie.Router
using Genie.Renderer.Json
using GRIB

grib_file = "/home/neo/Projekte/Predly/weather/gfs_workflow/download/data/gfs.20210109/06/gfs.t06z.pgrb2.0p25.f002"
f = GribFile(grib_file)

function getVals(band, lat, lon)
  for (index, message) in enumerate(f)

    start_time = time()
    data_array = []
    parameter_name = message["parameterName"]
    parameter_unit = message["parameterUnits"]
    parameter_ecmf = message["nameECMF"]
    units_ecmf     = message["unitsECMF"]
    name           = message["name"]
    units          = message["units"]
    level          = message["level"]
    type_of_level  = message["typeOfLevel"]
    short_name     = message["shortName"]

      if index == band

          near = Nearest(message)
          lons, lats, values, distances = find(near, message, lat, lon)
          seek(f, 0)
          println("Took ", time()-start_time)
          println(lons, lats, values, distances)
          destroy(near)
          return [lons,lats,values,distances]
      end
  end
end

route("/") do
  serve_static_file("welcome.html")
end

function plus(a, b)
  return string(a+b , a , b)
end

route("/hello") do

  band = haskey(@params, :band) ? @params(:band) : 450
  lat = haskey(@params, :lat) ? @params(:lat) : 50.0
  lon = haskey(@params, :lon) ? @params(:lon) : 8.25

  start_time = time()
  resp = Dict()

  println( typeof(band), typeof(lat), typeof(lon))

  getVals(band, lat, lon)

  """
  for n in 1:length(lats)
    resp[n] = Dict()
    resp[n]["lat"] = lats[n]
    resp[n]["lon"] = lons[n]
    resp[n]["val"] = vals[n]
  end
  """

  return( "hi")

end


text:
#=
surface_test:
- Julia version: 1.1.0
- Author: Jirik
=#
using Test
using Logging
using SparseArrays
using Pio3d
using Distributed
if nprocs() == 1
    addprocs(3)
end
# # using Revise
# ENV["JULIA_DEBUG"] = "surface_extraction_parallel"
# ENV["JULIA_DEBUG"] = "all"
# using ViewerGL
# using Plasm
# using LarSurf
@everywhere using LarSurf
# @everywhere using Distributed
# global_logger(SimpleLogger(stdout, Logging.Debug))
# # set logger on all workers
# for wid in workers()
#     @spawnat wid global_logger(SimpleLogger(stdout, Logging.Debug))
# end
@testset "Init and deinit" begin
    block_size = [2, 2, 2]
    pth = Pio3d.datasets_join_path("medical/orig/sample_data/nrn4.pklz")
	datap = Pio3d.read3d(pth)
	data3d = datap["data3d"]
	voxelsize_mm = datap["voxelsize_mm"]
    @info "File path " pth
    data = LarSurf.Experiments.report_init_row(@__FILE__)
    V1, FVtri = LarSurf.Experiments.experiment_make_surf_extraction_and_smoothing(
    data3d, voxelsize_mm;
    data=data,
	output_csv_file="test_nrn4.csv"
    )
    @test length(V1) > 1
    @info "size V1 = " size(V1)
    @info "size FVtri = " size(FVtri)
    @test size(V1, 1) == 3
    @test size(V1, 2) > 10
    @test size(FVtri[1], 1) == 3
    @test size(FVtri, 1) > 10

end