text:
# # Rate Constants and T-Dependency
#
# Here, we will just briefly summarize how temperature affects enzymatic rate constants and how much uncertainty there actually is in literature (in terms of how to best define them and how plastic some of the variables are)

## Add usual tools we use:
using BenchmarkTools
using PyPlot

#----------------------------------------------------------------------------

## load Photosynthesis module:
using Land.Photosynthesis
#----------------------------------------------------------------------------

## Specify Field Type
const FT = Float32;
#----------------------------------------------------------------------------

# ## T-dependence
# In literature, there are different implementations of the temperature dependence of enzymatic reaction rates. Common among all of them is a typical Arrhenius formulation with Activation energy $E_a$, which leads to a temperature dependence of a quantity $V$ given the standard value defined at 25 °C as $T_{ref}$.
#
# $$V(T)= V(T_{ref})\, \underbrace{\exp\left(\frac{E_a}{RT_{ref}}-\frac{E_a}{RT}\right)}_{\text{Activation}}$$
#
# Other formulations add de-activation of proteins due to denaturalization at higher temperatures:
#
# $$V(T)= V(T_{ref})\, \exp\left(\frac{E_a}{RT_{ref}}-\frac{E_a}{RT}\right) \underbrace{\frac{1+\exp\left[(\Delta ST_{ref}-H_d)/RT_{ref}\right]}{1+\exp\left[(\Delta ST-H_d)/RT\right]}}_{\text{de-activation}}$$
#
# which includes an entropy term $\Delta S$ and the energy for deactivation $H_d$.
#
# To illustrate the differences, we show two different implementations, one only using the activation part based on Bernacchi et al 2001[^1] and one with the typical CLM5 implementation.
#
# The T-Dependency applies to $V_{max}$, $J_{max}$, $K_{c}$, $K_{o}$, $K_{pep}$, Respiration, and $\Gamma^{*}$, and there are multiple parameter sets for each of the parameters. Below, we show the differences of the sets for each of $V_{max}$, $J_{max}$, $K_{c}$, $K_{o}$, $K_{pep}$, Respiration, and $\Gamma^{*}$.
#
# [^1]: Bernacchi, C.J., Singsaas, E.L., Pimentel, C., Portis Jr, A.R. and Long, S.P., 2001. Improved temperature response functions for models of Rubisco‐limited photosynthesis. Plant, Cell & Environment, 24(2), pp.253-259.
#
# ### $V_{cmax}$
#

## Bernacchi set
vc_td_bernacchi = Photosynthesis.VcmaxTDBernacchi(FT)
## CLM set
vc_td_clm       = Photosynthesis.VcmaxTDCLM(FT)
## Leuning set
vc_td_leuning   = Photosynthesis.VcmaxTDLeuning(FT)

## an array of temperature from 10 to 50 degree C
t_array = collect(FT, 10:50) .+ FT(273.15)
vcm_bernacchi = photo_TD_from_val(vc_td_bernacchi, FT(60), t_array)
vcm_clm       = photo_TD_from_val(vc_td_clm      , FT(60), t_array)
vcm_leuning   = photo_TD_from_val(vc_td_leuning  , FT(60), t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, vcm_bernacchi, "r-", label="Bernacchi")
plot(t_array, vcm_clm      , "g-", label="CLM5"     )
plot(t_array, vcm_leuning  , "b-", label="Leuning"  )
xlabel("Temperature (K)"            , fontsize=16)
ylabel("\$V_{cmax}\$ (μmol m⁻² s⁻¹)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $V_{pmax}$
#

## Boyd set
vp_td_boyd = Photosynthesis.VpmaxTDBoyd(FT)

## an array of temperature from 10 to 50 degree C
vpm_boyd = photo_TD_from_val(vp_td_boyd, FT(60), t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, vpm_boyd, "r-", label="Boyd")
xlabel("Temperature (K)"            , fontsize=16)
ylabel("\$V_{pmax}\$ (μmol m⁻² s⁻¹)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $J_{max}$
#

## Bernacchi set
j_td_bernacchi = Photosynthesis.JmaxTDBernacchi(FT)
## CLM set
j_td_clm       = Photosynthesis.JmaxTDCLM(FT)
## Leuning set
j_td_leuning   = Photosynthesis.JmaxTDLeuning(FT)

## an array of temperature from 10 to 50 degree C
jm_bernacchi = photo_TD_from_val(j_td_bernacchi, FT(120), t_array)
jm_clm       = photo_TD_from_val(j_td_clm      , FT(120),  t_array)
jm_leuning   = photo_TD_from_val(j_td_leuning  , FT(120),  t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, jm_bernacchi, "r-", label="Bernacchi")
plot(t_array, jm_clm      , "g-", label="CLM5"     )
plot(t_array, jm_leuning  , "b-", label="Leuning"  )
xlabel("Temperature (K)"           , fontsize=16)
ylabel("\$J_{max}\$ (μmol m⁻² s⁻¹)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $K_{c}$
#

## Bernacchi set
kc_td_bernacchi = Photosynthesis.KcTDBernacchi(FT)
## CLM set
kc_td_clm       = Photosynthesis.KcTDCLM(FT)

## an array of temperature from 10 to 50 degree C
kc_bernacchi = photo_TD_from_set(kc_td_bernacchi, t_array)
kc_clm       = photo_TD_from_set(kc_td_clm      , t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, kc_bernacchi, "r-", label="Bernacchi")
plot(t_array, kc_clm      , "g-", label="CLM5"     )
xlabel("Temperature (K)", fontsize=16)
ylabel("\$K_{c}\$ (Pa)" , fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $K_{o}$
#

## Bernacchi set
ko_td_bernacchi = Photosynthesis.KoTDBernacchi(FT)
## CLM set
ko_td_clm       = Photosynthesis.KoTDCLM(FT)

## an array of temperature from 10 to 50 degree C
ko_bernacchi = photo_TD_from_set(ko_td_bernacchi, t_array)
ko_clm       = photo_TD_from_set(ko_td_clm      , t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, ko_bernacchi, "r-", label="Bernacchi")
plot(t_array, ko_clm      , "g-", label="CLM5"     )
xlabel("Temperature (K)", fontsize=16)
ylabel("\$K_{o}\$ (Pa)" , fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $K_{pep}$
#

## Boyd set
kpep_td_boyd = Photosynthesis.KpepTDBoyd(FT)
## CLM set
kpep_td_clm  = Photosynthesis.KpepTDCLM(FT)

## an array of temperature from 10 to 50 degree C
kpep_boyd = photo_TD_from_set(kpep_td_boyd, t_array)
kpep_clm  = photo_TD_from_set(kpep_td_clm , t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, kpep_boyd, "r-", label="Boyd")
plot(t_array, kpep_clm , "g-", label="CLM5")
xlabel("Temperature (K)" , fontsize=16)
ylabel("\$K_{pep}\$ (Pa)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### Respiration
#

## Bernacchi set
re_td_bernacchi = Photosynthesis.RespirationTDBernacchi(FT)
## CLM set
re_td_clm       = Photosynthesis.RespirationTDCLM(FT)

## an array of temperature from 10 to 50 degree C
re_bernacchi = photo_TD_from_val(re_td_bernacchi, FT(1.0), t_array)
re_clm       = photo_TD_from_val(re_td_clm      , FT(1.0), t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, re_bernacchi, "r-", label="Bernacchi")
plot(t_array, re_clm      , "g-", label="CLM5"     )
xlabel("Temperature (K)"                , fontsize=16)
ylabel("Respiration Rate (μmol m⁻² s⁻¹)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ### $\Gamma^{*}$
#

## Bernacchi set
Γs_td_bernacchi = Photosynthesis.ΓStarTDBernacchi(FT)
## CLM set
Γs_td_clm       = Photosynthesis.ΓStarTDCLM(FT)

## an array of temperature from 10 to 50 degree C
Γs_bernacchi = photo_TD_from_set(Γs_td_bernacchi, t_array)
Γs_clm       = photo_TD_from_set(Γs_td_clm      , t_array)

## plot the TD curves
figure(figsize=(8,6), dpi=100)
tight_layout(true)
plot(t_array, Γs_bernacchi, "r-", label="Bernacchi")
plot(t_array, Γs_clm      , "g-", label="CLM5"     )
xlabel("Temperature (K)"     , fontsize=16)
ylabel("\$\\Gamma^{*}\$ (Pa)", fontsize=16)
legend(loc="upper left")
gcf()
#----------------------------------------------------------------------------

# ## Summary
# The way land surface models implement temperature variations in $V_{c,max}$ is still rather variable and there is no clear consensus in what model is best. In addition, gas exchange datasets performed at different temperatures might actually alias different confounding factors into the $V_{c,max}$ determination, so that the actually fitted $V_{c,max}$ strongly depends on the model formulation. E.g. mesophyll conductance is usually ignored in deriving $V_{c,max}$ but might be temperature dependent, which can cause T-dependent errors in $V_{c,max}$. Also, the T-dependencies may change with time, known as acclimation of optimal temperature.
#
# That said, having the option to use different model forumations (as well as parameter settings) and actually using the ones that are implemented in land surface models directly when fitting leaf level measurements seems to be a prudent way to improve models in the future. More leaf level data might help in the future, especially looking at temperature dependencies, mesophyll conductance, and fluorescence yields. So many chemical and phsyical aspects change with temperature that it is hard to isolate impacting factors. Now, at least, we can run the models directly in the REPL, which is great for prototyping anc validating the model on the individual module level without any effort in terms of additional coding.

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# Generate completions.json

open(joinpath(dirname(@__FILE__), "completions.json"), "w") do io
  println(io, "{")
  for (word, char) in Base.REPLCompletions.latex_symbols
    println(io, "  \"$(word[2:end])\": \"$char\",")
  end
  for (word, char) in Base.REPLCompletions.emoji_symbols
    println(io, "  \"$(word[2:end])\": \"$char\",")
  end
  skip(io, -2)
  println(io)
  println(io, "}")
end

ext: jl

max_stars_repo_name: one-more-minute/atom-latex-completions

max_issues_repo_name: one-more-minute/atom-latex-completions


text:
using NeuroCore, Test, Unitful

using NeuroCore: F64Sec, F64Tesla, F64Hz, IntDeg, OneF64Sec, OneF64Tesla, OneF64Hz, OneIntDeg
using NeuroCore: UnkownContrast, IODINE, ipos, jpos

include("bids_tests.jl")

#= TODO
    :stream_offset,
    :auxfiles,
    :srcfile,
    :calmax,
    :calmin,
    :freqdim,
    :phasedim,
    :slicedim,
    :slice_start,
=#


text:
module Examples
using LinearAlgebra
using Random
using DataFrames
using CSV

srcDir = dirname(@__FILE__)*"\\"

include(srcDir *"Utils_states.jl")
include(srcDir *"Examples_sep.jl")
include(srcDir *"Examples_ent.jl")

# using .Utils
using .Utils_states, .Examples_sep, .Examples_ent

export get_examples,
       get_example_overview

"""Returns dict with three keys: ["rand", "sep", "ent"] each entry is also a dict."""
get_examples() = merge(Examples_sep.get_sep_example(),Examples_ent.get_ent_example())
get_example(ex::String) = get_examples()[ex]
get_example(ex::Int) = get_examples()[[keys(get_examples())...][ex]]

siz_dic = Dict(  4 => (2,2),
                 6 => (2,3),
                 8 => (2,4),
                 9 => (3,3),
                 16 => (4,4))

"""Returns a DataFrame object containing the"""
function get_example_meta(isWrite = false)
    df  = DataFrame( Example = String[], Size = Any[], isSeparable = String[], isReal = Bool[], Bi_rank = Any[], hasPPᵀ = Bool[]);
    for sep in ["sep","ent"]
        ρ_dict = (sep =="sep") ? Examples_sep.get_sep_example() : Examples_ent.get_ent_example()
        for k in keys(ρ_dict)
            ρ      =  ρ_dict[k]
            ex     = k
            siz    = siz_dic[size(ρ)[1]]
            isr    = isreal(ρ)
            rank1  = rank(ρ)
            rank2  = rank(take_Pᵀ(ρ, 2, siz))
            isPPt  = Utils_states.isPPᵀ(ρ,siz)

            push!(df, [ex,siz,sep ,isr,(rank1,rank2),isPPt] )
        end
    end
    sort!(df,[:Example,:isSeparable])
    isWrite ? (CSV.write(dirname(dirname(@__FILE__))*"\\ExSummary.csv", df , delim = "&")) : (return df)
end
#
get_example_meta(ex::Int) = get_example_meta()[ex,:]
get_example_meta(ex::Int) = get_example_meta()[ex,:]

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
# conditional neighbor probability for generally asymmetric SNE
# generalizes kernel and distance
function conditional_neighbor(X::AbstractMatrix, σ²::AbstractVector{<:Real},
                    kernel = gaussian, distance = euclidean)
    x = [c for c in eachcol(X)]
    xt, vt = permutedims(x), permutedims(σ²)
    K = @. kernel(distance(x, xt) / sqrt(vt)) # lazy?
    for i in diagind(K) # "Because we are only interested in modeling pairwise similarities, we set the value of pi|i to zero."
        K[i] = 0
    end
    sumK = sum(K, dims = 1) # this is an O(n²) operation!
    @. K /= sumK # columns are normalized probability distributions now
end

# adjusts variances of conditional neighbor distribution to achieve perplexity u
function optimize_variances(X::AbstractMatrix, u::Real, kernel = gaussian, distance = euclidean;
                            tol::Real = 1e-6, max_iter::Int = 32)
    n = size(X, 2)
    σ² = fill(10one(eltype(X)), n) # if this is too small initially, Newton has problems
    @threads for i in 1:n
        σ²[i] = optimize_variances(X, u, i, σ²[i], kernel, distance, tol = tol, max_iter = max_iter)
    end
    any(isnan, σ²) && throw(DomainError("Warning: variance optimization yielded NaN: $σ²"))
    return σ²
end

# optimizing variances using dense operations
function optimize_variances(X::AbstractMatrix, u::Real, i::Int, σ²::Real,
                            kernel = gaussian, distance = euclidean;
                            tol::Real = 1e-6, max_iter::Int = 32)
    x = [c for c in eachcol(X)]
    n = length(x)
    function objective(log_σ²::Real) # allocating version, is faster for high d, equal for low
        local σ² = exp(log_σ²)
        p = @. kernel(distance(x[i:i], x) / sqrt(σ²)) # this has sparse approximation (see sparse.jl)
        p[i] = 0
        p ./= sum(p)
        entropy(p) - log2(u)
    end
    log_σ² = newtons_method(objective, log(σ²), tol = tol, max_iter = max_iter)
    return exp(log_σ²)
end


text:
module CoordinateSystem

using LinearAlgebra

export CSys

struct CSys
    O::Vector{Float64}
    x::Vector{Float64}
    y::Vector{Float64}
    z::Vector{Float64}
    T::Matrix{Float64} #transorm_matrix
end

function CSys(o::Vector{Float64},p₁::Vector{Float64},p₂::Vector{Float64})
    v₁=p₁-o
    v₂=p₂-o
    if abs(v₁⋅v₂/norm(v₁)/norm(v₂))==1
        error("Two vectors should not be parallel!")
    end
    x=v₁/norm(v₁)
    z=v₁×v₂
    z=z/norm(z)
    y=z×x
    T=[reshape(x,1,3);
        reshape(y,1,3);
        reshape(z,1,3)]
    CSys(o,x,y,z,T)
end

end


text:
using Base.Test

require("XtraMaps.jl");

@everywhere f(x) = begin; sleep(0.5); 2x; end
a, b = ones(20), sones(20);
print("map!: "); @time map!(f,a);
print("pmap!: "); @time pmap!(f,b);
for i=1:20; @test a[i] == b[i]; end

@everywhere g(x, y) = begin; sleep(0.1); 2x-y; end
c, d = ones(50), sones(50);
print("map_with_indices: "); @time x = map_with_indices(g,c);
print("pmap_with_indices: "); @time y = pmap_with_indices(g,c);
print("map_with_indicies!: "); @time map_with_indices!(g,c);
print("pmap_with_indicies!: "); @time pmap_with_indices!(g,d);
for i=1:50
  @test c[i] == d[i] && c[i] == x[i] && c[i] == y[i];
end

@everywhere h(x, y, z) = begin; sleep(0.08); x^2 + 119y - z^2; end
A, B = fill(1.3, (7, 6)), sfill(1.3, (7, 6));
print("mmap_with_indices: "); @time C = mmap_with_indices(h,A);
print("pmmap_with_indices: "); @time D = pmmap_with_indices(h,A);
print("mmap_with_indices!: "); @time mmap_with_indices!(h,A);
print("pmmap_with_indices!: "); @time pmmap_with_indices!(h,B);
for j=1:6, i=1:7
  @test_approx_eq A[i,j] B[i,j]
  @test_approx_eq A[i,j] C[i,j]
  @test_approx_eq A[i,j] D[i,j]
end


text:
using Tokenize

@testset "tree_functions" begin

    @testset "parse_tree format 1" begin
        tree = parse_tree("data/S10593.nex")
        @test length(tree) == 4

        nodes, taxa_labels, character_labels, title = tree
        @test typeof(nodes) == Array{Dict{String,Any},1}
        @test typeof(taxa_labels) == Dict{String,String}
        @test typeof(character_labels) == Dict{String,String}
        @test typeof(title) == String

        @test length(nodes[1]["Taxa"]) == 62
        @test nodes[1]["Descendents"][1] == 2
        @test nodes[1]["Groups"] == 2

        @test taxa_labels["32"] == "Tubificoides_parapectinatus_CE2490"
        @test taxa_labels["10"] == "Tubificoides_amplivasatus_CE1731"

        @test character_labels["Tubificoides_amplivasatus_CE1731"] == """---------------------CTGC-----GGAAGATCATT------------CCCGTTTCCACACACGATCGACAC-ACGTGCTTTAACGGAG--------GTACAGATCGCGTCGTG--------------------ATTACATTTGAACGCAATTCCGTTCGGCTGTCC-GACAACAGCAGGGCAGTCGT---------------------------------------------------GAGTGGTATGTGGGCCTGGGCAACCCACAAGGGCAGCTAG-------------------------------------------------------------------------------------------------CCTAGCAACGGCTACGATGCTCCCTCCCAACAGTGTCGAGTGATGGTGGAGAA------CGAGGAGC----------GTAACATGCACGAAAGAGTGCA-AGAACGAG-----------------------------TGCACAGACCAAGGAAAGTCAGAGGCAAGCTGATACGAAAGGAAACTGCGAAAGGAAGTACGGAGACTGTGCGAATGA-------------------------------------------ATTATCCCATGTTACCGAATA-----------CGATGAA-------TAC---TCGTTGAAA-CTTCGGGGTACCTGTCGT--------------------------------------GTCGACT-------CTTGC--------------------------------------------ACGATGAGCTGGTCGCGGGCGATC------CCCGCGGCTTCGCGAAGGTTCAAAGAAGCGCGGTT-----------------------CGTCGCTCGCCGT--------CGGTCGGCC-TAAATGCCGGC---------GTGCGAGCAGCGTGCC-------------GCGCTCGCCCCGAAGCCATTTGT-TCGTTTTTTTTCGTACGGTTGCTA--TCGGTCACCTGCCGGTCGACGGA---------------------------------------------------------------------------GGCC-------------------GCCCTTGTTGAACGACTGGGGCTGG--------CCGAAGTCGCGC--------------------------------------------------------------CGCGTA-CGGTGCTGCCCA------GTGTGGAATTTTG--TACAACTCTAAGCGGTGGATCACTCGGCTCGTGTGTCGATGAAGAGCGCAGCCAGCTGCGTGACCTAATGTGAATTGCAGGACACATTGAACATCGATATCTTGAACGCACATTGCGGCCTCGGGTATTCCCGAGGCCACGCCTGTCTCAGGG-TCGGTTGAAACTATCAATCG-TCGGTCGACT---TTTGCCGACGCATTGGAT-GTCGCGGT--GC-GCCGAACTATTA-------------------AGTCGTCGTCGTTCGCGATCGGTCGAAC-TAAATGCCGATCGCGGCGGCGG-----------------CGCA-----GCGCGCCACGTCGTCCGAAGTACAGACGGGT------C-------------------------TGC-GCCTAAGCCTTGCCGACCGCACGTCCT---CAGTCCGTCATT--GCG---------------GAAGGGGACGGC-GACTCGGAAACGT---TGAGAAGCCGTCGCCT---A-TTGCT--TTG-----GCGGTACGGT--------------------CG-GCGTAGGCTGTGTGGCAC--GAC-----------CCAAG-------------------------------AACGAA-------------CGAACGACCGC---------TTGAACGGGGTTGAGCGTCG------ACTCTTGCGTCTTCCTAATTACTCGCGC-------------------------------------------------------------------------------------------------------------------------GTCTTCAT------TCTTCTTCGAC------------------CTGAGATCAGACGAGATTACCCGCTGA---ATTTAA"""

        @test title == "Tubificoides_ITS"

        tree = parse_tree("data/S10593.nex", taxa_to_remove=["Tubificoides_imajimai_CE536", "Tubificoides_benedii_II_CE2692"])
        @test length(tree) == 4

        nodes, taxa_labels, character_labels, title = tree
        @test !(["30", "45"] in nodes[1]["Taxa"])
        @test !haskey(taxa_labels, "30")
        @test !haskey(taxa_labels, "45")
        @test !haskey(character_labels, "Tubificoides_imajimai_CE536")
        @test !haskey(character_labels, "Tubificoides_benedii_II_CE2692")
    end

    @testset "parse_tree format 2" begin
        tree = parse_tree("data/E1E2L1.nex")
        nodes, taxa_labels, character_labels, title = tree

        @test length(nodes) == 421

    end

    @testset "remove_from_tree" begin
        tree = "((1,2,3),((4,5),6));"
        tree_tokens = [untokenize(token) for token in collect(tokenize(tree))]
        remove_from_tree!(tree_tokens,["1"])

        @test !occursin("1", join(tree_tokens))

        tree = "((apple,pinapple,pear),((tomatoe,pepper),potatoe)));"
        tree_tokens = [untokenize(token) for token in collect(tokenize(tree))]
        remove_from_tree!(tree_tokens,["pepper"])
        @test join(tree_tokens) == "((apple,pinapple,pear),((tomatoe),potatoe)));"
        @test !occursin("pepper", join(tree_tokens))

        remove_from_tree!(tree_tokens,["pinapple", "pear"])
        @test join(tree_tokens) == "((apple),((tomatoe),potatoe)));"
        @test !occursin("pinapple", join(tree_tokens))
        @test !occursin("pear", join(tree_tokens))

        tree = "(,(apple,pinapple,pear),((tomatoe,pepper),potatoe)));"
        tree_tokens = [untokenize(token) for token in collect(tokenize(tree))]
        remove_from_tree!(tree_tokens,["pinapple", "pear"])
        @test !occursin("pear", join(tree_tokens))

        tree = "((1,2,3),(a),((4,5),6));"
        tree_tokens = [untokenize(token) for token in collect(tokenize(tree))]
        remove_from_tree!(tree_tokens,["a"])

        @test join(tree_tokens) == "((1,2,3),((4,5),6));"
        @test !occursin("a", join(tree_tokens))

    end

    @testset "get_nodes" begin
        tree = "(1,((2,3),((4,5),6)));"

        nodes = get_nodes(tree)

        @test length(nodes) == 5
        @test nodes[1]["Taxa"] == ["1", "2", "3", "4", "5", "6"]
        @test nodes[2]["Taxa"] == ["2", "3", "4", "5", "6"]
        @test nodes[3]["Taxa"] == ["2", "3"]
        @test nodes[4]["Taxa"] == ["4", "5", "6"]
        @test nodes[5]["Taxa"] == ["4", "5"]

        @test nodes[1]["Descendents"] == [2]
        @test nodes[2]["Descendents"] == [3,4]
        @test nodes[3]["Descendents"] == []
        @test nodes[4]["Descendents"] == [5]
        @test nodes[5]["Descendents"] == []

        @test nodes[1]["Groups"] == 2
        @test nodes[2]["Groups"] == 2
        @test nodes[3]["Groups"] == 2
        @test nodes[4]["Groups"] == 2
        @test nodes[5]["Groups"] == 2

        tree = "((1,2,3),((4,5),6)));"
        nodes = get_nodes(tree)

        @test length(nodes) == 4
        @test nodes[1]["Taxa"] == ["1", "2", "3", "4", "5", "6"]
        @test nodes[2]["Taxa"] == ["1", "2", "3"]
        @test nodes[3]["Taxa"] == ["4", "5", "6"]
        @test nodes[4]["Taxa"] == ["4", "5"]

        @test nodes[1]["Descendents"] == [2,3]
        @test nodes[2]["Descendents"] == []
        @test nodes[3]["Descendents"] == [4]
        @test nodes[4]["Descendents"] == []

        @test nodes[1]["Groups"] == 2
        @test nodes[2]["Groups"] == 3
        @test nodes[3]["Groups"] == 2
        @test nodes[4]["Groups"] == 2

        tree = "((1,2,3),((4,5),6)));"
        nodes = get_nodes(tree, taxa_to_remove=["2"])

        @test length(nodes) == 4
        @test nodes[1]["Taxa"] == ["1", "3", "4", "5", "6"]
        @test nodes[2]["Taxa"] == ["1", "3"]
        @test nodes[3]["Taxa"] == ["4", "5", "6"]
        @test nodes[4]["Taxa"] == ["4", "5"]

        @test nodes[1]["Descendents"] == [2,3]
        @test nodes[2]["Descendents"] == []
        @test nodes[3]["Descendents"] == [4]
        @test nodes[4]["Descendents"] == []

        @test nodes[1]["Groups"] == 2
        @test nodes[2]["Groups"] == 2
        @test nodes[3]["Groups"] == 2
        @test nodes[4]["Groups"] == 2

        tree = "((1,2,3),((4,5),6)));"
        nodes = get_nodes(tree, taxa_to_remove=["1", "2", "3"])

        @test length(nodes) == 3
        @test nodes[1]["Taxa"] == []
        @test nodes[2]["Taxa"] == ["4", "5", "6"]
        @test nodes[3]["Taxa"] == ["4", "5"]

        @test nodes[1]["Descendents"] == []
        @test nodes[2]["Descendents"] == [3]
        @test nodes[3]["Descendents"] == []

        @test nodes[1]["Groups"] == 1
        @test nodes[2]["Groups"] == 2
        @test nodes[3]["Groups"] == 2

    end

    @testset "add_nodes!" begin

    end
end


text:
function uploadCommandList(rhd::FPGA,commandList, auxCommandSlot, bank)

    #error checking goes here

    for i=1:length(commandList)

        SetWireInValue(rhd,WireInCmdRamData, commandList[i])
        SetWireInValue(rhd,WireInCmdRamAddr, i-1)
        SetWireInValue(rhd,WireInCmdRamBank, bank)
        UpdateWireIns(rhd)
        if auxCommandSlot == "AuxCmd1"
		    if (OPEN_EPHYS)
                ActivateTriggerIn(rhd,TrigInRamWrite,0)
		    else
	            ActivateTriggerIn(rhd,TrigInConfig,1)
	        end
        elseif auxCommandSlot == "AuxCmd2"
            if (OPEN_EPHYS)
                ActivateTriggerIn(rhd,TrigInRamWrite,1)
            else
                ActivateTriggerIn(rhd,TrigInConfig,2)
            end
        elseif auxCommandSlot == "AuxCmd3"
            if (OPEN_EPHYS)
                ActivateTriggerIn(rhd,TrigInRamWrite,2)
            else
                ActivateTriggerIn(rhd,TrigInConfig,3)
            end
        end
    end

    nothing
end

function selectAuxCommandBank(rhd::FPGA,port, commandslot, bank)

    #Error checking goes here

    if port=="PortA"
        bitShift=0
    elseif port=="PortB"
        bitShift=4
    elseif port=="PortC"
        bitShift=8
    elseif port=="PortD"
        bitShift=12
    end

    if commandslot=="AuxCmd1"
        SetWireInValue(rhd,WireInAuxCmdBank1,(bank<<bitShift),(0x000f<<bitShift))
    elseif commandslot=="AuxCmd2"
        SetWireInValue(rhd,WireInAuxCmdBank2,(bank<<bitShift),(0x000f<<bitShift))
    elseif commandslot=="AuxCmd3"
        SetWireInValue(rhd,WireInAuxCmdBank3,(bank<<bitShift),(0x000f<<bitShift))
    end

    UpdateWireIns(rhd)

    nothing
end

function selectAuxCommandLength(rhd::FPGA,commandslot,loopIndex,endIndex)

    #Error checking goes here

    if commandslot=="AuxCmd1"
        if (OPEN_EPHYS)
            SetWireInValue(rhd,WireInAuxCmdLoop1,loopIndex)
            SetWireInValue(rhd,WireInAuxCmdLength1,endIndex)
        else
            SetWireInValue(rhd,WireInAuxCmdLoop, loopIndex, 0x000003ff);
            SetWireInValue(rhd,WireInAuxCmdLength, endIndex, 0x000003ff);
        end
    elseif commandslot=="AuxCmd2"
        if (OPEN_EPHYS)
            SetWireInValue(rhd,WireInAuxCmdLoop2,loopIndex)
            SetWireInValue(rhd,WireInAuxCmdLength2,endIndex)
        else
            SetWireInValue(rhd, WireInAuxCmdLoop, loopIndex << 10, 0x000003ff << 10);
            SetWireInValue(rhd, WireInAuxCmdLength, endIndex << 10, 0x000003ff << 10);
        end
    elseif commandslot=="AuxCmd3"
        if (OPEN_EPHYS)
            SetWireInValue(rhd,WireInAuxCmdLoop3,loopIndex)
            SetWireInValue(rhd,WireInAuxCmdLength3,endIndex)
        else
            SetWireInValue(rhd,WireInAuxCmdLoop, loopIndex << 20, 0x000003ff << 20);
            SetWireInValue(rhd,WireInAuxCmdLength, endIndex << 20, 0x000003ff << 20);
        end
    end

    UpdateWireIns(rhd)

    nothing
end


text:
#####################################################################################
# Hamiltonians
#####################################################################################

"Build the unbounded superconducting hamiltonian (leads)"
function hS(p = Params())
    modelS = models_device(p).modelS
    lat = lat_infinite(p)
    return hamiltonian(lat, modelS, orbitals= (:eup,:edown,:hup,:hdown))
end

function nanoribbon_irregular_edge(p = Params())
    @unpack edges, normal = p
    m = models_device(p)
    hterms =  m.modelNtop + m.modelNbot + m.model0
    lat = unitcell(lat_infinite(p), (1, 3); region = r -> 0 < r[2] < 20)
    return lat |> hamiltonian(hterms, orbitals=  normal ? (:eup, :edown) : (:eup,:edown,:hup,:hdown)) |>
        parametric(m.onsiteNtop_ribbon!, m.onsiteNbot_ribbon!)
end

" Constructor returns a SNS junction y infinite (infinite SC-N armchair interface)
and a normal unbounded parametric hamiltonian as 2D ham in 3D space.
IMPORTANT SMOOTH MUST BE SET TO FALSE or dir = :y "
function nanoribbonperiodic(p = Params())
    @unpack normal = p
    lat = lat_infinite(p)
    m = models_device(p)
    devicereg, = xyregs(p)
    br1, br2, br3 = nanoribbonbravais(p)
    hterms =  m.modelNtop + m.modelNbot + m.model00
    hamN = lat |>
        hamiltonian(hterms, orbitals = normal ? (:eup, :edown) : (:eup,:edown,:hup,:hdown))
    phamN = hamN |> parametric(m.rashba!, m.onsiteNtop_inf!, m.onsiteNbot_inf!,
        m.Δdimer_top!, m.Δdimer_bot!,m.cantingNAtBb!, m.cantingNAbBt!)
    hamSNS = hamN |> unitcell(br1, br2) |> unitcell(br3, region = devicereg)
    phamSNS = hamSNS |>
        parametric(m.onsiteNtop_ribbon!, m.onsiteNbot_ribbon!, m.field!,
            m.Δdimer_top!,m.Δdimer_bot!, m.cantingNAtBb!, m.cantingNAbBt!, m.onsiteS!,
            m.rashba!, m.interhopSCt1!, m.interhopSCt3!, m.peierlsnoBperp!,
            m.transparency!, m.localselfenergy!, m.analyticalselfenergyon!,
            m.scphasesnoBperp!, m.analyticalselfenergyhop!)
    return phamSNS, phamN
end

"Build a nanoribon of width W with params P::Params and infinite edge (armchair or zigzag)
 to the vacuum"
function nanoribbon(p = Params())
    @unpack edges, normal = p
    m = models_device(p)
    hterms =  m.modelNtop + m.modelNbot + m.model0
    nanoribreg = xyregs(p)[2]
    lat = unitcell(lat_infinite(p), ifelse(edges == :armchair_zigzag,(-1,1),(1,1)), region = nanoribreg)
    return lat |> hamiltonian(hterms, orbitals=  normal ? (:eup, :edown) : (:eup,:edown,:hup,:hdown)) |>
        parametric(m.onsiteNtop_ribbon!, m.onsiteNbot_ribbon!)
end

"Build the bounded parametric hamiltonian given `p::Params`, free of dangling bonds.
if `mask = true` we create a rectangular hole in the center of the device whose dimensions
are specified by `Δx_mask, Δy_mask`"
function paramhams(p = Params(); flat = true)
    @unpack t0, t1, t3, λ, normal, mask , Δx_mask, Δy_mask, Ln, W, a0, scale = p
    m = models_device(p)
    lat = lat_infinite(p)
    resetonsites! = @onsite!(o -> 0*o)
    resethoppings! = @hopping!(t -> 0*t)
    reg = Quantica.combine(xyregs(p)[1], xyregs(p)[2])
    mask == false ?  devicereg = reg : devicereg = xor(reg, RegionPresets.rectangle((Ln-2*Δx_mask, W-2*Δy_mask), (0., W/2 + a0 *scale)))
    ham = lat |> hamiltonian(m.model0, orbitals =  normal ? (:eup, :edown) : (:eup,:edown,:hup,:hdown)) |>
            unitcell(mincoordination = coordinationfunc(t0,t1,t3,λ), region = devicereg)
            #note when λ = 0 coordinationfunc stops working correctly, put an infinitesimal but finite value
    pham = ham |>
        parametric(m.onsiteNtop!, m.onsiteNbot!, m.Δdimer_top!, m.Δdimer_bot!, m.cantingNAtBb!,
        m.cantingNAbBt!, m.onsiteS!, m.interhopSCt1!, m.interhopSCt3!, m.peierls!,
        m.transparency!, m.localselfenergy!, m.analyticalselfenergyon!, m.scphases!,
        m.analyticalselfenergyhop!)
    pdham = ham |> parametric(resetonsites!, resethoppings!, m.pairing!, m.dpeierls!) # old: m.onsiteS!, before dpeierls
    return pham, pdham
end

"Implementation of a monolayer graphene model with Kane-Mele SOC + SC -> Kitaev instead of
Bilayer + Rashba + Zeeman + SC -> Oreg"
function Fukane end

"Determines the minimum connectivity in the lattice so that no dangling bonds form. It reads
how many nonzero hopping elements there are (up to a given precision) to return the minimum
coordination number in such a way we avoid dangling sites in the corners"
coordinationfunc(t0::Number, t1::Number, t3::Number, λ::Number) = coordinationfunc([t0, t1, t3, λ])
coordinationfunc(hopvect::Array) = length(hopvect[abs.(hopvect) .> 1e-10]) + 1


text:
module GeometricStructures

include("points.jl")
include("linesegments.jl")
include("convexhull.jl")

export
    AbstractPoint,
    AbstractPoint2D,
    Point2D,
    LineSegment2D,
    convexhull

end # module

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
"""
    AbstractOptimizer

Base type for all optimizers.
"""
abstract AbstractOptimizer

"""
    AbstractLearningRateScheduler

Base type for all learning rate scheduler.
"""
abstract AbstractLearningRateScheduler

"""
    AbstractMomentumScheduler

Base type for all momentum scheduler.
"""
abstract AbstractMomentumScheduler

"""
    OptimizationState

# Attributes:
* `batch_size`: The size of the mini-batch used in stochastic training.
* `curr_epoch`:
  The current epoch count. Epoch 0 means no training yet, during the first
  pass through the data, the epoch will be 1; during the second pass, the
  epoch count will be 1, and so on.
* `curr_batch`:
  The current mini-batch count. The batch count is reset during every epoch.
  The batch count 0 means the beginning of each epoch, with no mini-batch
  seen yet. During the first mini-batch, the mini-batch count will be 1.
* `curr_iter`:
  The current iteration count. One iteration corresponds to one mini-batch,
  but unlike the mini-batch count, the iteration count does **not** reset
  in each epoch. So it track the *total* number of mini-batches seen so far.
"""
type OptimizationState
  batch_size :: Int
  curr_epoch :: Int
  curr_batch :: Int
  curr_iter  :: Int
end
OptimizationState(batch_size::Int) = OptimizationState(batch_size, 0, 0, 0)

"""
    get_learning_rate(scheduler, state)

# Arguments
* `scheduler::AbstractLearningRateScheduler`: a learning rate scheduler.
* `state::OptimizationState`: the current state about epoch, mini-batch and iteration count.

Returns the current learning rate.
"""
function get_learning_rate end

################################################################################
# The learning rate module
module LearningRate
import ..mx: AbstractLearningRateScheduler, OptimizationState, get_learning_rate

"""
    LearningRate.Fixed

Fixed learning rate scheduler always return the same learning rate.
"""
type Fixed <: AbstractLearningRateScheduler
  learning_rate :: Float64
end
get_learning_rate(self :: Fixed, state :: OptimizationState) = self.learning_rate

"""
    LearningRate.Exp

``\eta_t = \eta_0\gamma^t``. Here ``t`` is the epoch count, or the iteration
count if `decay_on_iteration` is set to true.
"""
type Exp <: AbstractLearningRateScheduler
  learning_rate :: Float64
  gamma         :: Float64
  on_iteration  :: Bool
end
function Exp(base_lr::Real; gamma::Real=0.9, decay_on_iteration::Bool=false)
  @assert(0 < gamma < 1)
  Exp(Float64(base_lr), Float64(gamma), decay_on_iteration)
end
get_learning_rate(self :: Exp, state :: OptimizationState) =
    self.learning_rate * self.gamma ^ (self.on_iteration ? state.curr_iter : state.curr_epoch)
"""
    LearningRate.Inv

``\eta_t = \eta_0 * (1 + \gamma * t)^(-power)``.
Here ``t`` is the epoch count, or the iteration count if `decay_on_iteration`
is set to true.
"""
type Inv <: AbstractLearningRateScheduler
  learning_rate :: Float64
  gamma         :: Float64
  power         :: Float64
  on_iteration  :: Bool
end
function Inv(base_lr :: Real; gamma::Real=0.9, power::Real=0.5, decay_on_iteration::Bool=false)
  @assert(0 < gamma < 1)
  @assert(0 <= power)
  Inv(Float64(base_lr), Float64(gamma), Float64(power), decay_on_iteration)
end
get_learning_rate(self :: Inv, state :: OptimizationState) =
  self.learning_rate * ( 1 + self.gamma * (self.on_iteration ? state.curr_iter : state.curr_epoch)) ^ (-self.power)
end# module LearningRate
################################################################################
function get_lr_scheduler(scheduler :: Any, lr :: Real)
  if isa(scheduler, AbstractLearningRateScheduler)
    return scheduler
  else
    return LearningRate.Fixed(lr)
  end
end

"""
    get_momentum(scheduler, state)

* `scheduler::AbstractMomentumScheduler`: the momentum scheduler.
* `state::OptimizationState`: the state about current epoch, mini-batch and iteration count.

Returns the current momentum.
"""
function get_momentum
end

################################################################################
# The Momentum module
module Momentum
import ..mx: AbstractMomentumScheduler, OptimizationState, get_momentum

"""
    Momentum.Null

The null momentum scheduler always returns 0 for momentum. It is also used to
explicitly indicate momentum should not be used.
"""
type Null <: AbstractMomentumScheduler
end
get_momentum(self :: Null, state :: OptimizationState) = 0.0

"""
    Momentum.Fixed

Fixed momentum scheduler always returns the same value.
"""
type Fixed <: AbstractMomentumScheduler
  momentum :: Float64
end
get_momentum(self :: Fixed, state :: OptimizationState) = self.momentum
end # module Momentum
################################################################################
function get_momentum_scheduler(scheduler :: Any, momentum :: Real)
  if isa(scheduler, AbstractMomentumScheduler)
    return scheduler
  elseif momentum == 0
    return Momentum.Null()
  else
    return Momentum.Fixed(momentum)
  end
end

"""
    get_updater(optimizer)

A utility function to create an updater function, that uses its closure to
store all the states needed for each weights.

* `optimizer::AbstractOptimizer`: the underlying optimizer.
"""
function get_updater(optimizer :: AbstractOptimizer)
  states = Dict{Int,Any}()
  function updater(index :: Int, grad :: NDArray, weight :: NDArray)
    if !haskey(states, index)
      states[index] = create_state(optimizer, index, weight)
    end
    update(optimizer, index, weight, grad, states[index])
  end
  return updater
end

"""
    AbstractOptimizerOptions

Base class for all optimizer options.
"""
abstract AbstractOptimizerOptions

"""
    normalized_gradient(opts, state, grad)

* `opts::AbstractOptimizerOptions`: options for the optimizer, should contain the field
          `grad_scale`, `grad_clip` and `weight_decay`.
* `state::OptimizationState`: the current optimization state.
* `weight::NDArray`: the trainable weights.
* `grad::NDArray`: the original gradient of the weights.

   Get the properly normalized gradient (re-scaled and clipped if necessary).
"""
function normalized_gradient(opts::AbstractOptimizerOptions, state::OptimizationState,
                             weight::NDArray, grad::NDArray)
  grad_scale = 1.0 / state.batch_size

  grad = grad_scale * grad
  if opts.grad_clip > 0
    grad = clip(grad, -opts.grad_clip, opts.grad_clip)
  end
  @inplace grad += opts.weight_decay * weight

  return grad
end

include("optimizers/sgd.jl")
include("optimizers/adam.jl")


text:
function unit(::Type{T}, i::Integer) where {T <: StaticVector}
    tup = ntuple(Val(length(T))) do j
        return ifelse(i == j, 1, 0)
    end
    return T(tup)
end

macro fixed_vector(name, parent)
    expr = quote
        struct $(name){S,T} <: $(parent){S,T}
            data::NTuple{S,T}

            function $(name){S,T}(x::NTuple{S,T}) where {S,T}
                return new{S,T}(x)
            end

            function $(name){S,T}(x::NTuple{S,Any}) where {S,T}
                return new{S,T}(StaticArrays.convert_ntuple(T, x))
            end
        end

        size_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name){S,T} where S}, or) where {T} = T
        eltype_or(::Type{$(name){S,T} where T}, or) where {S} = or
        eltype_or(::Type{$(name){S,T}}, or) where {S,T} = T

        size_or(::Type{$(name){S,T} where S}, or) where {T} = or
        size_or(::Type{$(name){S,T} where T}, or) where {S} = Size{(S,)}()
        size_or(::Type{$(name){S,T}}, or) where {S,T} = (S,)

        # Array constructor
        function $(name){S}(x::AbstractVector{T}) where {S,T}
            @assert S <= length(x)
            return $(name){S,T}(ntuple(i -> x[i], Val(S)))
        end

        function $(name){S,T1}(x::AbstractVector{T2}) where {S,T1,T2}
            @assert S <= length(x)
            return $(name){S,T1}(ntuple(i -> convert(T1, x[i]), Val(S)))
        end

        function $(name){S,T}(x) where {S,T}
            return $(name){S,T}(ntuple(i -> convert(T, x), Val(S)))
        end

        $(name){S}(x::T) where {S,T} = $(name){S,T}(ntuple(i -> x, Val(S)))
        $(name){1,T}(x::T) where {T} = $(name){1,T}((x,))
        $(name)(x::NTuple{S}) where {S} = $(name){S}(x)
        function $(name)(x::T) where {S,T <: Tuple{Vararg{Any,S}}}
            return $(name){S,StaticArrays.promote_tuple_eltype(T)}(x)
        end

        function $(name){S}(x::T) where {S,T <: Tuple}
            return $(name){S,StaticArrays.promote_tuple_eltype(T)}(x)
        end
        $(name){S,T}(x::StaticVector) where {S,T} = $(name){S,T}(Tuple(x))

        @generated function (::Type{$(name){S,T}})(x::$(name)) where {S,T}
            idx = [:(x[$i]) for i in 1:S]
            return quote
                $($(name)){S,T}($(idx...))
            end
        end

        @generated function Base.convert(::Type{$(name){S,T}}, x::$(name)) where {S,T}
            idx = [:(x[$i]) for i in 1:S]
            return quote
                $($(name)){S,T}($(idx...))
            end
        end

        @generated function (::Type{SV})(x::StaticVector) where {SV <: $(name)}
            len = size_or(SV, size(x))[1]
            return if length(x) == len
                :(SV(Tuple(x)))
            elseif length(x) > len
                elems = [:(x[$i]) for i in 1:len]
                :(SV($(Expr(:tuple, elems...))))
            else
                error("Static Vector too short: $x, target type: $SV")
            end
        end

        Base.@pure StaticArrays.Size(::Type{$(name){S,Any}}) where {S} = Size(S)
        Base.@pure StaticArrays.Size(::Type{$(name){S,T}}) where {S,T} = Size(S)

        Base.@propagate_inbounds function Base.getindex(v::$(name){S,T}, i::Int) where {S,T}
            return v.data[i]
        end

        Base.Tuple(v::$(name)) = v.data
        function Base.convert(::Type{$(name){S,T}}, x::NTuple{S,T}) where {S,T}
            return $(name){S,T}(x)
        end
        function Base.convert(::Type{$(name){S,T}}, x::Tuple) where {S,T}
            return $(name){S,T}(convert(NTuple{S,T}, x))
        end

        @generated function StaticArrays.similar_type(::Type{SV}, ::Type{T},
                                                      s::Size{S}) where {SV <: $(name),T,S}
            return if length(S) === 1
                $(name){S[1],T}
            else
                StaticArrays.default_similar_type(T, s(), Val{length(S)})
            end
        end

        Base.:(*)(a::$name, b::$name) = a .* b
        function Base.broadcasted(f, a::AbstractArray{T}, b::$name) where {T <: $name}
            return Base.broadcasted(f, a, (b,))
        end
    end
    return esc(expr)
end

abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
@fixed_vector Point AbstractPoint
@fixed_vector Vec StaticVector

const Mat = SMatrix
const VecTypes{N,T} = Union{StaticVector{N,T},NTuple{N,T}}
const Vecf0{N} = Vec{N,Float32}
const Pointf0{N} = Point{N,Float32}
Base.isnan(p::Union{AbstractPoint,Vec}) = any(x -> isnan(x), p)

# Create constes like Mat4f0, Point2, Point2f0
for i in 1:4
    for T in [:Point, :Vec]
        name = Symbol("$T$i")
        namef0 = Symbol("$T$(i)f0")
        @eval begin
            const $name = $T{$i}
            const $namef0 = $T{$i,Float32}
            export $name
            export $namef0
        end
    end
    name = Symbol("Mat$i")
    namef0 = Symbol("Mat$(i)f0")
    @eval begin
        const $name{T} = $Mat{$i,$i,T,$(i * i)}
        const $namef0 = $name{Float32}
        export $name
        export $namef0
    end
end

export Mat, Vec, Point, unit
export Vecf0, Pointf0


text:
import MitosisStochasticDiffEq as MSDE
using Mitosis, DiffEqNoiseProcess, StochasticDiffEq
using Test, Random
using LinearAlgebra

@testset "Internal solver tests" begin
  @testset "EulerMaruyama == EM test" begin
    @testset "n=1, m=1 (univariate, scalar) tests" begin
      # define SDE function
      f(u,p,t) = p[1]*u + p[2]
      g(u,p,t) = p[3]*u

      # time span
      tstart = 0.0
      tend = 1.0
      dt = 0.02
      trange = tstart:dt:tend

      B, β, σ̃ = -0.1, 0.2, 1.3
      p = [B, β, σ̃]

      u0 = rand(1)

      kernel = MSDE.SDEKernel(Mitosis.AffineMap(B, β), Mitosis.ConstantMap(σ̃), trange, p)
      # pass noise process and compare with EM()
      Ws = cumsum([[zero(u0)];[sqrt(trange[i+1]-ti)*randn(size(u0))
              for (i,ti) in enumerate(trange[1:end-1])]])
      NG = NoiseGrid(trange,Ws)

      uendEM, (tsEM, uEM, noiseEM) = MSDE.sample(kernel, u0, EM(false), NG)
      uend1, (ts1, u1, noise1) = MSDE.sample(kernel, u0, MSDE.EulerMaruyama!(), Ws)

      @test tsEM ≈ ts1 rtol=1e-12
      @test uEM ≈ u1 rtol=1e-12
      @test uendEM ≈ uend1 rtol=1e-12
    end
  end

  @testset "strong convergence tests" begin
    @testset "n=1, m=1 (univariate, scalar) tests" begin
      #..
    end
  end

  @testset "weak convergence tests" begin
    @testset "n=1, m=1 (univariate, scalar) tests" begin
      # ..
    end
  end

end


text:
import AQFED.TermStructure: ConstantBlackModel
import AQFED.Black: blackScholesFormula

export AndersenLakeNRepresentationJH
function AndersenLakeNRepresentationJH(
    model::ConstantBlackModel,
    tauMax::Float64,
    atol::Float64,
    nC::Int,
    nIter::Int,
    nTS1::Int,
    nTS2::Int;
    isCall::Bool = false,
)
    if iseven(nTS1)
        throw(DomainError(string("nTS1 must be odd but was ", nTS1)))
    end
    if iseven(nTS2)
        throw(DomainError(string("nTS2 must be odd but was ", nTS2)))
    end
    K = 1.0
    avec = zeros(nC + 1)
    qvec = zeros(nC + 1)
    avecD = zeros(nC + 1)
    qvecD = zeros(nC + 1)
    wvec = zeros(nTS1)
    yvec = zeros(nTS1)
    ndiv2 = trunc(Int, (nTS1 + 1) / 2)
    hn = tanhsinhStep(nTS1)
    hvec = hn .* (ndiv2-1:-1:1)
    svec = @. pi * sinh(hvec) / 2
    @. @view(yvec[1:ndiv2-1]) = tanh(svec) #may be more precise to store y+1 directly instead
    @. @view(wvec[1:ndiv2-1]) = hn * pi * cosh(hvec) / (2 * cosh(svec)^2)
    for i = 1:ndiv2-1
        yvec[nTS1+1-i] = -yvec[i]
        wvec[nTS1+1-i] = wvec[i]
    end
    yvec[ndiv2] = 0
    wvec[ndiv2] = pi * hn / 2
    r = model.r
    q = model.q
    modelB = model
    if isCall  #use McDonald and Schroder symmetry
        r, q = q, r
        modelB = ConstantBlackModel(vol, r, q)
    end
    vol = model.vol

    capX = K
    capXD = K * r / q
    logCapX = log(capX)
    logCapXD = log(capXD)
    local fprev = capX
    local fprevD = capXD
    qvec[nC+1] = 0
    qvecD[nC+1] = 0
    for i = nC:-1:1
        zi = cos((i - 1) * pi / nC)
        taui = tauMax / 4 * (1 + zi)^2
        fi = americanBoundaryPutQDP(false, modelB, fprev, K, taui, atol)
        fprev = fi
        qvec[i] = (log(fi / capX))^2
        fi = americanBoundaryPutQDP(true, modelB, fprevD, K, taui, atol)
        fprevD = fi
        qvecD[i] = (log(fi / capXD))^2
    end
    B = computeBoundaryFromQVec(-1, qvec, capX)
    Bdown = computeBoundaryFromQVec(1, qvecD, capXD)
    crossingIndex = updateCrossingBoundariesALUp!(nC, B, capX, Bdown, capXD)
    tauMaxOrig = tauMax
    if crossingIndex >= 1
        #refine crossingIndex
        i = crossingIndex
        zi = cos(i * pi / nC)
        taui = tauMax / 4 * (1 + zi)^2
        tauStart = taui
        fStart = B[i]
        fStartD = Bdown[i]
        i += 1
        zi = cos(i * pi / nC)
        taui = tauMax / 4 * (1 + zi)^2
        tauEnd = taui
        fEnd = B[i]
        fEndD = Bdown[i]

        tauCrossing = searchCrossingTimeQDPN(
            tauEnd,
            fEnd,
            fEndD,
            tauStart,
            fStart,
            fStartD,
            K,
            model,
            atol,
        )

        tauMax = tauCrossing
        fprev = capX
        fprevD = capXD
        for i = nC:-1:1
            zi = cos((i - 1) * pi / nC)
            taui = tauMax / 4 * (1 + zi)^2
            fi = americanBoundaryPutQDP(false, modelB, fprev, K, taui, atol)
            fprev = fi
            qvec[i] = (log(fi / capX))^2
            fi = americanBoundaryPutQDP(true, modelB, fprevD, K, taui, atol)
            fprevD = fi
            qvecD[i] = (log(fi / capXD))^2
        end
        B = computeBoundaryFromQVec(-1, qvec, capX)
        Bdown = computeBoundaryFromQVec(1, qvecD, capXD)
    end ##of crossing refinement
    @. qvec = log(B / capX)^2
    @. qvecD = log(Bdown / capXD)^2
    tauVector = zeros(nTS1)
    k1 = zeros(nTS1)
    k2 = zeros(nTS1)
    d1Vector = zeros(nTS1)
    d2Vector = zeros(nTS1)
    #println("AL init ", "Up ", computeBoundaryFromQVec(-1, qvec, capX), " Down Boundary ", computeBoundaryFromQVec(1, qvecD, capXD))
    for j = 1:nIter
        updateAvec!(avec, nC, qvec)
        ##up boundary
        updateAvec!(avecD, nC, qvecD)
        for i = 1:nC
            zi = cos((i - 1) * pi / nC)
            taui = tauMax / 4 * (1 + zi)^2
            Kstari = K * exp(-(r - q) * taui)
            lnBtaui = logCapX - sqrt(qvec[i])
            sum1k = 0.0
            sum2k = 0.0
            @. tauVector = taui / 4 * (1 + yvec)^2
            crossIndex = 0
            @inbounds for sk1 = 1:nTS1
                if yvec[sk1] != -1
                    tauk = tauVector[sk1]
                    zck = 2 * sqrt((taui - tauk) / tauMax) - 1
                    qckD = chebQck(avecD, zck)
                    qck = chebQck(avec, zck)
                    lnBtaukD = logCapXD + sqrt(qckD)
                    lnBtauk = logCapX - sqrt(qck)
                    if (lnBtaukD > lnBtauk) && crossIndex == 0
                        crossIndex = sk1
                    end
                    if crossIndex == 0
                        sqrtv = sqrt(tauk) * vol
                        d1k = ((lnBtaui - lnBtaukD) + (r - q) * tauk) / sqrtv + sqrtv / 2
                        d2k = d1k - sqrtv
                        sum1k -= wvec[sk1] * exp(-q * tauk) * (yvec[sk1] + 1) * normcdf(d1k)
                        sum2k -= wvec[sk1] * exp(-r * tauk) * (yvec[sk1] + 1) * normcdf(d2k)

                        d1k = ((lnBtaui - lnBtauk) + (r - q) * tauk) / sqrtv + sqrtv / 2
                        d2k = d1k - sqrtv
                        sum1k += wvec[sk1] * exp(-q * tauk) * (yvec[sk1] + 1) * normcdf(d1k)
                        sum2k += wvec[sk1] * exp(-r * tauk) * (yvec[sk1] + 1) * normcdf(d2k)
                    end
                end
            end

            sum1k = exp(q * taui) / 2 * taui * sum1k
            sum2k = exp(r * taui) / 2 * taui * sum2k
            sqrtv = sqrt(taui) * vol
            d1i = ((lnBtaui - log(K)) + (r - q) * taui) / sqrtv + sqrtv / 2
            d2i = d1i - sqrtv

            Ni = -1 + exp(r * taui) + normcdf(d2i) + r * sum2k
            Di = -1 + exp(q * taui) + normcdf(d1i) + q * sum1k
            NiOverDi = Ni / Di
            if Di == 0.0 && Ni == 0.0
                #use asymptotic expansion cdf = erfc(-x/sqrt2)/2 and erfc(x) = e^{-x^2}/(x*sqrtpi)*(1-1/(2*x^2))
                NiOverDi = exp(-(d2i^2 - d1i^2) / 2) * (d1i / d2i)
            end
            fi = Kstari * NiOverDi
            if fi < 0
                fi = 1e-4 * Kstari
            end
            lfc = log(fi / capX)
            if isnan(lfc)
                throw(DomainError(
                    fi,
                    string("Nan qvec ", capX, " ", lnBtaui, " ", qvec[i]),
                ))
            end
            qvec[i] = lfc^2
        end
        qvec[nC+1] = 0

        ##Down boundary
        updateAvec!(avec, nC, qvec)
        for i = 1:nC
            zi = cos((i - 1) * pi / nC)
            taui = tauMax / 4 * (1 + zi)^2
            Kstari = K * exp(-(r - q) * taui)
            lnBtauiD = logCapXD + sqrt(qvecD[i])
            sum1kD = 0.0
            sum2kD = 0.0
            @. tauVector = taui / 4 * (1 + yvec)^2
            crossIndex = 0
            @inbounds for sk1 = 1:nTS1
                if yvec[sk1] != -1
                    tauk = tauVector[sk1]
                    zck = 2 * sqrt((taui - tauk) / tauMax) - 1
                    qckD = chebQck(avecD, zck)
                    qck = chebQck(avec, zck)
                    lnBtaukD = logCapXD + sqrt(qckD)
                    lnBtauk = logCapX - sqrt(qck)
                    if (lnBtaukD > lnBtauk) && crossIndex == 0
                        crossIndex = sk1
                    end
                    if crossIndex == 0
                        sqrtv = sqrt(tauk) * vol
                        d1k = ((lnBtauiD - lnBtaukD) + (r - q) * tauk) / sqrtv + sqrtv / 2
                        d2k = d1k - sqrtv
                        sum1kD -=
                            wvec[sk1] * exp(-q * tauk) * (yvec[sk1] + 1) * normcdf(d1k)
                        sum2kD -=
                            wvec[sk1] * exp(-r * tauk) * (yvec[sk1] + 1) * normcdf(d2k)

                        d1k = ((lnBtauiD - lnBtauk) + (r - q) * tauk) / sqrtv + sqrtv / 2
                        d2k = d1k - sqrtv
                        sum1kD +=
                            wvec[sk1] * exp(-q * tauk) * (yvec[sk1] + 1) * normcdf(d1k)
                        sum2kD +=
                            wvec[sk1] * exp(-r * tauk) * (yvec[sk1] + 1) * normcdf(d2k)
                    end
                end
            end

            sum1kD = exp(q * taui) / 2 * taui * sum1kD
            sum2kD = exp(r * taui) / 2 * taui * sum2kD
            sqrtv = sqrt(taui) * vol
            d1i = ((lnBtauiD - log(K)) + (r - q) * taui) / sqrtv + sqrtv / 2
            d2i = d1i - sqrtv
            fiold = capXD * exp(sqrt(qvecD[i]))
            Ni =
                -1 + exp(r * taui) + normcdf(d2i) + r * sum2kD -
                fiold / Kstari * q * (sum1kD)
            Di = -1 + exp(q * taui) + normcdf(d1i)
            NiOverDi = Ni / Di
            if Di == 0.0 && Ni == 0.0
                #use asymptotic expansion cdf = erfc(-x/sqrt2)/2 and erfc(x) = e^{-x^2}/(x*sqrtpi)*(1-1/(2*x^2))
                NiOverDi = exp(-(d2i^2 - d1i^2) / 2) * (d1i / d2i)
            end
            fi = Kstari * NiOverDi
            if fi < 0
                fi = 1e-4 * Kstari
            end
            lfc = log(fi / capXD)
            if isnan(lfc)
                throw(DomainError(
                    fi,
                    string("Nan qvec ", capX, " ", lnBtaui, " ", qvec[i]),
                ))
            end
            qvecD[i] = lfc^2
        end
        qvecD[nC+1] = 0
        ##update boundaries
        B = computeBoundaryFromQVec(-1, qvec, capX)
        Bdown = computeBoundaryFromQVec(1, qvecD, capXD)
        updateCrossingBoundariesALUp!(nC, B, capX, Bdown, capXD)
        @. qvec = (log(B / capX))^2
        @. qvecD = (log(Bdown / capXD))^2
        # updateAvec!(avec, nC, qvec)
        # updateAvec!(avecD, nC, qvecD)
        # println("AL ",j, " Up ", computeBoundaryFromQVec(-1, qvec, capX), " Down Boundary ", computeBoundaryFromQVec(1, qvecD, capXD))
    end

    updateAvec!(avecD, nC, qvecD)
    updateAvec!(avec, nC, qvec)
    if nTS2 != nTS1
        wvec = zeros(nTS2)
        yvec = zeros(nTS2)
        ndiv2 = trunc(Int, (nTS2 + 1) / 2)
        hn = tanhsinhStep(nTS2)
        hi = (ndiv2 - 1) * hn
        for i = 1:ndiv2-1
            p2s = pi * sinh(hi) / 2
            yvec[i] = tanh(p2s)
            p2c = pi * cosh(hi) / 2
            cp2s = cosh(p2s)
            wvec[i] = hn * p2c / (cp2s^2)
            hi = hi - hn
        end
        yvec[ndiv2] = 0
        wvec[ndiv2] = pi * hn / 2
        for i = 1:ndiv2-1
            yvec[nTS2+1-i] = -yvec[i]
            wvec[nTS2+1-i] = wvec[i]
        end
    end
    return AndersenLakeNRepresentation(
        isCall,
        model,
        tauMax,
        tauMaxOrig,
        tauMaxOrig,
        nC,
        nTS1,
        nTS2,
        capX,
        capXD,
        avec,
        avecD,
        qvec,
        qvecD,
        wvec,
        yvec,
    )

    ##up

end

computeBoundaryFromQVec(sign::Int, qvec::Vector{Float64}, capX::Float64) =
    @. exp(sign * sqrt(qvec)) * capX

function updateCrossingBoundariesALUp!(
    m::Int,
    B::Vector{Float64},
    capX::Float64,
    Bdown::Vector{Float64},
    capXD::Float64,
)::Int
    for i = m-1:-1:1
        if B[i+1] < B[i]
            B[i] = B[i+1]
        end
        if Bdown[i] < Bdown[i+1]
            Bdown[i] = Bdown[i+1]
        end
        if B[i] > capX
            B[i] = capX
        end
        if Bdown[i] < capXD
            Bdown[i] = capXD
        end
    end
    icross = m
    while icross > 0 && Bdown[icross] < B[icross]
        icross -= 1
    end
    if icross >= 1
        if Bdown[icross] > capX
            Bdown[icross] = capX
        end
        if B[icross] < capXD
            B[icross] = capXD
        end
        level = B[icross]
        level = max(Bdown[icross+1], level)
        level = min(B[icross+1], level)
        for i = icross:-1:1
            B[i] = level
            Bdown[i] = B[i]
        end
    end
    return icross
end


text:
# SCIP defines for character, string, and numeric constants
const ARTIFICIALVARNAMEPREFIX = "andresultant_"
const TRUE = 1
const FALSE = 0
const SCIP_VERSION = 310
const SCIP_SUBVERSION = 0
const SCIP_COPYRIGHT = "Copyright (c) 2002-2014 Konrad-Zuse-Zentrum fuer Informationstechnik Berlin (ZIB)"
const SCIP_VARTYPE_BINARY_CHAR = 'B'
const SCIP_VARTYPE_INTEGER_CHAR = 'I'
const SCIP_VARTYPE_IMPLINT_CHAR = 'M'
const SCIP_VARTYPE_CONTINUOUS_CHAR = 'C'
const SCIP_LONGINT_FORMAT = "lld"
const SCIP_REAL_FORMAT = "lf"
const SCIP_MAXSTRLEN = 1024
const SCIP_HASHSIZE_PARAMS = 4099
const SCIP_HASHSIZE_NAMES = 131101
const SCIP_HASHSIZE_CUTPOOLS = 131101
const SCIP_HASHSIZE_CLIQUES = 131101
const SCIP_HASHSIZE_NAMES_SMALL = 8011
const SCIP_HASHSIZE_CUTPOOLS_SMALL = 8011
const SCIP_HASHSIZE_CLIQUES_SMALL = 8011
const SCIP_HASHSIZE_VBC = 131101
const SCIP_EXPR_DEGREEINFINITY = 65535

# SCIP defines for known C types and pointers to SCIP structs.
typealias SCIP Void
typealias SCIP_CONS Void
typealias SCIP_VAR Void
typealias SCIP_Real Float64
typealias SCIP_Bool Uint
typealias SCIP_NLROW Void
typealias SCIP_SOL Void
typealias SCIP_EXPRTREE Void
typealias SCIP_PROFILE Void
typealias SCIP_CONSHDLR Void
typealias SCIP_ROW Void
typealias SCIP_SEPA Void
typealias SCIP_HASHMAP Void
typealias SCIP_EXPRGRAPHNODE Void
typealias SCIP_EXPRGRAPH Void
typealias SCIP_QUADVARTERM Void
typealias SCIP_BILINTERM Void
typealias SCIP_NLPI Void
typealias SCIP_NLPIPROBLEM Void
typealias SCIP_BRANCHRULEDATA Void
typealias SCIP_BRANCHRULE Void
typealias SCIP_CONFLICTHDLRDATA Void
typealias SCIP_CONFLICTHDLR Void
typealias SCIP_CONSHDLRDATA Void
typealias SCIP_PROPTIMING Void
typealias SCIP_CONSDATA Void
typealias SCIP_CUT Void
typealias SCIP_CUTPOOL Void
typealias SCIP_DIALOG Void
typealias SCIP_DIALOGHDLR Void
typealias SCIP_DIALOGDATA Void
typealias SCIP_DISPDATA Void
typealias SCIP_DISP Void
typealias SCIP_MESSAGEHDLR Void
typealias SCIP_EVENTHDLR Void
typealias SCIP_EVENTHDLRDATA Void
typealias SCIP_EVENTTYPE Void
typealias SCIP_EVENT Void
typealias SCIP_NODE Void
typealias SCIP_COL Void
typealias SCIP_INTERVAL Void
typealias SCIP_EXPR Void
typealias SCIP_QUADELEM Void
typealias SCIP_EXPRDATA_MONOMIAL Void
typealias SCIP_EXPRINTDATA Void
typealias SCIP_HEURDATA Void
typealias SCIP_HEUR Void
typealias SCIP_HEURTIMING Void
typealias SCIP_HISTORY Void
typealias SCIP_VALUEHISTORY Void
typealias SCIP_CLIQUE Void
typealias SCIP_BASESTAT Void
typealias SCIP_MESSAGEHDLRDATA Void
typealias SCIP_SPARSESOL Void
typealias SCIP_QUEUE Void
typealias SCIP_PQUEUE Void
typealias SCIP_HASHTABLE Void
typealias SCIP_HASHMAPLIST Void
typealias SCIP_RESOURCEACTIVITY Void
typealias SCIP_DIGRAPH Void
typealias SCIP_BT Void
typealias SCIP_BTNODE Void
typealias SCIP_NODESEL Void
typealias SCIP_NODESELDATA Void
typealias SCIP_PARAM Void
typealias SCIP_PARAMDATA Void
typealias SCIP_PRESOLDATA Void
typealias SCIP_PRESOL Void
typealias SCIP_PRICERDATA Void
typealias SCIP_PRICER Void
typealias SCIP_PROPDATA Void
typealias SCIP_PROP Void
typealias SCIP_READERDATA Void
typealias SCIP_READER Void
typealias SCIP_RELAXDATA Void
typealias SCIP_RELAX Void
typealias SCIP_SEPADATA Void
typealias SCIP_DOMCHG Void
typealias SCIP_VARDATA Void
typealias SCIP_HOLELIST Void
typealias SCIP_BDCHGINFO Void
typealias SCIP_BDCHGIDX Void
typealias SCIP_BOUNDCHG Void
typealias SCIP_PROBDATA Void
typealias SCIP_REALARRAY Void
typealias SCIP_LPI Void
typealias SCIP_NLPSTATISTICS Void
typealias SCIP_EVENTDATA Void
typealias SCIP_CLOCK Void
typealias SCIP_INTARRAY Void
typealias SCIP_BOOLARRAY Void
typealias SCIP_PTRARRAY Void


text:
const GITHUB_ORGANIZATION = "UnofficialJuliaMirror"
const GITHUB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
const GITHUB_BOT_PERSONAL_ACCESS_TOKEN = ENV["GITHUB_BOT_PERSONAL_ACCESS_TOKEN"]


text:
# Defines the configuration of the Slave.
# used in the creation of the slave side algorithm.
# defines the path of the slave's local dataset to be loaded

struct SlaveConfiguration
    # for future implementations
    # host::String
    # port::Int64
    filepath::String

    SlaveConfiguration() = new()
    SlaveConfiguration(s::String) = new("")
end


text:
export ε_init, εₛ, εₛ⁻¹, ε_tensor, test_εs, εₘₐₓ, ñₘₐₓ, nₘₐₓ

εᵥ = SHM3( [	1. 	0. 	0.
                0. 	1. 	0.
                0. 	0. 	1.  ]
)

function ε_tensor(n::Float64)
    n² = n^2
    ε = SHM3( [	n²      0. 	    0.
                0. 	    n² 	    0.
                0. 	    0. 	    n²  ]
            )
end

function test_εs(n₁::Float64,n₂::Float64,n₃::Float64)
    ε₁ = ε_tensor(n₁)
    ε₂ = ε_tensor(n₂)
    ε₃ = ε_tensor(n₃)
    return ε₁, ε₂, ε₃
end

function ε_init(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D;Δx=6.,Δy=4.,Nx=64,Ny=64)::Array{SHM3,3}
    g = MaxwellGrid(Δx,Δy,Nx,Ny)
    tree = KDTree(shapes)
    return Float64[ isnothing(findfirst([xx,yy],tree)) ? εᵥ[i,j] : findfirst([xx,yy],tree).data[i,j] for xx=g.x,yy=g.y,i=1:3,j=1:3 ]
end

function ε_init(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D,g::MaxwellData)::Array{SHM3,3}
    tree = KDTree(shapes)
    return Float64[ isnothing(findfirst([xx,yy],tree)) ? εᵥ[i,j] : findfirst([xx,yy],tree).data[i,j] for xx=g.x,yy=g.y,i=1:3,j=1:3 ]
end

function ε_init(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D, x::AbstractVector{Float64}, y::AbstractVector{Float64})::Array{SHM3,3}
    tree = KDTree(shapes)
    return Float64[ isnothing(findfirst([xx,yy],tree)) ? εᵥ[i,j] : findfirst([xx,yy],tree).data[i,j] for xx=x,yy=y,i=1:3,j=1:3 ]
end

function εₛ(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D,tree::KDTree,x::Real,y::Real,δx::Real,δy::Real)::SHM3
    x1,y1 = x+δx/2.,y+δy/2
    x2,y2 = x+δx/2.,y-δy/2
    x3,y3 = x-δx/2.,y-δy/2
    x4,y4 = x-δx/2.,y+δy/2

    s1 = findfirst([x1,y1],tree)
    s2 = findfirst([x2,y2],tree)
    s3 = findfirst([x3,y3],tree)
    s4 = findfirst([x4,y4],tree)

    ε1 = isnothing(s1) ? εᵥ : s1.data
    ε2 = isnothing(s2) ? εᵥ : s2.data
    ε3 = isnothing(s3) ? εᵥ : s3.data
    ε4 = isnothing(s4) ? εᵥ : s4.data

    if (ε1==ε2==ε3==ε4)
        return ε1
    else
        sinds = [ isnothing(ss) ? length(shapes)+1 : findfirst(isequal(ss),shapes) for ss in [s1,s2,s3,s4]]
        s_fg = shapes[min(sinds...)]
        r₀,nout = surfpt_nearby([x, y], s_fg)
        # bndry_pxl[i,j] = 1
        # nouts[i,j,:] = nout
        vxl = (SVector{2,Float64}(x3,y3), SVector{2,Float64}(x1,y1))
        rvol = volfrac(vxl,nout,r₀)
        sind_bg = max(sinds...)
        ε_bg = sind_bg > length(shapes) ? εᵥ : shapes[sind_bg].data
        return SHM3(kottke_avg_param(
                SHM3(s_fg.data),
                SHM3(ε_bg),
                SVector{3,Float64}(nout[1],nout[2],0),
                rvol,))
    end
end

function εₛ(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D,Δx=6.,Δy=4.,Nx=64,Ny=64)::Array{SHM3,3}
    g=MaxwellGrid(Δx,Δy,Nx,Ny)
    tree = KDTree(shapes)
    ε_sm = zeros(Float64,g.Nx,g.Ny,3,3)
    for i=1:g.Nx, j=1:g.Ny
        ε_sm[i,j,:,:] = εₛ(shapes,tree,g.x[i],g.y[j],g.δx,g.δy)
    end
    return ε_sm
end

function εₛ(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D,g::MaxwellGrid)::Array{SHM3,3}
    tree = KDTree(shapes)
    ε_sm = copy(reshape(  [εₛ(shapes,tree,g.x[i],g.y[j],g.δx,g.δy) for i=1:g.Nx,j=1:g.Ny] , (g.Nx,g.Ny,1)) )
end

# function εₛ⁻¹(shapes::AbstractVector{T},g::MaxwellGrid) where T <: GeometryPrimitives.Shape{2,4,D} where D
function εₛ⁻¹(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D,g::MaxwellGrid)::Array{SHM3,3}
    tree = KDTree(shapes)
    ε_sm_inv = copy(reshape( [SHM3(inv(εₛ(shapes,tree,g.x[i],g.y[j],g.δx,g.δy))) for i=1:g.Nx,j=1:g.Ny], (g.Nx,g.Ny,1)) )
end

function εₘₐₓ(shapes::AbstractVector{T} where T <: GeometryPrimitives.Shape{2,4,D} where D)
    maximum(vec([shapes[i].data[j,j] for j=1:3,i=1:size(shapes)[1]]))
end

ñₘₐₓ(ε⁻¹::AbstractArray{SHM3,3})::Float64 = √(maximum(3 ./ tr.(ε⁻¹)))
nₘₐₓ(ε::AbstractArray{SHM3,3})::Float64 = √(maximum(reinterpret(Float64,ε)))


text:
@testset "test key with value" begin
    d = Dict("foo" => "bar")
    @test_throws ErrorException PSI.find_key_with_value(d, "fake")
end

@testset "Test ProgressMeter" begin
    @test !PSI.progress_meter_enabled()
end

@testset "Axis Array to DataFrame" begin
    # The axis_array_to_dataframe test the use of the `to_matrix` and `get_column_names` methods
    one = PSI.DenseAxisArray{Float64}(undef, 1:2)
    fill!(one, 1.0)
    mock_key = PSI.VariableKey(ActivePowerVariable, ThermalStandard)
    one_df = PSI.axis_array_to_dataframe(one, mock_key)
    test_df = DataFrames.DataFrame(PSI.encode_key(mock_key) => [1.0, 1.0])
    @test one_df == test_df

    two = PSI.DenseAxisArray{Float64}(undef, [:a], 1:2)
    fill!(two, 1.0)
    two_df = PSI.axis_array_to_dataframe(two, mock_key)
    test_df = DataFrames.DataFrame(:a => [1.0, 1.0])
    @test two_df == test_df

    three = PSI.DenseAxisArray{Float64}(undef, [:a], 1:2, 1:3)
    @test_throws ErrorException PSI.axis_array_to_dataframe(three, mock_key)

    four = PSI.DenseAxisArray{Float64}(undef, [:a], 1:2, 1:3, 1:5)
    @test_throws ErrorException PSI.axis_array_to_dataframe(four, mock_key)

    sparse_num =
        JuMP.Containers.@container([i = 1:10, j = (i + 1):10, t = 1:24], 0.0 + i + j + t)
    @test_throws MethodError PSI.axis_array_to_dataframe(sparse_num, mock_key)

    i_num = 1:10
    j_vals = Dict("$i" => string.((i + 1):11) for i in i_num)
    sparse_valid =
        JuMP.Containers.@container([i = string.(i_num), j = j_vals[i], t = 1:24], rand())
    df = PSI.axis_array_to_dataframe(sparse_valid, mock_key)
    @test size(df) == (24, 55)
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ProfileViewGtk

using Gtk.ShortNames, GtkUtilities, Colors, FileIO
import Cairo
if VERSION < v"0.4.0-dev+3275"
    using Base.Graphics
else
    using Graphics
end

type ZoomCanvas
    bb::BoundingBox  # in user-coordinates
    c::Canvas
end

function __init__()
    eval(Expr(:import, :ProfileView))
end

function view(data = Profile.fetch(); lidict=nothing, kwargs...)
    bt, uip, counts, lidict, lkup = ProfileView.prepare_data(data, lidict)
    # Display in a window
    c = @Canvas()
    setproperty!(c, :expand, true)
    f = @Frame(c)
    tb = @Toolbar()
    bx = @Box(:v)
    push!(bx, tb)
    push!(bx, f)
    tb_open = @ToolButton("gtk-open")
    tb_save_as = @ToolButton("gtk-save-as")
    push!(tb, tb_open)
    push!(tb, tb_save_as)
    signal_connect(open_cb, tb_open, "clicked", Void, (), false, (c,kwargs))
    signal_connect(save_as_cb, tb_save_as, "clicked", Void, (), false, (c,data,lidict,kwargs))
    win = @Window(bx, "Profile")
    if data != nothing && !isempty(data)
        viewprof(c, bt, uip, counts, lidict, lkup; kwargs...)
    end
    showall(win)
end

function viewprof(c, bt, uip, counts, lidict, lkup; C = false, colorgc = true, fontsize = 12, combine = true, pruned=[])
    img, lidict, imgtags = ProfileView.prepare_image(bt, uip, counts, lidict, lkup, C, colorgc, combine, pruned)
    img24 = UInt32[convert(UInt32, convert(RGB24, img[i,j])) for i = 1:size(img,1), j = size(img,2):-1:1]'
    surf = Cairo.CairoRGBSurface(img24)
    imw = size(img24,2)
    imh = size(img24,1)
    panzoom(c, (0,imw), (0,imh))
    panzoom_mouse(c)
    panzoom_key(c)
    lasttextbb = BoundingBox(1,0,1,0)
    standard_motion = function (c, event)
        # Repair image from ovewritten text
        ctx = getgc(c)
        w = width(c)
        if width(lasttextbb) > 0
            h = height(c)
            xview, yview = guidata[c, :xview], guidata[c, :yview]
            set_coords(ctx, xview, yview)
            rectangle(ctx, lasttextbb)
            set_source(ctx, surf)
            p = Cairo.get_source(ctx)
            Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
            fill(ctx)
        end
        # Write the info
        xd, yd = event.x, event.y
        xu, yu = device_to_user(ctx, xd, yd)
        tag = gettag(xu, yu)
        if tag != ProfileView.TAGNONE
            li = lidict[tag.ip]
            if VERSION < v"0.5.0-dev+4192"
                str = string(basename(string(li.file)), ", ", li.func, ": line ", li.line)
            else
                str = ""
                for l in li
                    str = string(str, string(basename(string(l.file)), ", ", l.func, ": line ", l.line), "; ")
                end
            end
            set_source(ctx, ProfileView.fontcolor)
            Cairo.set_font_face(ctx, "sans-serif $(fontsize)px")
            lasttextbb = deform(Cairo.text(ctx, xu, yu, str, halign = xd < w/3 ? "left" : xd < 2w/3 ? "center" : "right"), -2, 2, -2, 2)
        end
        reveal(c)
    end
    c.mouse.motion = standard_motion
    draw(c) do widget
        ctx = getgc(widget)
        w = width(widget)
        h = height(widget)
        xview, yview = guidata[c, :xview], guidata[c, :yview]
        set_coords(ctx, xview, yview)
        rectangle(ctx, xview.min, yview.min, width(xview), width(yview))
        set_source(ctx, surf)
        p = Cairo.get_source(ctx)
        Cairo.pattern_set_filter(p, Cairo.FILTER_NEAREST)
        fill(ctx)
    end
    # From a given position, find the underlying tag
    function gettag(xu, yu)
        x = ceil(Int, xu)
        y = ceil(Int, yu)
        Y = size(imgtags, 2)
        x = max(1, min(x, size(imgtags, 1)))
        y = max(1, min(y, Y))
        imgtags[x,Y-y+1]
    end
    # Hover over a block and see the source line
    # Right-click prints the full path, function, and line to the console
    c.mouse.button3press = function (c, event)
        ctx = getgc(c)
        xd, yd = event.x, event.y
        xu, yu = device_to_user(ctx, xd, yd)
        tag = gettag(xu, yu)
        if tag != ProfileView.TAGNONE
            li = lidict[tag.ip]
            if VERSION < v"0.5.0-dev+4192"
                println(li.file, ", ", li.func, ": line ", li.line)
            else
                firstline = true
                for l in li
                    if !firstline
                        print("  ")
                    else
                        firstline = false
                    end
                    println(l.file, ", ", l.func, ": line ", l.line)
                end
            end
        end
    end
    nothing
end

function open_cb(::Ptr, settings::Tuple)
    c, kwargs = settings
    selection = open_dialog("Load profile data", toplevel(c), ("*.jlprof","*"))
    isempty(selection) && return nothing
    data, lidict = load(File(format"JLD", selection), "li", "lidict")
    bt, uip, counts, lidict, lkup = ProfileView.prepare_data(data, lidict)
    viewprof(c, bt, uip, counts, lidict, lkup; kwargs...)
    nothing
end

function save_as_cb(::Ptr, profdata::Tuple)
    c, data, lidict = profdata
    selection = save_dialog("Save profile data as JLD file", toplevel(c), ("*.jlprof",))
    isempty(selection) && return nothing
    FileIO.save(File(format"JLD", selection), "li", data, "lidict", lidict)
    nothing
end

end


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ 754d60ac-4cb0-443a-af2a-25235e544a31
begin
  using Pkg
  Pkg.activate(@__DIR__)
  Pkg.instantiate()
end

# ╔═╡ 0256dc40-e8b0-40a8-a6bf-47194d7e8e12
using Unitful

# ╔═╡ 0f76a79f-8675-4ec1-ba1f-6363f43ec697
using Measurements

# ╔═╡ 4dbb0300-d546-48dd-9bca-7c69b794f8ce
md"# Basic demonstration of Julia package composability"

# ╔═╡ eb23f2c6-61f4-4a82-b795-033f6f2a0674
md"Notebook from [HelloJulia.jl](https://github.com/ablaom/HelloJulia.jl)"

# ╔═╡ ec6ad8e5-c854-41db-935f-ddf6a6bfbbdd
md"Instantiate package environment:"

# ╔═╡ b35b6a58-0335-4dab-8af5-148d95ea2900
md"The Unitiful package allows you to **bind physical units** to numerical data:"

# ╔═╡ a2808946-3e6a-4c31-828a-4cc485149963
a = 5.0u"m^2/s^2"

# ╔═╡ 92fe3d8a-f4c0-4de4-9e54-3f7b3ca87ecb
md"""
The using Measurements package allows you to **propogate uncertainties**
in numerical computations:
"""

# ╔═╡ bf15f7d4-5ef4-4184-95ea-2955abd45275
b = 5.0 ± 1.2 # or measurement(5.0, 1.2)

# ╔═╡ 4c307469-e7bd-4041-b1b3-dd47e367ba54
sqrt(b)

# ╔═╡ 9e262366-a19c-4d4c-8d7e-8e449afd1c48
md"""
The Unitful and Measurements packages are blissfully ignorant of one
another. That, is neither package is a dependency of the other. And yet the
following "just works", as if by magic:
"""

# ╔═╡ cddcd722-2cd4-4d95-a949-7f3bd292fe31
c = 5.0u"m^2/s^2" ± 0.1u"m^2/s^2"

# ╔═╡ 76009ff9-a8ba-4cba-8514-99938a2932db
sqrt(c)

# ╔═╡ 135dac9b-0bd9-4e1d-a0de-1721c1bc2df2
md"""
---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
"""

# ╔═╡ Cell order:
# ╟─4dbb0300-d546-48dd-9bca-7c69b794f8ce
# ╟─eb23f2c6-61f4-4a82-b795-033f6f2a0674
# ╟─ec6ad8e5-c854-41db-935f-ddf6a6bfbbdd
# ╠═754d60ac-4cb0-443a-af2a-25235e544a31
# ╟─b35b6a58-0335-4dab-8af5-148d95ea2900
# ╠═0256dc40-e8b0-40a8-a6bf-47194d7e8e12
# ╠═a2808946-3e6a-4c31-828a-4cc485149963
# ╟─92fe3d8a-f4c0-4de4-9e54-3f7b3ca87ecb
# ╠═0f76a79f-8675-4ec1-ba1f-6363f43ec697
# ╠═bf15f7d4-5ef4-4184-95ea-2955abd45275
# ╠═4c307469-e7bd-4041-b1b3-dd47e367ba54
# ╟─9e262366-a19c-4d4c-8d7e-8e449afd1c48
# ╠═cddcd722-2cd4-4d95-a949-7f3bd292fe31
# ╠═76009ff9-a8ba-4cba-8514-99938a2932db
# ╟─135dac9b-0bd9-4e1d-a0de-1721c1bc2df2


text:
using jInv.Mesh;
using FactoredEikonalFastMarching;
using PyPlot
close("all")

include("runAccuracyExperiment.jl");
include("getAnalyticalMediums.jl");
include("getWorkUnit.jl");

function runPaperExperiments()
numOfRefinements = 4;

println("******************** 2D experiments ***************************");

I = [4,8];
n = zeros(Int64,2);

h0 = [0.05,0.05];
WU = zeros(numOfRefinements)
for k = 1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	WU[k] = GetWorkunitForLoop(n,h);
end

println("******************** Gradient Slowness ***************************");
println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	tc="testcase1";
	println("*******************************************************************************")
	println("Testcase 1")
    println("*******************************************************************************")
	(kappaSquared,src,T_exact) = getAnalyticalConstGrad2D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
end

println("******************** Gradient Velocity ***************************");

println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	tc="testcase2";
	println("*******************************************************************************")
	println("*******************************************************************************")
	println("Testcase 2")
    println("*******************************************************************************")
	(kappaSquared,src,T_exact) = getAnalyticalConstGradInv2D(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
end

println("******************** Gaussian factor ***************************");

println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
h0 = [0.05,0.05];
for k=1:numOfRefinements
	h = (0.5^k)*h0;
	n[1] = round(Int64,(I[1]/h[1])+1);
	n[2] = round(Int64,(I[2]/h[2])+1);
	tc="testcase3";
	println("*******************************************************************************")
    println("*******************************************************************************")
	println("Testcase 3")
    println("*******************************************************************************")
	(kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
	runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
end

#=println("*******************************************************************************")=#
# println("*******************************************************************************")
# println("************************* 3-D EXPERIMENTS *************************************")
# println("*******************************************************************************")
# println("*******************************************************************************")

# ##########################################################################################
# ##########################################################################################
# # 3D experiments:
# ##########################################################################################

# I = [1.6,1.6,0.8];
# n = zeros(Int64,3);

# numOfRefinements = 4;

# h0 = [0.1,0.1,0.1];
# WU = zeros(numOfRefinements)
# for k = 1:numOfRefinements
	# h = (0.5^k)*h0;
	# n[1] = round(Int64,(I[1]/h[1])+1);
	# n[2] = round(Int64,(I[2]/h[2])+1);
	# n[3] = round(Int64,(I[3]/h[3])+1);
	# WU[k] = GetWorkunitForLoop(n,h);
# end

# println("******************** Gradient Slowness ***************************");
# println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

# for k=1:numOfRefinements
	# h = (0.5^k)*h0;
	# n[1] = round(Int64,(I[1]/h[1])+1);
	# n[2] = round(Int64,(I[2]/h[2])+1);
	# n[3] = round(Int64,(I[3]/h[3])+1);
	# (kappaSquared,src,T_exact) = getAnalyticalConstGrad3D(n,h);
	# runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
# end

# println("******************** Gradient Velocity ***************************");
# println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
# for k=1:numOfRefinements
	# h = (0.5^k)*h0;
	# n[1] = round(Int64,(I[1]/h[1])+1);
	# n[2] = round(Int64,(I[2]/h[2])+1);
	# n[3] = round(Int64,(I[3]/h[3])+1);
	# (kappaSquared,src,T_exact) = getAnalyticalConstGradInv3D(n,h);
	# runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
# end

# println("******************** Gaussian factor ***************************");
# println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
# for k=1:numOfRefinements
	# h = (0.5^k)*h0;
	# n[1] = round(Int64,(I[1]/h[1])+1);
	# n[2] = round(Int64,(I[2]/h[2])+1);
	# n[3] = round(Int64,(I[3]/h[3])+1);
	# (kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
	# runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k]);
#=end=#
return;
end
############################################################################################

runPaperExperiments();


text:
# test of solution stability
using Distributed;
addprocs(20);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();

N = 5;
dataList = Dict();

caseList = [13,33,123];
Δt = 0.25;
iterMax = 20;
trainNo = 20;
NN = 5000;

for ci in 1:length(caseList)
    TList = [24,36,48,72,96];
    pathListDRaw = load("pathHist_5000.jld");
    pathDictA = pathListDRaw["pathDict"][ci];
    dataList[ci] = Dict();

    for T in TList
        τ = Int64(1/6*T);
        for j in procs()
            remotecall_fetch(readInData,j,ci,caseList,T,τ);
        end
        pathDict = pathDictA[T];
        dataList[ci][T] = Dict();

        # train trainNo of strategies and obtain the 95% PI
        for tn in 1:trainNo
            # train the stochastic programming strategy
            cutDictPG = preGen(T, Δt, N, iterMax);
            #cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, true, false, 30, 30, cutDictPG);
            cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(T, Δt, N, false, false,
                max(Int64(round(500/N)),20), max(Int64(round(500/N)),20),cutDictPG);
            for j in procs()
                remotecall_fetch(cutIni,j,cutDict);
            end

            # test against 5000 samples
            solSDDP, LBSDDP, costSDDP = exeForward(T, Δt, NN, false, pathDict);
            listSDDP = [costSDDP[i] for i in 1:NN];
            meanSDDP = mean(listSDDP);
            sigmaSDDP = std(listSDDP);
            println(round(meanSDDP,digits = 2)," ",round(meanSDDP - 1.96*sigmaSDDP,digits = 2)," ",round(meanSDDP + 1.96*sigmaSDDP,digits = 2));
            dataList[ci][T][tn] = [LBSDDP, costSDDP, listSDDP, meanSDDP, sigmaSDDP];
            save("stabilityResults.jld","data",dataList);
        end
    end
end


text:
JSON.lower(pa::PostAggregator) = non_nothing_dict(pa)

"""
    Arithmetic(name::String, fn::String, fields::Vector{<:PostAggregator}; ordering=nothing)

Apply the `fn` aggregation to the `fields`.

ordering should be a String if provided.
"""
struct Arithmetic <: PostAggregator
    type::String
    name::String
    fn::String
    fields::Vector{<:PostAggregator}
    ordering
    function Arithmetic(name, fn, fields; ordering=nothing)
        fn = string(fn)
        fn ∈ ["+", "-", "*", "/", "quotient"] || error("Unknown arithmetic function")
        ordering === nothing || ordering == "numericFirst" || error("Invalid ordering")
        new("arithmetic", name, fn, fields, ordering)
    end
end

"""
    FieldAccess(name::String, aggregator::String)
    FieldAccess(aggregator::String)

Access the `aggregator` and return as `name`.
"""
struct FieldAccess <: PostAggregator
    type::String
    name::String
    fieldName::String
    FieldAccess(name, aggregator) = new("fieldAccess", name, aggregator)
end
FieldAccess(aggregator) = FieldAccess(aggregator, aggregator)

"""
    FinalizingFieldAccess(name::String, aggregator::String)
    FinalizingFieldAccess(aggregator::String)

Access the `aggregator` and return as `name`.
"""
struct FinalizingFieldAccess <: PostAggregator
    type::String
    name::String
    fieldName::String
    FinalizingFieldAccess(name, aggregator) = new("finalizingFieldAccess", name, aggregator)
end
FinalizingFieldAccess(aggregator) = FinalizingFieldAccess(aggregator, aggregator)

"""
    ConstantPA(name::String, value::Real)

Return a constant value as `name`.
"""
struct ConstantPA <: PostAggregator
    type::String
    name::String
    value::Real
    ConstantPA(name, value) = new("constant", name, value)
end

"""
    Greatest(name::String, fields::Vector{<:PostAggregator}, dtype::String)

Find the largest value per row across the `fields` of type `dtype` (double,
long) and return as `name`.
"""
struct Greatest <: PostAggregator
    type::String
    name::String
    fields::Vector{<:PostAggregator}
    function Greatest(name, fields, dtype)
        dtype = lowercase(dtype)
        dtype ∈ ["double", "long"] || error("Invalid data type")
        new(dtype * "Greatest", name, fields)
    end
end

"""
    Least(name::String, fields::Vector{<:PostAggregator}, dtype::String)

Find the smallest value per row across the `fields` of type `dtype` (double,
long) and return as `name`.
"""
struct Least <: PostAggregator
    type::String
    name::String
    fields::Vector{<:PostAggregator}
    function Least(name, fields, dtype)
        dtype = lowercase(dtype)
        dtype ∈ ["double", "long"] || error("Invalid data type")
        new(dtype * "Least", name, fields)
    end
end

"""
    JavaScriptPA(name::String, fieldNames::Vector{String}, jsfunction::String)

Post aggregate the `fieldNames` aggregators using a custom `jsfunction`
JavaScript function.
"""
struct JavaScriptPA <: PostAggregator
    type::String
    name::String
    fieldNames::Vector{String}
    jsfunction::String
    JavaScriptPA(name, fields, jsfunction) = new("javascript", name, fields, jsfunction)
end

JSON.lower(pa::JavaScriptPA) = Dict("type" => pa.type, "name" => pa.name, "fieldNames" => pa.fieldNames, "function" => pa.jsfunction)

"""
    HyperUniqueCardinality(name::String, field::String)
    HyperUniqueCardinality(field::String)

Post aggregator of a hyper unique aggregation.
"""
struct HyperUniqueCardinality <: PostAggregator
    type::String
    name::String
    fieldName::String
    HyperUniqueCardinality(name, field) = new("hyperUniqueCardinality", name, field)
end
HyperUniqueCardinality(field) = HyperUniqueCardinality(field, field)


text:
X = collect(Float64, reshape(1:40, 10, 4))
x = rand(10) * 10

D = DataFrame(A=rand(10), B=collect(1:10), C=[hex(x) for x in 11:20])
D_NA = deepcopy(D)
D_NA[1, :A] = NA

@testset "Array" begin
    # Rescale Vector
    xx = deepcopy(x)
    mu, sigma = standardize!(xx)
    @test mu ≈ mean(x)
    @test sigma ≈ std(x)
    @test abs(mean(xx)) <= 10e-10
    @test std(xx) ≈ 1

    xx = deepcopy(x)
    mu, sigma = standardize!(xx, mu, sigma)
    @test abs(mean(xx)) <= 10e-10
    @test std(xx) ≈ 1

    xx = deepcopy(x)
    mu, sigma = standardize!(xx, mu, sigma, obsdim=1)
    @test abs(mean(xx)) <= 10e-10
    @test std(xx) ≈ 1

    xx = deepcopy(x)
    mu = deepcopy(x) .- 1
    sigma = ones(x)
    mu, sigma = standardize!(xx, mu, sigma, obsdim=1)
    @test mean(xx) ≈ 1

    # Rescale Matrix
    XX = deepcopy(X)
    standardize!(XX)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    XX = deepcopy(X)
    standardize!(XX, obsdim=2)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    XX = deepcopy(X)
    standardize!(XX, obsdim=1)
    @test abs(sum(mean(XX, 1))) <= 10e-10
    @test vec(std(XX, 1)) ≈ ones(size(X, 2))

    XX = deepcopy(X)
    mu = vec(mean(XX, 1))
    sigma = vec(std(XX, 1))
    standardize!(XX, mu, sigma, obsdim=1)
    @test abs(sum(mean(XX, 1))) <= 10e-10

    XX = deepcopy(X)
    mu = vec(mean(XX, 2))
    sigma = vec(std(XX, 2))
    standardize!(XX, mu, sigma, obsdim=2)
    @test abs(sum(mean(XX, 2))) <= 10e-10

    XX = deepcopy(X)
    flt = [1,2]
    standardize!(XX, obsdim=1, operate_on=flt)
    @test abs(sum(mean(XX[:,flt], 1))) <= 10e-10
    @test vec(std(XX[:,flt], 1)) ≈ ones(2)
    @test all(X[:,[3,4]] .== XX[:,[3,4]])

    XX = deepcopy(X)
    flt = [2,8]
    mu = vec(mean(XX, 2))
    sigma = vec(std(XX, 2))
    standardize!(XX, mu[flt], sigma[flt], obsdim=2, operate_on=flt)
    @test abs(sum(mean(XX[flt,:], 2))) <= 10e-10

    scaler = StandardScaler(X)
    XX = transform(X, scaler)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    scaler = fit(StandardScaler, X)
    XX = transform(X, scaler)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    Xi = round.(Int, X)
    XX = transform(Xi, scaler)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    XX, scaler = fit_transform(StandardScaler, X)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    XX = deepcopy(X)
    scaler = fit_transform!(StandardScaler, XX)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    scaler = fit(StandardScaler, X, obsdim=2)
    XX = transform(X, scaler)
    @test abs(sum(mean(XX, 2))) <= 10e-10
    @test std(XX, 2) ≈ ones(size(X, 1))

    scaler = fit(StandardScaler, X, obsdim=1)
    XX = transform(X, scaler)
    @test abs(sum(mean(XX, 1))) <= 10e-10
    @test vec(std(XX, 1)) ≈ ones(size(X, 2))

    flt = [1,4]
    scaler = fit(StandardScaler, X, obsdim=1, operate_on=flt)
    XX = transform(X, scaler)
    xx = transform(vec(X[1,:]), scaler)
    @test abs(sum(mean(XX[:,flt], 1))) <= 10e-10
    @test vec(std(XX[:,flt], 1)) ≈ ones(size(X[:,flt], 2))
    @test all(xx .== XX[1,:])

    XX = deepcopy(X)
    xx = vec(X[1,:])
    flt = [1,4]
    scaler = fit(StandardScaler, X, obsdim=1, operate_on=flt)
    transform!(XX, scaler)
    transform!(xx, scaler)
    @test abs(sum(mean(XX[:,flt], 1))) <= 10e-10
    @test vec(std(XX[:,flt], 1)) ≈ ones(size(X[:,flt], 2))
    @test all(xx .== XX[1,:])
end

@testset "DataFrame" begin
    DD = deepcopy(D)
    mu, sigma = standardize!(DD)
    @test abs(sum([mean(DD[colname]) for colname in names(DD)[1:2]])) <= 10e-10
    @test mean([std(DD[colname]) for colname in names(DD)[1:2]]) - 1 <= 10e-10

    DD = deepcopy(D)
    mu, sigma = standardize!(DD, operate_on=[:A,:B,:C])
    @test abs(sum([mean(DD[colname]) for colname in names(DD)[1:2]])) <= 10e-10
    @test mean([std(DD[colname]) for colname in names(DD)[1:2]]) - 1 <= 10e-10

    DD = deepcopy(D_NA)
    m, s = standardize!(DD, [:A,:B,:C])
    @test all(DD[2:end, :A] .== D_NA[2:end, :A])
    @test all(DD[:C] .== D_NA[:C])
    @test mean(DD[:B]) ≈ 0.0

    DD = deepcopy(D_NA)
    m = 0.0
    s = 1.0
    m, s = standardize!(DD, m, s, :A)
    @test all(DD[2:end, :A] .== D_NA[2:end, :A])
    @test all(DD[:C] .== D_NA[:C])

    DD = deepcopy(D)
    mu, sigma = standardize!(DD, mu, sigma, operate_on=[:A,:B])
    @test abs(sum([mean(DD[colname]) for colname in names(DD)[1:2]])) <= 10e-10
    @test mean([std(DD[colname]) for colname in names(DD)[1:2]]) - 1 <= 10e-10

    # skip columns that contain NA values
    DD = deepcopy(D_NA)
    mu, sigma = standardize!(DD)
    @test isna(DD[1, :A])
    @test all(DD[2:end, :A] .== D_NA[2:end, :A])
    @test abs(mean(DD[:B])) < 10e-10
    @test abs(std(DD[:B])) - 1 < 10e-10

    scaler = StandardScaler(D)
    DD = transform(D, scaler)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test mean(DD[:B]) <= 10e-10
    @test std(DD[:B]) - 1  <= 10e-10
    @test all(DD[:C] .== D[:C])

    scaler = fit(StandardScaler, D)
    DD = transform(D, scaler)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test mean(DD[:B]) <= 10e-10
    @test std(DD[:B]) - 1  <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD, scaler = fit_transform(StandardScaler, D)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test mean(DD[:B]) <= 10e-10
    @test std(DD[:B]) - 1  <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    scaler = fit_transform!(StandardScaler, DD)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test mean(DD[:B]) <= 10e-10
    @test std(DD[:B]) - 1  <= 10e-10
    @test all(DD[:C] .== D[:C])

    colnames = [:A, :B]
    offset = Float64[mean(D[colname]) for colname in colnames]
    scale = Float64[std(D[colname]) for colname in colnames]
    scaler = StandardScaler(D, offset, scale)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test mean(DD[:B]) <= 10e-10
    @test std(DD[:B]) - 1  <= 10e-10
    @test all(DD[:C] .== D[:C])

    scaler = fit(StandardScaler, D, operate_on=[:A, :C])
    DD = transform(D, scaler)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test all(DD[:B] .== D[:B])
    @test all(DD[:C] .== D[:C])
    @test mean(D[:A]) != mean(DD[:A])

    DD = deepcopy(D)
    scaler = fit(StandardScaler, DD, operate_on=[:A, :C])
    transform!(DD, scaler)
    @test mean(DD[:A]) <= 10e-10
    @test std(DD[:A]) - 1  <= 10e-10
    @test all(DD[:B] .== D[:B])
    @test all(DD[:C] .== D[:C])
    @test mean(D[:A]) != mean(DD[:A])
end