text:
export DecisionTreeClassifier, DecisionTreeRegressor

using MLJScientificTypes

using CategoricalArrays

import DecisionTree

## DESCRIPTIONS

const DTC_DESCR = "Decision Tree Classifier."
const DTR_DESCR = "Decision Tree Regressor."

## CLASSIFIER

struct TreePrinter{T}
    tree::T
end
(c::TreePrinter)(depth) = DecisionTree.print_tree(c.tree, depth)
(c::TreePrinter)() = DecisionTree.print_tree(c.tree, 5)

Base.show(stream::IO, c::TreePrinter) =
    print(stream, "TreePrinter object (call with display depth)")

"""
DecisionTreeClassifer(; kwargs...)

A variation on the CART decision tree classifier from [https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md](https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md).

Inputs are tables with ordinal columns. That is, the element scitype
of each column can be `Continuous`, `Count` or `OrderedFactor`.

Instead of predicting the mode class at each leaf, a UnivariateFinite
distribution is fit to the leaf training classes, with smoothing
controlled by an additional hyperparameter `pdf_smoothing`: If `n` is
the number of observed classes, then each class probability is
replaced by `pdf_smoothing/n`, if it falls below that ratio, and the
resulting vector of probabilities is renormalized. Smoothing is only
applied to classes actually observed in training. Unseen classes
retain zero-probability predictions.

To visualize the fitted tree in the REPL, set `verbosity=2` when
fitting, or call `report(mach).print_tree(display_depth)` where `mach`
is the fitted machine, and `display_depth` the desired
depth. Interpretting the results will require a knowledge of the
internal integer encodings of classes, which are given in
`fitted_params(mach)` (which also stores the raw learned tree object
from the DecisionTree.jl algorithm).

For post-fit pruning, set `post-prune=true` and set
`min_purity_threshold` appropriately. Other hyperparameters as per
package documentation cited above.

"""
@mlj_model mutable struct DecisionTreeClassifier <: MLJBase.Probabilistic
    pruning_purity::Float64         = 1.0::(_ ≤ 1)
    max_depth::Int                  = (-)(1)::(_ ≥ -1)
    min_samples_leaf::Int           = 1::(_ ≥ 0)
    min_samples_split::Int          = 2::(_ ≥ 2)
    min_purity_increase::Float64    = 0.0::(_ ≥ 0)
    n_subfeatures::Int              = 0::(_ ≥ 0)
    display_depth::Int              = 5::(_ ≥ 1)
    post_prune::Bool                = false
    merge_purity_threshold::Float64 = 0.9::(0 ≤ _ ≤ 1)
    pdf_smoothing::Float64          = 0.05::(0 ≤ _ ≤ 1)
end

#> A required `fit` method returns `fitresult, cache, report`. (Return
#> `cache=nothing` unless you are overloading `update`)
function MLJBase.fit(model::DecisionTreeClassifier, verbosity::Int, X, y)
    Xmatrix = MLJBase.matrix(X)
    yplain  = MLJBase.int(y)
    classes_seen = filter(in(unique(y)), MLJBase.classes(y[1]))
    integers_seen = MLJBase.int(classes_seen) #unique(yplain)

    tree = DecisionTree.build_tree(yplain, Xmatrix,
                                   model.n_subfeatures,
                                   model.max_depth,
                                   model.min_samples_leaf,
                                   model.min_samples_split,
                                   model.min_purity_increase)
    if model.post_prune
        tree = DecisionTree.prune_tree(tree, model.merge_purity_threshold)
    end

    verbosity < 2 || DecisionTree.print_tree(tree, model.display_depth)

    fitresult = (tree, classes_seen, integers_seen)

    #> return package-specific statistics (eg, feature rankings,
    #> internal estimates of generalization error) in `report`, which
    #> should be a named tuple with the same type every call (can have
    #> empty values):

    cache = nothing
    report = (classes_seen=classes_seen,
              print_tree=TreePrinter(tree))

    return fitresult, cache, report
end

function get_encoding(classes_seen)
    a_cat_element = classes_seen[1]
    return Dict(c => MLJBase.int(c) for c in MLJBase.classes(a_cat_element))
end

MLJBase.fitted_params(::DecisionTreeClassifier, fitresult) =
    (tree_or_leaf=fitresult[1], encoding=get_encoding(fitresult[2]))

function smooth(prob_vector, smoothing)
    threshold = smoothing/length(prob_vector)
    smoothed_vector = map(prob_vector) do p
        p < threshold ? threshold : p
    end
    smoothed_vector = smoothed_vector/sum(smoothed_vector)
    return smoothed_vector
end

function MLJBase.predict(model::DecisionTreeClassifier
                     , fitresult
                     , Xnew)
    Xmatrix = MLJBase.matrix(Xnew)

    tree, classes_seen, integers_seen = fitresult

    y_probabilities =
        DecisionTree.apply_tree_proba(tree, Xmatrix, integers_seen)

    return [MLJBase.UnivariateFinite(classes_seen,
                                     smooth(y_probabilities[i,:],
                                            model.pdf_smoothing))
            for i in 1:size(y_probabilities, 1)]
end

## REGRESSOR

"""
DecisionTreeRegressor(; kwargs...)

CART decision tree classifier from
[https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md](https://github.com/bensadeghi/DecisionTree.jl/blob/master/README.md). Predictions
are Deterministic.

Inputs are tables with ordinal columns. That is, the element scitype
of each column can be `Continuous`, `Count` or `OrderedFactor`.

For post-fit pruning, set `post-prune=true` and set
`pruning_purity_threshold` appropriately. Other hyperparameters as per
package documentation cited above.

"""
@mlj_model mutable struct DecisionTreeRegressor <: MLJBase.Deterministic
    pruning_purity_threshold::Float64 = 0.0::(0 ≤ _ ≤ 1)
    max_depth::Int					  = (-)(1)::(_ ≥ -1)
    min_samples_leaf::Int			  = 5::(_ ≥ 0)
    min_samples_split::Int			  = 2::(_ ≥ 2)
    min_purity_increase::Float64	  = 0.0::(_ ≥ 0)
    n_subfeatures::Int				  = 0::(_ ≥ 0)
    post_prune::Bool				  = false
end

function MLJBase.fit(model::DecisionTreeRegressor, verbosity::Int, X, y)
    Xmatrix   = MLJBase.matrix(X)
    fitresult = DecisionTree.build_tree(float.(y), Xmatrix
       , model.n_subfeatures
       , model.max_depth
       , model.min_samples_leaf
       , model.min_samples_split
       , model.min_purity_increase)

    if model.post_prune
        fitresult = DecisionTree.prune_tree(fitresult,
                                            model.pruning_purity_threshold)
    end
    cache = nothing
    report = nothing

    return fitresult, cache, report
end

MLJBase.fitted_params(::DecisionTreeRegressor, fitresult) =
    (tree_or_leaf = fitresult,)

function MLJBase.predict(model::DecisionTreeRegressor
                     , fitresult
                     , Xnew)
    Xmatrix = MLJBase.matrix(Xnew)
    return DecisionTree.apply_tree(fitresult,Xmatrix)
end

##
## METADATA
##

metadata_pkg.((DecisionTreeClassifier, DecisionTreeRegressor),
              name="DecisionTree",
              uuid="7806a523-6efd-50cb-b5f6-3fa6f1930dbb",
              url="https://github.com/bensadeghi/DecisionTree.jl",
              julia=true,
              license="MIT",
              is_wrapper=false)

metadata_model(DecisionTreeClassifier,
               input=MLJBase.Table(Continuous, Count, OrderedFactor),
               target=AbstractVector{<:MLJBase.Finite},
               weights=false,
               descr=DTC_DESCR)

metadata_model(DecisionTreeRegressor,
               input=MLJBase.Table(Continuous, Count, OrderedFactor),
               target=AbstractVector{MLJBase.Continuous},
               weights=false,
               descr=DTR_DESCR)

ext: jl


text:
@testset "Snapshot reading tests" begin

# Bad data filter procedures.
uploadSnapshot( mpSim, "snapshot/test foto 1", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )
@test ( mpSim.orgSize == 50 ) && !mpSim.isVirgin

uploadSnapshot( mpSim, "snapshot/test foto 2", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )
@test ( mpSim.orgSize == 45 ) && !mpSim.isVirgin

uploadSnapshot( mpSim, "snapshot/test foto 3", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )
@test ( mpSim.orgSize == 45 ) && !mpSim.isVirgin

uploadSnapshot( mpSim, "snapshot/test foto 4", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )
@test ( mpSim.orgSize == 50 ) && !mpSim.isVirgin

uploadSnapshot( mpSim, "snapshot/test foto 5", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H"), refDate=Date( 2021, 1, 1 ) )
@test ( mpSim.orgSize == 50 ) && !mpSim.isVirgin

# Generate data procedure.
clearSimulationAttributes!( mpSim )
clearSimulationBaseNodes!( mpSim )
uploadSnapshot( mpSim, "snapshot/test foto 1", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H"), generateSimData=true )
@test ( mpSim.orgSize == 50 ) && !mpSim.isVirgin
@test length( mpSim.attributeList ) == 3
@test length( mpSim.baseNodeList ) == 5

uploadSnapshot( mpSim, "snapshot/test foto 1", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", 0), generateSimData=true )
@test ( mpSim.orgSize == 50 ) && !mpSim.isVirgin
nodes = get.( Ref(mpSim.baseNodeList), ["A", "B", "C", "D", "E"], nothing )
@test all( length.( getfield.( nodes, :inNodeSince ) ) .== [17, 13, 10, 5, 5] )

# Attrition generation procedure.
setNodeAttritionScheme!( mpSim.baseNodeList["A"], "A/B" )
setNodeAttritionScheme!( mpSim.baseNodeList["B"], "A/B" )
setNodeAttritionScheme!( mpSim.baseNodeList["C"], "C/D" )
setNodeAttritionScheme!( mpSim.baseNodeList["D"], "C/D" )
setNodeAttritionScheme!( mpSim.baseNodeList["E"], "E" )

uploadSnapshot( mpSim, "snapshot/test foto 1", "Sheet1", ["B", "C", "D"],
    ("A", "G", "F", "E", "H") )

for node in keys( mpSim.baseNodeList )
    setNodeAttritionScheme!( mpSim.baseNodeList[node], "default" )
end  # for node in keys( mpSim.baseNodeList )

end  # @testset "Snapshot reading tests"


text:
module NeuralALU

include("layers.jl")

export NALU, NAC

end #module


text:
# ---
# title: 66. Plus One
# id: problem66
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/plus-one/description/>
# hidden: true
# ---
#
# Given a **non-empty** array of decimal digits representing a non-negative
# integer, increment one to the integer.
#
# The digits are stored such that the most significant digit is at the head of
# the list, and each element in the array contains a single digit.
#
# You may assume the integer does not contain any leading zero, except the
# number 0 itself.
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = [1,2,3]
#     Output: [1,2,4]
#     Explanation: The array represents the integer 123.
#
#
# **Example 2:**
#
#
#
#     Input: digits = [4,3,2,1]
#     Output: [4,3,2,2]
#     Explanation: The array represents the integer 4321.
#
#
# **Example 3:**
#
#
#
#     Input: digits = [0]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= digits.length <= 100`
#   * `0 <= digits[i] <= 9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DrWatson
@quickactivate :TimeProbeSeismic

close("all")
ee = (0, .8*25, .206*25, 0)

n, d, m, m0 = h5read(datadir("models", "overthrust_model.h5"), "n", "d", "m", "m0")
m0[:, 20:end] = imfilter(m0[:, 20:end] ,Kernel.gaussian(5));
n = Tuple(n)
d = Tuple(d)
vp_t = m'.^(-.5);
vp_0 = m0'.^(-.5);

inds = Dict(j=>i for (i,j)=enumerate([2^k for k=1:8]))

plt_dict = (cmap=:gist_ncar, vmin=1.5, vmax=6.0, extent=ee, aspect=:auto)

files = readdir(datadir("fwi_overthrust"))

dfts = [f for f=files if occursin("dft", f)]
probed = [f for f=files if occursin("ps", f)]

trinds = [200, 480, 680]
# True and initial model
figure(figsize=(12, 9))
subplot(2,1,1)
imshow(vp_t; plt_dict...)
# title("True")
xlabel("X (km)")
ylabel("Depth (km)")
vlines(x=[t*d[1]/1000 for t=trinds], colors=:k, ymin=0, ymax=ee[3])
subplot(2,1,2)
imshow(vp_0; plt_dict...)
# title("Initial")
xlabel("X (km)")
ylabel("Depth (km)")
tight_layout()
wsave(plotsdir("fwi_overthrust", "Init.png"), gcf())

# Load true
res = wload(datadir("fwi_overthrust", "fwi_std.bson"))
vp_std = reshape(res[:x], n)'.^(-.5)
ϕ_std = res[:ϕ_trace]

# DFT
ϕ_dft = Array{Any}(undef, 8)
vp_dft = Array{Any}(undef, 8)
figure(figsize=(12, 9))
subplot(3,3,1)
imshow(vp_std; plt_dict...)
title("FWI")
xlabel("X (km)")
ylabel("Depth (km)")
for f in dfts
    res = wload(datadir("fwi_overthrust", f))
    num = parse(Int, split(split(f, "_")[2], ".")[1][4:end])
    vp_dft[inds[num]] = reshape(res[:x], n)'.^(-.5)
    subplot(3, 3, inds[num]+1)
    imshow(vp_dft[inds[num]]; plt_dict...)
    title("OTDFT $(num)")
    xlabel("X (km)")
    ylabel("Depth (km)")
    ϕ_dft[inds[num]] = res[:ϕ_trace]
end
tight_layout()
wsave(plotsdir("fwi_overthrust", "DFT_fwi.png"), gcf())

# probed
ϕ_ps = Array{Any}(undef, 8)
vp_ps = Array{Any}(undef, 8)
figure(figsize=(12, 9))
subplot(3,3,1)
imshow(vp_std; plt_dict...)
title("Standard")
xlabel("X (km)")
ylabel("Depth (km)")
for f in probed
    res = wload(datadir("fwi_overthrust", f))
    num = parse(Int, split(split(f, "_")[2], ".")[1][3:end])
    vp_ps[inds[num]] = reshape(res[:x], n)'.^(-.5)
    subplot(3, 3, inds[num]+1)
    imshow(vp_ps[inds[num]]; plt_dict...)
    title("PFWI $(num)")
    xlabel("X (km)")
    ylabel("Depth (km)")
    ϕ_ps[inds[num]] = res[:ϕ_trace]
end
tight_layout()
wsave(plotsdir("fwi_overthrust", "probed_fwi.png"), gcf())
# Traces

depth = range(0, ee[3], length=n[2])

for t in trinds
    fig, axs = subplots(nrows=2, ncols=4, figsize=(9, 5), sharex=true, sharey=true)
    fig.subplots_adjust(wspace=0, hspace=0)
    for i=1:8
        axs[i].plot(vp_t[:, t], depth, label="True")
        axs[i].plot(vp_std[:, t], depth, label="FWI")
        axs[i].plot(vp_ps[i][:, t], depth, label="PFWI $(2^i)")
        axs[i].plot(vp_dft[i][:, t], depth, label="OTDFT $(2^i)")
        axs[i].set_ylim(depth[end], depth[1])
        axs[i].legend(loc="lower left")
        i%2 == 0 && axs[i].set_xlabel(L"$V_p$")
        i < 3 && axs[i].set_ylabel("Depth (km)")
    end
    title("X=$((t-1)*25/1000) km")
    # tight_layout()
    wsave(plotsdir("fwi_overthrust", "vertical_trace_$(t).png"), gcf())
end

for t in trinds
    fig, axs = subplots(nrows=1, ncols=3, figsize=(12, 9), sharey=true)
    fig.subplots_adjust(wspace=0)
    for (n,i)=enumerate([1, 4, 5])
        axs[n].plot(vp_t[:, t], depth, label="True")
        axs[n].plot(vp_std[:, t], depth, label="FWI")
        axs[n].plot(vp_ps[i][:, t], depth, label="PFWI $(2^i)")
        axs[n].plot(vp_dft[i][:, t], depth, label="OTDFT $(2^i)")
        axs[n].set_ylim(depth[end], depth[1])
        axs[n].legend(loc="lower left")
        n==1 && axs[n].set_ylabel("Depth (km)")
        n==2 && axs[n].set_xlabel(L"$V_p$")
        n==2 && axs[n].set_title("X=$((t-1)*25/1000) km")
    end
    # tight_layout()
    wsave(plotsdir("fwi_overthrust", "vertical_trace_$(t)_select.png"), gcf())
end

# Convergence
cols = ["b", "r", "y", "c", "g", "m", :lime, :indigo]

figure(figsize=(12, 9))
for i=1:8
    plot(ϕ_ps[i]/ϕ_ps[i][1], color=cols[i], linestyle=:dashed, marker="o", label="PFWI $(2^i)")
end
plot(ϕ_std/ϕ_std[1], label="STD")
for i=1:8
    plot(ϕ_dft[i]/ϕ_dft[i][1], color=cols[i], linestyle=:dashed, marker="^", label="OTDFT $(2^i)")
end
legend(loc="upper right", ncol=2)
xlabel("Iteration")
ylabel("Normalized objective")

wsave(plotsdir("fwi_overthrust", "convergence.png"), gcf())

# FWI plot for SEG abstract
figure(figsize=(12, 9))

subplot(3,3,1)
imshow(vp_t; plt_dict...)
title("True")
xlabel("X (km)")
ylabel("Depth (km)")
vlines(x=[t*d[1]/1000 for t=trinds], colors=:k, ymin=0, ymax=ee[3])

subplot(3,3,2)
imshow(vp_0; plt_dict...)
title("Initial")
xlabel("X (km)")
ylabel("Depth (km)")

subplot(3,3,3)
imshow(vp_std; plt_dict...)
title("FWI")
xlabel("X (km)")
ylabel("Depth (km)")

for (j, ps)=zip([4,5,6], [4, 16, 64])
    i = inds[ps]
    subplot(3, 3, j)
    imshow(vp_ps[i]; plt_dict...)
    title("PFWI $(ps)")
    xlabel("X (km)")
    ylabel("Depth (km)")

    subplot(3, 3, j+3)
    imshow(vp_dft[i]; plt_dict...)
    title("OTDFT $(ps)")
    xlabel("X (km)")
    ylabel("Depth (km)")
end
tight_layout()

wsave(plotsdir("fwi_overthrust", "probed_vs_dft.png"), gcf())

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
"A type for matched waves."
mutable struct MatchPlaneWaveMode{T<:AbstractFloat,Dim}
    PlaneWaveModes::Vector{EffectivePlaneWaveMode{T,Dim}}
    discrete_wave::DiscretePlaneWaveMode{T}
    x_match::Vector{T} # waves are matched between discrete_wave.x_match
end

"""
    match_error(m_wave::MatchPlaneWaveMode{T}, shape::Shape

Calculates the difference between the match of MatchPlaneWaveMode.PlaneWaveModes and MatchPlaneWaveMode.discrete_wave. When this is small, we know that the method has likely converged.
"""
function match_error(m_wave::MatchPlaneWaveMode{T}, shape::Shape; apply_norm::Function=norm) where T<:AbstractFloat
    avg_eff = DiscretePlaneWaveMode(m_wave.x_match, m_wave.PlaneWaveModes,shape)
    j0 = findmin(abs.(m_wave.discrete_wave.x .- m_wave.x_match[1]))[2]
    len = length(m_wave.x_match)

    return apply_norm(m_wave.discrete_wave.amplitudes[j0:end,:,:][:] - avg_eff.amplitudes[:])/len
end

function MatchPlaneWaveMode(ω::T, source::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
        tol::T = T(1e-5),
        max_size::Int = 200,
        wave_effs::Vector{EffectivePlaneWaveMode{T,2}} = EffectivePlaneWaveMode{T,2}[],
        x::AbstractVector{T} = [-one(T)],
        L_match::Int = 0,
        kws...
    ) where T<:Number

    θin = transmission_angle(source,material)
    k = real(ω / source.medium.c)

    if isempty(wave_effs)
        wave_effs = WaveModes(ω, source, material;
            extinction_rescale = false,
            tol = T(10)*tol,
            kws...)
    end
    basis_order = wave_effs[1].basis_order

    # use non-dimensional effective waves
    wave_non_effs = map(wave_effs) do w
        EffectivePlaneWaveMode(ω, w.wavenumber / k, w.basis_order, w.direction, w.eigenvectors)
    end

    if first(x) == - one(T)
        # using non-dimensional wave_non_effs and a12k results in non-dimensional mesh X
        a12k = T(2) * k * material.species[1].exclusion_distance * outer_radius(material.species[1])
        L_match, X =  x_mesh_match(wave_non_effs; a12 = a12k, tol = tol, max_size=max_size)
    else
        X = x .* real(k)
        if L_match == 0
            L_match = Int(round(length(X)/2))
        end
    end

    # Normalise the amplitude of the effective waves
    # discrete_waves = [
    #     DiscretePlaneWaveMode(X[L_match:L_match+1], w, material.shape)
    # for w in wave_non_effs]

    # We just scale w.eigenvectors to avoid rounding errors
    # wave_non_effs = map(eachindex(wave_non_effs)) do i
    #     w = wave_non_effs[i]
    #     # amps = w.eigenvectors / norm(discrete_waves[i].eigenvectors[1,:,1])
    #     amps = w.eigenvectors * abs(exp(- im * w.wavenumber * X[L_match]))
    #     EffectivePlaneWaveMode(ω, w.wavenumber, w.basis_order, w.direction, amps)
    # end;

    J = length(collect(X)) - 1
    len = (J + 1)  * (2basis_order + 1)

    (MM_quad,b_mat) = discrete_wave_system(ω, X, source, material; tol=tol, basis_order=basis_order, kws...);
    MM_mat = reshape(MM_quad, (len, len));
    b = reshape(b_mat, (len));

    (LT_mat, ER_mat, b_eff) = match_arrays(ω, wave_non_effs, L_match, X, source, material);

    B = b - ER_mat*b_eff
    As = (ER_mat*LT_mat + MM_mat)\B
    As_mat = reshape(As, (J+1, 2basis_order+1, 1))

    αs = LT_mat*As + b_eff
    # use these αs to correct the magnitude of the amplitudes of the effective waves
    # and re-dimensionalise the effective wavenumbers
    wave_non_effs = map(eachindex(wave_effs)) do i
        w = wave_non_effs[i]
        amps = αs[i] .* w.eigenvectors
        EffectivePlaneWaveMode(ω, k * w.wavenumber, w.basis_order, w.direction, amps)
    end

    mwave = MatchPlaneWaveMode(
        wave_non_effs,
        DiscretePlaneWaveMode(basis_order, real.( collect(X) ./ k), As_mat),
        real.(collect(X[L_match:end]) ./ k)
    )

    println("The relative match error was: ", match_error(mwave, material.shape))

    # return MatchPlaneWaveMode(wave_effs, DiscretePlaneWaveMode(basis_order, collect(X)./k, As_mat), collect(X[L_match:end])./k)
    return mwave
end

"Returns (x,L), where x[L:end] is the mesh used to match with wave_effs."
function x_mesh_match(wave_non_effs::Vector{EffectivePlaneWaveMode{T,Dim}}; kws... ) where {T<:AbstractFloat,Dim}
    # wave_effs[end] establishes how long X should be, while wave_effs[1] establishes how fine the mesh should be.

   # If there is only one wave, then it doesn't make sense to extend the mesh until it decays.
   # Instead we choose, arbitrarily, a quarter of the wavelength.
   # Note, having only one wave is very unusual, but tends to happen in the very low frequency limit.
    x = if length(wave_non_effs) > 1
            x_mesh(wave_non_effs[end], wave_non_effs[1]; kws...)
        else
            x_mesh(wave_non_effs[1]; max_x = (pi/2) / abs(wave_non_effs[1].wavenumber), kws...)
        end
    #NOTE previously used abs(cos(wave_effs[1].θ_eff)*abs(wave_effs[1].wavenumber)) instead of abs(wave_effs[1].wavenumber)

    x_match = x[end]
    x_max = (length(x) < length(wave_non_effs)*T(1.5)) ?
         x[end] + T(1.5)*(x[2] - x[1])*length(wave_non_effs) :
         T(2) * x[end]

    x = 0.0:(x[2] - x[1]):x_max

    L = findmin(abs.(x .- x_match))[2]

    return L, x
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
immutable HyperboloidIntercept{T} <: Model
    β₀::T
    a::T
    α::T
    μ::T
end

function Base.show(io::IO, ::Type{HyperboloidIntercept})
    return print(io, "Intercept Hyperboloid")
end

function Base.show(io::IO, m::HyperboloidIntercept)
    @printf(
        io,
        "Intercept Hyperboloid (β₀ = %.3f, a = %.3f, α = %.3f, μ = %.3f)",
        m.β₀,
        m.a,
        m.α,
        m.μ,
    )
end

function initial_parameters(::Type{HyperboloidIntercept})
    return Float64[0.0, log(1.0), log(1.0), log(1.0)]
end

canonical_parameters(m::HyperboloidIntercept) = (m.a, m.α, m.μ, m.β₀)

function HyperboloidIntercept(parameters::Vector)
    β₀ = parameters[1]
    a = exp(parameters[2])
    α = exp(parameters[3])
    μ = exp(parameters[4])
    return HyperboloidIntercept(β₀, a, α, μ)
end

@inline function predict(
    m::HyperboloidIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, α, μ = m.β₀, m.a, m.α, m.μ
    u1 = x1 * (1 + α * t1)^(-μ)
    u2 = x2 * (1 + α * t2)^(-μ)
    z = β₀ + a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::HyperboloidIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, α, μ = m.β₀, m.a, m.α, m.μ
    u1 = x1 * (1 + α * t1)^(-μ)
    u2 = x2 * (1 + α * t2)^(-μ)
    gr = (
        1.0,
        (u2 - u1) * a,
        a * (x2 * -μ * (1 + α * t2)^(-μ - 1) * t2 - x1 * -μ * (1 + α * t1)^(-μ - 1) * t1) * α,
        a * (-x2 * (1 + α * t2)^(-μ) * log(1 + α * t2) + x1 * (1 + α * t1)^(-μ) * log(1 + α * t1)) * μ,
    )
    return gr
end


text:
using Compat
using Test
using MPI

MPI.Init()

function scatter_array(A, root)
    comm = MPI.COMM_WORLD
    T = eltype(A)
    if MPI.Comm_rank(comm) == root
        global B = copy(A)
    else
        global B = Array{T}(undef,1)
    end
    C = MPI.Scatter(B, 1, root, comm)
end

comm = MPI.COMM_WORLD
root = 0

A = collect(1:MPI.Comm_size(comm))
B = scatter_array(A, root)
@test B[1] == MPI.Comm_rank(comm) + 1
for typ in Base.uniontypes(MPI.MPIDatatype)
    global A = convert(Vector{typ},collect(1:MPI.Comm_size(comm)))
    global B = scatter_array(A, root)
    @test B[1] == convert(typ,MPI.Comm_rank(comm) + 1)
end

MPI.Finalize()


text:
"""
    HasAnatomicalStructures{T}

"""
struct HasAnatomicalStructures{T} end

struct HasVertices end

struct HasFaces end

struct HasCoordinates end

"""
    AnatomicalAtlas

Serves as a metadata to assist other types (arrays, meshes, etc.) in defining
atlases.
"""
struct AnatomicalAtlas{A,I,C,P,H}
    anatomy::A
    identifier::I
    colors::C
    positions::P
    hierarchy::H
end

# This is the GL type we try to immitate
# - GLNormalAttributeMesh
# - GLNormalMesh
# - GLNormalUVMesh
# - GLNormalVertexcolorMesh
# - GLPlainMesh
#


text:
""" Model{G,S,D,L,P}

Structure for representing the petri net model

represented by grounding, states, transition function, transition rates, and predicates
"""
struct Model{G,S,D,L,P}
  g::G  # grounding
  S::S  # states
  Δ::D  # transition function
  Λ::L  # transition rate
  Φ::P  # if state should happen
end

""" Model(s::S, δ::D, λ::L, ϕ::P)

Constructor to initialize a Petri net
"""
Model(s::S, δ::D, λ::L, ϕ::P) where {S,D,L,P} = Model{Any,S,D,L,P}(missing, s, δ, λ, ϕ)

""" Model(s::S, δ::D)

Constructor to initialize a Petri net with just states and transition functions
"""
Model(s::S, δ::D) where {S<:Vector,D<:Vector{Tuple{Operation, Operation}}} = Model(s, δ, [],[])

samelist(x,y) = length(x) == length(y) && all(isequal.(x, y))

function ==(x::Petri.Model,y::Petri.Model)
    return samelist(x.S, y.S) && samelist(x.Δ, y.Δ) && samelist(x.Λ, y.Λ) && samelist(x.Φ, y.Φ)
end

abstract type AbstractProblem end

""" Problem{M<:Model, S, N}

Structure for representing a petri net problem

represented by a petri net model, initial state, and number of steps
"""
struct Problem{M<:Model, S, N} <: AbstractProblem
  m::M
  initial::S
  steps::N
end

struct ParamProblem{M<:Model, S, V, N} <: AbstractProblem
  m::M
  initial::S
  param::V
  steps::N
end


text:
module Connectomes

using LightXML
using DataFrames
using SparseArrays
using SimpleWeightedGraphs
using Graphs
using Makie
using DelimitedFiles
using FileIO
using LinearAlgebra
using Colors
using Serialization
import Graphs: degree

include("graphs.jl")
export Connectome
export graph_filter
export degree
export Connectome2FS
export FS2Connectome
export read_cmtk_parcellation

include("plotting.jl")
export plot_cortex
export plot_cortex!
export plot_mesh
export testplot
export plot_roi
export plot_roi!
export plot_connectome
export plot_parc
export plot_vertex
export plot_vertex!
export plot_edges!

end


text:
@testset "function as actor" begin
    A = @actor () -> @info "Hello World!"
    @test isnothing(A[])
    @test isnothing(A()[])

    A = @actor msg -> @info "Hello " * msg
    @test isnothing(A["World!"])
    @test isnothing(A("World!")[])

    A = @actor msg -> "Hello " * msg
    @test isnothing(A["World!"])
    @test A("World!")[] == "Hello World!"
end

@testset "struct as actor" begin
    Base.@kwdef mutable struct Counter
        n::Int = 0
    end

    (c::Counter)() = c.n += 1

    A = @actor Counter() name="counter"

    @test nameof(A) == "counter"

    Threads.@threads for _ in 1:1_000
        A[]
    end

    @test A.n[] == 1_000

    Threads.@threads for _ in 1:1_000
        A()
    end

    @test A.n[] == 2_000
end

@testset "SysMsg is processed immediately" begin
    Base.@kwdef struct TestSysMsgArrival
        msgs::Vector{Int} = []
    end

    OL.handle(a::TestSysMsgArrival, ::OL.StartMsg) = push!(a.msgs, 0)
    function (a::TestSysMsgArrival)(i::Int)
        sleep(0.1)
        push!(a.msgs, i)
    end

    A = @actor TestSysMsgArrival()

    begin
        A[1]
        A[2]
        A[3]
        A[OL.StartMsg(nothing)]
    end
    @test A.msgs[] == [1,0,2,3]
end

@testset "Error handling" begin
    A = @actor iseven
    @test A(2)[] == true
    @test_throws MethodError A(:x)[]
    # the actor now should have restarted
    @test A(2)[] == true
end


text:
const RANKS = "A23456789TJQKA" # wrapped to create the straights
const SUITS = "cdhs"
const DECK = [string(f, s) for f in collect(RANKS[1:13]), s in collect(SUITS)] |> vec |> sort
const STARTING_HANDS = collect(combinations(DECK, 2))


text:
imu = 6 ± 0.4

u = complex(5 ± 0.3, imu)
v = complex(imu, 8 ± 0.9)

# Addition and subtraction
test_approx_eq(2u - u - 3v + 2v + pi, u - v + pi)

# Multiplication and power
test_approx_eq(v*v*u + u*u*v, u^2*v + v^2*u)

# Division
test_approx_eq(v/v + u/u, complex(2 ± 0))

# Inverse
test_approx_eq(inv(v) + inv(u), (v + u)/(v*u))

# Cosine
for a in (u, v); test_approx_eq(cos(acos(a)), a); end
for a in (u, v); test_approx_eq(cosh(acosh(a)), a); end

# Sine
for a in (u, v); test_approx_eq(sin(asin(a)), a); end
for a in (u, v); test_approx_eq(sinh(asinh(a)), a); end

# Tangent
for a in (v, u); test_approx_eq_eps(tan(a), sin(a)/cos(a), 3e-16); end
for a in (u, v); test_approx_eq(tan(atan(a)), a); end
for a in (u, v); test_approx_eq_eps(tanh(a), sinh(a)/cosh(a), 3e-16); end
for a in (u, v); test_approx_eq(tanh(atanh(a)), a); end

# Exponential
for a in (u, v); test_approx_eq(exp(a), e^a); end
for a in (u, v); test_approx_eq(expm1(a), exp(a) - 1); end
for a in (u, v); test_approx_eq(exp10(a), 10^a); end

# Logarithm
for a in (u, v); test_approx_eq(exp(log(a)), a); end
for a in (u, v); test_approx_eq(log(e, a), log(a)); end
for a in (u, v); test_approx_eq(log(10, a), log10(a)); end
for a in (u, v); test_approx_eq(log1p(a), log(1 + a)); end

# Square root
for a in (u, v); test_approx_eq(sqrt(a), a^0.5); end
for a in (u, v); test_approx_eq(sqrt(a)*sqrt(a), a); end

# Absolute value
for a in (u, v); test_approx_eq(abs(a), sqrt(abs2(a))); end

# TODO: Define a method for lgamma with complex argument.
# # Factorial and gamma
# for a in (u, v); test_approx_eq(gamma(a), factorial(a - 1)); end
# for a in (u, v); test_approx_eq(gamma(a + 1), factorial(a)); end
# for a in (u, v); test_approx_eq(lgamma(a), log(gamma(a))); end


text:
module GitTools

struct GitConfig
    user_local::String
    user_email::String
    user_github::String
    def_branch::String
    GitConfig() = new("","","","")
    GitConfig(github::String ; user::String, email::String, branch::String) = new(user, email, github, branch)
end

# We do not use the StdLib LibGit2
# We actually want to check that Git is installed locally and can be run...

const requireversion = v"2.10"
const gitversion = Ref("<undef>")
function queryversion()
    try
        v = strip(read(`git version`, String))
        if !startswith(v, "git version ")
            @warn "git is installed but its version string does not start with `git version`, assume $v >= $requireversion."
            return v # we assume it is at least $requireversion...
        else
            v = SubString(v, length("git version ")+1)
            m = match(r"^\d+\.\d+", v) # major.minor
            if isnothing(m)
                @warn "git is installed with strange version string, $v >= $requireversion."
                return v # we assume it is at least $requireversion...
            else
                if VersionNumber(m.match) < requireversion
                    @error "git is installed at version $(m.match), we require >= $requireversion."
                    return "" # we known our config will fail
                end
                return v # cleaned version string
            end
        end
        return v
    catch # git not installed
        return ""
    end
end

function isinstalled()
    if gitversion[] == "<undef>"
        gitversion[] = queryversion()
    end
    return gitversion[] != ""
end

function isgitrepo(path::String)
    ispath(path) && ispath(joinpath(path, ".git"))
end

function readconfig(keyname::String)
    try
        String(strip(read(`git config $keyname`, String)))
    catch
         ""
    end
end

function getconfig(path::Union{String,Nothing}=nothing)
    conf = GitConfig()
    isinstalled() || return conf

    dir = nothing
    if !isnothing(path)
        dir = pwd()
        cd(path) # OK to fail if invalid
    end

    try
        conf = GitConfig(readconfig("github.user");
            user   = readconfig("user.name"),
            email  = readconfig("user.email"),
            branch = readconfig("init.defaultBranch"))
    finally
        isnothing(dir) || cd(dir)
    end

    return conf
end

function commit(gitrepo::String, msg::String)
    curr = pwd()
    cd(gitrepo)
    execute(`git commit -a -m "$msg"`)
    cd(curr)
end

# function checkconfig()

#     version = tryversion()
#     if isnothing(version)
#         @error "Git does not seem installed! The command `git version` failed.
#         Install Git or ensure that it is in the PATH, on Windows you may run
#         from the Git Bash shell"
#     end
#     @info "Git installed: $version"

#     user_git = user_github()
#     if isnothing(user_git)
#         @error "No GitHub user defined. Execute the following command
#         to configure the name of your GitHub account.
#                 git config --global github.user your_account_name
#         "
#     end
#     @info "GitHub account: $user_git"

#     user_loc = user_local()
#     if isnothing(user_loc)
#         @error "No local Git user defined. Execute the following command
#         to configure the local user. It is recommanded to user the same
#         user name as your GitHub user account.
#                 git config --global user.name \"your_user_name\"
#         "
#     end
#     @info "Git user: $user_loc"

#     user_mail = user_email()
#     if isnothing(user_mail)
#         @error "No user email defined. Execute the following command
#         to configure the name of your local git contact email.
#                 git config --global user.email name@example.com
#         "
#     end
#     @info "Git email: $user_mail"

#     branch = defaultbranch()
#     if isnothing(branch)
#         @error "No default branch name defined. Execute
#                 git config --global init.defaultBranch main"
#     end

#     return GitConfig(user_loc, user_mail, user_git, branch)

# end

end # module GitTools


text:
module RefCellGallery

using UnstructuredGrids.Core

export SEGMENT
export TRIANGLE
export SQUARE
export LEX_SQUARE
export LEX_HEXAHEDRON
export TETRAHEDRON

const SEGMENT = RefCell(
  ndims = 1,
  faces = [ [[1],[2]] ],
  facetypes = [ [1,1] ],
  reffaces = [ [VERTEX] ],
  coordinates = Float64[-1 1;],
  vtkid = 3,
  vtknodes = [1,2] )

const TRIANGLE = RefCell(
  ndims = 2,
  faces = [ [[1],[2],[3]], [[1,2],[2,3],[3,1]] ],
  facetypes = [ [1,1,1], [1,1,1] ],
  reffaces = [ [VERTEX], [SEGMENT] ],
  coordinates = Float64[ 0 1 0; 0 0 1],
  vtkid = 5,
  vtknodes = [1,2,3])

const TETRAHEDRON = RefCell(
  ndims = 3,
  faces = [
    [[1],[2],[3],[4]],
    [[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]],
    [[1,3,2],[1,2,4],[1,4,3],[2,3,4]] ],
  facetypes = [ [1,1,1,1], [1,1,1,1,1,1], [1,1,1,1] ],
  reffaces = [ [VERTEX], [SEGMENT], [TRIANGLE] ],
  coordinates = Float64[ 0 1 0 0; 0 0 1 0; 0 0 0 1],
  vtkid = 10,
  vtknodes = [1,2,3,4])

const SQUARE = RefCell(
  ndims = 2,
  faces = [ [[1],[2],[3],[4]], [[1,2],[2,3],[3,4],[4,1]] ],
  facetypes = [ [1,1,1,1], [1,1,1,1] ],
  reffaces = [ [VERTEX], [SEGMENT] ],
  coordinates = Float64[ -1 1 1 -1; -1 -1 1 1],
  vtkid = 9,
  vtknodes = [1,2,3,4])

const LEX_SQUARE = RefCell(
  ndims = 2,
  faces = [ [[1],[2],[3],[4]], [[1,2],[3,4],[1,3],[2,4]] ],
  facetypes = [ [1,1,1,1], [1,1,1,1] ],
  reffaces = [ [VERTEX], [SEGMENT] ],
  coordinates = Float64[ -1 1 -1 1; -1 -1 1 1],
  vtkid = 9,
  vtknodes = [1,2,4,3])

const LEX_HEXAHEDRON = RefCell(
  ndims = 3,
  faces = [
    [[1],[2],[3],[4],[5],[6],[7],[8]],
    [[1,2],[3,4],[1,3],[2,4],[5,6],[7,8],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]],
    [[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]]],
  facetypes = [ fill(1,8), fill(1,12), fill(1,6) ],
  reffaces = [ [VERTEX], [SEGMENT], [LEX_SQUARE] ],
  coordinates = Float64[ -1 1 -1 1 -1 1 -1 1; -1 -1 1 1 -1 -1 1 1; -1 -1 -1 -1 1 1 1 1],
  vtkid = 12,
  vtknodes = [1,2,4,3,5,6,8,7])

end # module RefCellGallery


text:
using Revise
using JuMP, EAGO

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 1000,
                                                    "iteration_limit" => 100000,

                                                    "cp_depth" => -1,
                                                    "cp_repetitions" => -1,
                                                    "cp_forward_reverse_limit" => 2,

                                                    "cut_min_iterations" => 3,
                                                    "cut_max_iterations" => 3,
                                                    "objective_cut_on" => true,
                                                    "subgrad_tighten" => true,

                                                    "obbt_depth" => 4,
                                                    "obbt_repetitions" => 4,
                                                    "obbt_aggressive_on" => true,

                                                    "upper_bounding_depth" => 4,

                                                    "fbbt_lp_depth" => 100000,
                                                    "fbbt_lp_repetitions" => 3))
# ----- Variables ----- #
x_Idx = Any[2, 3, 4]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[2], 1.0e-6)
JuMP.set_upper_bound(x[2], 1.0)
JuMP.set_lower_bound(x[3], 1.0e-6)
JuMP.set_upper_bound(x[3], 1.0)
JuMP.set_lower_bound(x[4], 1.0e-6)
JuMP.set_upper_bound(x[4], 1.0)

# ----- Constraints ----- #
@constraint(m, e2, x[2]+x[3]+x[4] == 1.0)

# ----- Objective ----- #
@NLobjective(m, Min, ((15.3261663216011*x[2]+23.2043471859416*x[3]+6.69678129464404*x[4])*log(2.1055*x[2]+3.1878*x[3]+0.92*x[4])+1.04055250396734*x[2]-2.24199441248417*x[3]+3.1618173099828*x[4]+6.4661663216011*x[2]*log(x[2]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4]))+12.2043471859416*x[3]*log(x[3]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4]))+0.696781294644034*x[4]*log(x[4]/(2.1055*x[2]+3.1878*x[3]+0.92*x[4]))+9.86*x[2]*log(x[2]/(1.972*x[2]+2.4*x[3]+1.4*x[4]))+12*x[3]*log(x[3]/(1.972*x[2]+2.4*x[3]+1.4*x[4]))+7*x[4]*log(x[4]/(1.972*x[2]+2.4*x[3]+1.4*x[4]))+(1.972*x[2]+2.4*x[3]+1.4*x[4])*log(1.972*x[2]+2.4*x[3]+1.4*x[4])+1.972*x[2]*log(x[2]/(1.972*x[2]+0.283910843616504*x[3]+3.02002220174195*x[4]))+2.4*x[3]*log(x[3]/(1.45991339466884*x[2]+2.4*x[3]+0.415073537580851*x[4]))+1.4*x[4]*log(x[4]/(0.602183324335333*x[2]+0.115623371371275*x[3]+1.4*x[4]))-17.2981663216011*x[2]*log(x[2])-25.6043471859416*x[3]*log(x[3])-8.09678129464404*x[4]*log(x[4])))

optimize!(m)


text:
for y in (DungBase, Dates), x in names(y, all = true)
    T = getproperty(y, x)
    if T isa DataType
        JSON3.StructType(::Type{T}) = JSON3.Struct()
    elseif T isa UnionAll
        JSON3.StructType(::Type{<:T}) = JSON3.Struct()
    end
end

file = joinpath(Databula.sourcefolder, "video.json")
videos = open(file, "r") do io
    JSON3.read(io, Vector{Union{WholeVideo, FragmentedVideo, DisjointVideo}})
end
serialize(joinpath(Databula.sourcefolder, "video"), videos)

#=file = joinpath(Databula.sourcefolder, "calibrations.json")
calibrations = open(file, "r") do io
    JSON3.read(io, Vector{Calibration})
end
serialize(joinpath(Databula.sourcefolder, "calibration"), calibrations)=#

rm(joinpath(Databula.sourcefolder, "calibrations.json"))


text:
struct MetaWeights{InnerMetaGraph,U<:Real} <: AbstractMatrix{U}
    meta_graph::InnerMetaGraph
end

show(io::IO, x::MetaWeights) = print(io, "metaweights")
show(io::IO, z::MIME"text/plain", x::MetaWeights) = show(io, x)

MetaWeights(g::MetaGraph) = MetaWeights{typeof(g),weighttype(g)}(g)

function getindex(w::MetaWeights, u::Int, v::Int)
    g = w.meta_graph
    metaindex = g.metaindex
    if has_edge(g, u, v)
        g.weightfunction(g[arrange(g, metaindex[u], metaindex[v], u, v)...])
    else
        g.defaultweight
    end
end

function size(d::MetaWeights)
    vertices = nv(d.meta_graph)
    (vertices, vertices)
end

weights(g::MetaGraph) = MetaWeights(g)

"""
    weighttype(g)

Return the weight type for metagraph `g`.
"""
weighttype(g::MetaGraph{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,Weight}) where {Weight} =
    Weight

"""
    weightfunction(g)

Return the weight function for metagraph `g`.
"""
weightfunction(g::MetaGraph) = g.weightfunction

"""
    defaultweight(g)

Return the default weight for metagraph `g`.
"""
defaultweight(g::MetaGraph) = g.defaultweight


text:
n = 20  # Number of observations
p = 5   # Number of parameters

info("Testing ", MOD.pca_components_eig!)
for T in FloatingPointTypes
    X = rand(T, n, p)
    μ = mean(X,1)
    H = (X .- μ)
    Σ = H'H
    for α in (zero(T), convert(T,0.5), one(T)), ϵ in (zero(T), one(T))
        V1, D1 = MOD.pca_components_eig!(copy(H), α, ϵ)
        V2, D2 = MOD.components_eig!(MOD.perturb!(MOD.regularize!(copy(Σ), α, trace(Σ)/p), ϵ))
        @test_approx_eq V1 V2
        @test_approx_eq D1 D2
    end
end

info("Testing ", MOD.pca_components_svd!)
for T in FloatingPointTypes
    X = rand(T, n, p)
    μ = mean(X,1)
    H = (X .- μ)
    V1, D1 = MOD.pca_components_svd!(copy(H))
    V2, D2 = MOD.components_svd!(copy(H))
    @test_approx_eq V1 V2
    @test_approx_eq D1 D2
end

info("Testing ", MOD.pca)
for T in FloatingPointTypes
    X = rand(T, n, p)
    μ = mean(X,1)
    H = (X .- μ)
    Σ = H'H/(n-1)
    for α in (zero(T), convert(T,0.5), one(T)), ϵ in (zero(T), one(T))
        V1, D1 = MOD.pca(copy(X), α, ϵ)
        D2, V2 = eig(MOD.perturb!(MOD.regularize!(copy(Σ), α, trace(Σ)/p), ϵ))
        @test_approx_eq abs(V1) abs(V2[:,end:-1:1])
        @test_approx_eq abs(D1) abs(D2[end:-1:1])
        @test_approx_eq abs(transform_pca(V1,X)) abs(X*(V2[:,end:-1:1]))
    end
end

info("Testing ", MOD.kpca)
for T in FloatingPointTypes
    X = rand(T, n, p)
    Xc = X .- mean(X,1)
    Σ = Xc'Xc/(n-1)
    κ = PolynomialKernel(one(T),zero(T),one(T))  # Equivalent to dot product
    for ϵ in (zero(T), one(T))
        V1, D1 = MOD.kpca(copy(X), κ, zero(T), ϵ)
        V2, D2 = MOD.pca(copy(X), zero(T), ϵ)
        @test_approx_eq (V1'Xc*Xc'V1/(n-1))[1:p,1:p] (V2'Σ*V2)
        @test_approx_eq abs(D1)[1:p] abs(D2)
    end
end


text:
using AugmentedGaussianProcesses
using Random: seed!
using MLKernels, Statistics, LinearAlgebra
using ForwardDiff
const AGP = AugmentedGaussianProcesses
seed!(42)
dims= 2; nA = 10; nB = 5
A = rand(nA,dims)
B = rand(nB,dims)

function adapt_AD(kernel::AGP.Kernel{T},X::AbstractMatrix{T2}) where {T,T2}
    T.(X)
end

function create_kernel(k,θ,name;ν=0)
    if name == "RBFKernel"
        return k(θ)
    elseif name == "MaternKernel"
        return k(θ,ν)
    end
end
#Compare kernel results with MLKernels package
@testset "Kernels" begin

    @testset "Kernel computation" begin
        θ = 0.5; ν = 2.0
        @testset "Iso Kernels" begin
            kernels_GP = [AGP.RBFKernel(θ),AGP.MaternKernel(θ,ν)]
            kernels_ML = [MLKernels.GaussianKernel(0.5/θ^2),MLKernels.MaternKernel(ν,θ)]
            test_name = ["RBFKernel","MaternKernel"]
            for (k_GP,k_ML,t_name) in zip(kernels_GP,kernels_ML,test_name)
                @testset "$t_name" begin
                    @test sum(abs.(MLKernels.kernelmatrix(k_ML,A)- AGP.kernelmatrix(A,k_GP))) ≈ 0 atol=1e-6
                    @test sum(abs.(MLKernels.kernelmatrix(k_ML,A,B)- AGP.kernelmatrix(A,B,k_GP))) ≈ 0 atol=1e-6
                    @test sum(abs.(diag(AGP.kernelmatrix(A,k_GP))-AGP.kerneldiagmatrix(A,k_GP))) ≈ 0 atol=1e-6
                end
            end
        end
        @testset "ARD Kernels" begin
            kernels_Iso = [AGP.RBFKernel(θ),AGP.MaternKernel(θ,ν)]
            kernels_ARD = [AGP.RBFKernel(fill(θ,dims)),AGP.MaternKernel(fill(θ,dims),ν)]
            test_name = ["RBFKernel","MaternKernel"]
            for (k_iso,k_ard,t_name) in zip(kernels_Iso,kernels_ARD,test_name)
                @testset "$t_name" begin
                    @test sum(abs.(AGP.kernelmatrix(A,k_iso)- AGP.kernelmatrix(A,k_ard))) ≈ 0 atol=1e-6
                    @test sum(abs.(AGP.kernelmatrix(A,B,k_iso)- AGP.kernelmatrix(A,B,k_ard))) ≈ 0 atol=1e-6
                    @test sum(abs.(diag(AGP.kernelmatrix(A,k_ard))-AGP.kerneldiagmatrix(A,k_ard))) ≈ 0 atol=1e-6
                end
            end
        end
    end

    @testset "Kernel derivatives" begin
        θ = 0.5; ν = 2.0
        @testset "Iso Kernels" begin
            kernels = [AGP.RBFKernel(θ)]#,AGP.MaternKernel(θ,ν)]
            kernels_AD = [AGP.RBFKernel]#,AGP.MaternKernel]
            test_name = ["RBFKernel"]#,"MaternKernel"]
            for (k,k_AD,t_name) in zip(kernels,kernels_AD,test_name)
                @testset "$t_name" begin
                    dA_ad = reshape(ForwardDiff.jacobian(x->begin;k_ = create_kernel(k_AD,x[1],t_name,ν=ν); AGP.kernelmatrix(adapt_AD(k_,A),k_); end,[θ]),nA,nA)
                    dAB_ad = reshape(ForwardDiff.jacobian(x->begin;k_ = create_kernel(k_AD,x[1],t_name,ν=ν); AGP.kernelmatrix(adapt_AD(k_,A),adapt_AD(k_,B),k_); end,[θ]),nA,nB)
                    @test sum(abs.(AGP.kernelderivativematrix(A,k).- dA_ad)) ≈ 0 atol=1e-6
                    @test sum(abs.(AGP.kernelderivativematrix(A,B,k)- dAB_ad)) ≈ 0 atol=1e-6
                    @test sum(abs.(AGP.kernelderivativediagmatrix(A,k)-diag(dA_ad))) ≈ 0 atol=1e-6
                end
            end
        end
        @testset "ARD Kernels" begin
            θ = rand(dims)
            kernels = [AGP.RBFKernel(θ)]#,AGP.MaternKernel(θ,ν)]
            kernels_AD = [AGP.RBFKernel]#,AGP.MaternKernel]
            test_name = ["RBFKernel"]#,"MaternKernel"]
            for (k,k_AD,t_name) in zip(kernels,kernels_AD,test_name)
                @testset "$t_name" begin
                    dA_ad = ForwardDiff.jacobian(x->begin;k_ = create_kernel(k_AD,x,t_name,ν=ν); AGP.kernelmatrix(adapt_AD(k_,A),k_); end,θ)
                    dAB_ad = ForwardDiff.jacobian(x->begin;k_ = create_kernel(k_AD,x,t_name,ν=ν); AGP.kernelmatrix(adapt_AD(k_,A),adapt_AD(k_,B),k_); end,θ)
                    @test sum(abs.(hcat(vec.(AGP.kernelderivativematrix(A,k))...).- dA_ad)) ≈ 0 atol=1e-6
                    @test sum(abs.(hcat(vec.(AGP.kernelderivativematrix(A,B,k))...)- dAB_ad)) ≈ 0 atol=1e-6
                    @test sum(abs.(sum(AGP.kernelderivativediagmatrix(A,k).-diag.(reshape.(eachcol(dA_ad),nA,nA))))) ≈ 0 atol=1e-6
                end
            end
        end

    end
end


text:
using Logging

example_scenario = joinpath(@__DIR__, "example_pool_qs.jl")
if !isempty(ARGS)
    length(ARGS) > 1 && error("Please only supply one scenario.")
    example_scenario = joinpath(@__DIR__, "example_$(ARGS[1]).jl")
    isfile(example_scenario) || error("Cannot find scenario file '$(example_scenario)'.")
end
@info "Running with '$example_scenario'"

using SVDD, OneClassActiveLearning, JSON, JuMP, Ipopt, Random, JLD
Random.seed!(0)

function run_experiment(experiment::Dict)
    Random.seed!(0)
    @info "host: $(gethostname()) worker: $(getpid()) exp_hash: $(experiment[:hash])"
    if isfile(experiment[:output_file])
        @warn "Aborting experiment because the output file already exists. Filename: $(experiment[:output_file])"
        return nothing
    end

    res = Result(experiment)
    errorfile = joinpath(experiment[:log_dir],"$(gethostname())_$(getpid()).error")
    try
        time_exp = @elapsed res = OneClassActiveLearning.active_learn(experiment)
        res.al_summary[:runtime] = Dict(:time_exp => time_exp)
    catch e
        res.status[:exit_code] = Symbol(typeof(e))
        @warn "Experiment $(experiment[:hash]) finished with unkown error."
        @warn e
    finally
        if res.status[:exit_code] != :success
            @info "Writing error hash to $errorfile."
            open(errorfile, "a") do f
                print(f, "$(experiment[:hash])\n")
            end
        end
        @info "Writing result to $(experiment[:output_file])."
        OneClassActiveLearning.write_result_to_file(experiment[:output_file], res)
    end
end

@info "Generating experiment settings."
include(example_scenario)
@info "Running experiment."
run_experiment(experiment)


text:
# Autogenerated wrapper script for Openresty_jll for x86_64-linux-musl
export luajit_dir, lualib_dir, nginx_dir, openresty, opm, resty

JLLWrappers.@generate_wrapper_header("Openresty")
JLLWrappers.@declare_file_product(luajit_dir)
JLLWrappers.@declare_file_product(lualib_dir)
JLLWrappers.@declare_file_product(nginx_dir)
JLLWrappers.@declare_executable_product(openresty)
JLLWrappers.@declare_executable_product(opm)
JLLWrappers.@declare_executable_product(resty)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        luajit_dir,
        "luajit",
    )

    JLLWrappers.@init_file_product(
        lualib_dir,
        "lualib",
    )

    JLLWrappers.@init_file_product(
        nginx_dir,
        "nginx",
    )

    JLLWrappers.@init_executable_product(
        openresty,
        "bin/openresty",
    )

    JLLWrappers.@init_executable_product(
        opm,
        "bin/opm",
    )

    JLLWrappers.@init_executable_product(
        resty,
        "bin/resty",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "pd" begin
    @testset "scalar" begin
        @test pd(PDGains(1, 2), 3, 4) == -11
    end

    @testset "vector" begin
        gains = PDGains(1, 2)
        e = [3; 4; 5]
        ė = [6; 7; 8]
        @test pd(gains, e, ė) == ((e, ė) -> pd(gains, e, ė)).(e, ė)
    end

    @testset "specifying desireds" begin
        gains = PDGains(5, 6)
        x = SVector(1, 2)
        ẋ = SVector(5, 6)
        @test pd(gains, x, ẋ) == pd(gains, x, zero(x), ẋ, zero(ẋ))
    end

    @testset "x-axis rotation" begin
        Random.seed!(56)
        gains = PDGains(2, 3)
        e = RotX(rand())
        ė = rand() * SVector(1, 0, 0)
        @test pd(gains, e, ė)[1] ≈ pd(gains, e.theta, ė[1])
    end

    @testset "orientation control" begin
        Random.seed!(57)
        mechanism = rand_floating_tree_mechanism(Float64) # single floating body
        joint = first(tree_joints(mechanism))
        body = successor(joint, mechanism)
        base = root_body(mechanism)

        state = MechanismState(mechanism)
        rand!(state)

        Rdes = rand(RotMatrix{3})
        ωdes = zero(SVector{3})
        gains = PDGains(100, 20)
        v̇des = similar(velocity(state))

        control_dynamics_result = DynamicsResult(mechanism)
        function control!(torques::SegmentedVector, t, state::MechanismState)
            H = transform_to_root(state, body)
            T = transform(twist_wrt_world(state, body), inv(H))
            R = rotation(H)
            ω = T.angular
            ωddes = pd(gains, R, Rdes, ω, ωdes)
            v̇desjoint = v̇des[joint]
            v̇desjoint .= SVector([ωddes; zero(ωddes)])
            wrenches = control_dynamics_result.jointwrenches
            accelerations = control_dynamics_result.accelerations
            inverse_dynamics!(torques, wrenches, accelerations, state, v̇des)
        end

        final_time = 3.
        simulate(state, final_time, control!; Δt = 1e-3)

        H = transform_to_root(state, body)
        T = transform(twist_wrt_world(state, body), inv(H))
        R = rotation(H)
        ω = T.angular

        @test isapprox(R * Rdes', one(typeof(Rdes)); atol = 1e-8)
        @test isapprox(ω, zero(ω); atol = 1e-8)
    end

    @testset "pose control" begin
        Random.seed!(58)
        mechanism = rand_floating_tree_mechanism(Float64) # single floating body
        joint = first(tree_joints(mechanism))
        body = successor(joint, mechanism)
        base = root_body(mechanism)

        baseframe = default_frame(base)
        desiredframe = CartesianFrame3D("desired")
        actualframe = frame_after(joint)

        state = MechanismState(mechanism)
        rand!(state)

        xdes = rand(Transform3D, desiredframe, baseframe)
        vdes = zero(Twist{Float64}, desiredframe, baseframe, actualframe)
        v̇des = similar(velocity(state))
        gains = SE3PDGains(baseframe, PDGains(100 * one(SMatrix{3, 3}), 20), PDGains(100., 20.)) # world-fixed gains

        control_dynamics_result = DynamicsResult(mechanism)
        function control!(torques::SegmentedVector, t, state::MechanismState)
            x = transform_to_root(state, body)
            invx = inv(x)
            v = transform(twist_wrt_world(state, body), invx)
            v̇desjoint = v̇des[joint]
            v̇desjoint .= SVector(pd(transform(gains, invx), x, xdes, v, vdes))
            wrenches = control_dynamics_result.jointwrenches
            accelerations = control_dynamics_result.accelerations
            inverse_dynamics!(torques, wrenches, accelerations, state, v̇des)
        end

        final_time = 3.
        simulate(state, final_time, control!; Δt = 1e-3)

        x = transform_to_root(state, body)
        v = transform(twist_wrt_world(state, body), inv(x))
        @test isapprox(x, xdes * one(Transform3D, actualframe, desiredframe), atol = 1e-6)
        @test isapprox(v, vdes + zero(Twist{Float64}, actualframe, desiredframe, actualframe), atol = 1e-6)
    end

    @testset "linearized SE(3) control" begin
        Random.seed!(59)
        for i = 1 : 100
            randpsd3() = (x = rand(SMatrix{3, 3}); x' * x)
            baseframe = CartesianFrame3D("base")
            bodyframe = CartesianFrame3D("body")
            gains = SE3PDGains(bodyframe, PDGains(randpsd3(), randpsd3()), PDGains(randpsd3(), randpsd3()))

            xdes = rand(Transform3D{Float64}, bodyframe, baseframe)
            Tdes = rand(Twist{Float64}, bodyframe, baseframe, bodyframe)

            ϵ = 1e-2
            x = xdes * Transform3D(bodyframe, bodyframe, AngleAxis(ϵ, randn(), randn(), randn()), ϵ * randn(SVector{3}))
            T = rand(Twist{Float64}, bodyframe, baseframe, bodyframe)

            accel_nonlin = pd(gains, x, xdes, T, Tdes)
            accel_lin = pd(gains, x, xdes, T, Tdes, SE3PDMethod{:Linearized}())
            @test isapprox(accel_nonlin, accel_lin; atol = 1e-6)

            T = Tdes
            accel_nonlin = pd(gains, x, xdes, T, Tdes)
            accel_lin = pd(gains, x, xdes, T, Tdes, SE3PDMethod{:Linearized}())
            @test isapprox(accel_nonlin, accel_lin; atol = 1e-6)
        end
    end
end


text:
println("Testing Slice..")

test_mesh = HomogenousMesh(Point{3,Float64}[Point{3,Float64}(0.0,0.0,10.0),
 Point{3,Float64}(0.0,10.0,10.0),
 Point{3,Float64}(0.0,0.0,0.0),
 Point{3,Float64}(0.0,10.0,0.0),
 Point{3,Float64}(10.0,0.0,10.0),
 Point{3,Float64}(10.0,10.0,10.0),
 Point{3,Float64}(10.0,0.0,0.0),
 Point{3,Float64}(10.0,10.0,0.0),
],Face{3,Int,0}[
 Face{3,Int,0}(3,7,5)
 Face{3,Int,0}(1,3,5)
 Face{3,Int,0}(1,2,3)
 Face{3,Int,0}(3,2,4)
 Face{3,Int,0}(1,5,6)
 Face{3,Int,0}(2,1,6)
 Face{3,Int,0}(4,8,3)
 Face{3,Int,0}(3,8,7)
 Face{3,Int,0}(7,8,6)
 Face{3,Int,0}(5,7,6)
 Face{3,Int,0}(2,6,4)
 Face{3,Int,0}(4,6,8)]
)
s = slice(test_mesh, [1.0, 2.0])

let
    @test length(s) == 2
    @test length(s[1]) == length(s[2]) == 1
    exp1 = [([0.0,1.0],[0.0,0.0]),([0.0,0.0],[1.0,0.0]),
                      ([1.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,1.0]),
                      ([10.0,1.0],[10.0,10.0]),([10.0,10.0],[1.0,10.0]),
                      ([1.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,1.0])]
    exp2 = [([0.0,2.0],[0.0,0.0]),([0.0,0.0],[2.0,0.0]),
                      ([2.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,2.0]),
                      ([10.0,2.0],[10.0,10.0]),([10.0,10.0],[2.0,10.0]),
                      ([2.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,2.0])]
    @test length(s[1][1]) == length(exp1)
    @test length(s[2][1]) == length(exp2)
end


text:
struct Dial
    str::String
    Dial(str::AbstractString) = new(chomp(str))
    Dial(str::Vector{String}) = new(join(str, "\n"))
end

abstract type FontSize end
struct Small <: FontSize end
struct Large <: FontSize end

@option "dials" struct Dials
    one::Vector{String}
    two::Vector{String}
    three::Vector{String}
    four::Vector{String}
    five::Vector{String}
    six::Vector{String}
    seven::Vector{String}
    eight::Vector{String}
    nine::Vector{String}
    zero::Vector{String}
    colon::Vector{String}
end

Base.split(d::Dial, dlm) = split(d.str, dlm)

function Base.hcat(dl::Dial, dr::Dial)
    buf = IOBuffer()
    n = min(length(split(dl, "\n")), length(split(dr, "\n")))
    for (i, (l1, l2)) in enumerate(zip(split(dl, "\n"), split(dr, "\n")))
        pfun = i != n ? println : print
        pfun(buf, l1 * " " * l2)
    end
    buf |> take! |> String |> Dial
end

Base.hcat(ds::Dial...) = reduce(hcat, ds)


text:
using GLVisualize, GLWindow, Colors, FixedPointNumbers

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
Direct RGBA rendering using an absorption model.
"""

function volume_data(N)
    vol     = Float32[sin(x/15.0)+sin(y/15.0)+sin(z/15.0) for x=1:N, y=1:N, z=1:N]
    max     = maximum(vol)
    min     = minimum(vol)
    vol     = (vol .- min) ./ (max .- min)
    vol     = map(x->RGBA{Float32}(x, 0, 1-x, 3.0), vol)
end

vol = visualize(volume_data(128), :absorption)
_view(vol, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using BlockArrays, Test

@testset "block dequeue" begin
    @testset "blockappend!(::BlockVector, _)" begin
        @testset for compatible in [false, true],
            srctype in [:BlockVector, :PseudoBlockVector, :PseudoBlockVector2, :Vector]

            dest = mortar([[1, 2, 3], [4, 5]])

            # Create `src` array:
            if compatible
                T = Int
            else
                T = Float64
            end
            if srctype === :BlockVector
                src = mortar([T[6, 7], T[8, 9]])
            elseif srctype === :PseudoBlockVector
                src = PseudoBlockVector(T[6:9;], [4])
            elseif srctype === :PseudoBlockVector2
                src = PseudoBlockVector(T[6:9;], [2, 2])
            elseif srctype === :Vector
                src = T[6:9;]
            else
                error("Unknown srctype = $srctype")
            end

            @test blockappend!(dest, src) === dest
            @test dest == 1:9

            @test dest[Block(1)] == [1, 2, 3]
            @test dest[Block(2)] == [4, 5]
            if blocklength(src) == 2
                @test dest[Block(3)] == [6, 7]
                @test dest[Block(4)] == [8, 9]
            elseif blocklength(src) == 1
                @test dest[Block(3)] == 6:9
            else
                error("Unexpected: blocklength(src) = ", blocklength(src))
            end

            src[1] = 666
            if compatible && srctype !== :PseudoBlockVector2
                @test dest[6] == 666
            else
                @test dest[6] == 6
            end
        end

        @testset "empty blocks" begin
            dest = mortar([[1, 2, 3], [4, 5]])
            @test blockappend!(dest, mortar([Int[]])) === dest == 1:5
            @test blocklength(dest) == 3
            @test blockappend!(dest, mortar([Int[], Int[]])) === dest == 1:5
            @test blocklength(dest) == 5
        end
    end

    @testset "blockpush!(::BlockVector, _)" begin
        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, [6]) === dest == 1:6
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, [6.0]) === dest == 1:6
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, Int[]) === dest == 1:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, (6, 7.0)) === dest == 1:7
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, (x for x in 6:7 if iseven(x))) === dest == 1:6
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpush!(dest, [6], Int[], 7:8) === dest == 1:8
        @test blocklength(dest) == 5
    end

    @testset "blockpushfirst!(::BlockVector, _)" begin
        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, [0]) === dest == 0:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, [0.0]) === dest == 0:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, Int[]) === dest == 1:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, (-1, 0.0)) === dest == -1:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, (x for x in -1:0 if iseven(x))) === dest == 0:5
        @test blocklength(dest) == 3

        dest = mortar([[1, 2, 3], [4, 5]])
        @test blockpushfirst!(dest, [-2], Int[], -1:0) === dest == -2:5
        @test blocklength(dest) == 5
    end

    @testset "blockpop!(::BlockVector, _)" begin
        A = mortar([[1, 2, 3], [4, 5]])
        @test blockpop!(A) == 4:5
        @test A == 1:3
        @test A[Block(1)] == 1:3
    end

    @testset "blockpopfirst!(::BlockVector, _)" begin
        A = mortar([[1, 2, 3], [4, 5]])
        @test blockpopfirst!(A) == 1:3
        @test A == 4:5
        @test A[Block(1)] == 4:5
    end

    @testset "append!(::BlockVector, _)" begin
        @testset "$label" for (label, itr) in [
            "UnitRange" => 6:9,
            "BlockVector" => mortar([[6, 7], [8, 9]]),
            "with length" => (x + 0 for x in 6:9),
            "no length" => (x for x in 6:9 if x > 0),
        ]
            dest = mortar([[1, 2, 3], [4, 5]])
            @test append!(dest, itr) === dest
            @test dest == 1:9
            @test dest[Block(1)] == [1, 2, 3]
            @test dest[Block(2)] == 4:9
        end
    end

    @testset "push!" begin
        A = mortar([[1, 2, 3], [4, 5]])
        push!(A, 6)
        push!(A, 7, 8, 9)
        @test A == 1:9
    end

    @testset "pushfirst!" begin
        A = mortar([[1, 2, 3], [4, 5]])
        pushfirst!(A, 0)
        pushfirst!(A, -3, -2, -1)
        @test A == -3:5
    end

    @testset "pop!" begin
        A = mortar([[1, 2, 3], [4, 5]])
        B = []
        while !isempty(A)
            push!(B, pop!(A))
        end
        @test A == []
        @test B == 5:-1:1
    end

    @testset "popfirst!" begin
        A = mortar([[1, 2, 3], [4, 5]])
        B = []
        while !isempty(A)
            push!(B, popfirst!(A))
        end
        @test A == []
        @test B == 1:5
    end
end


text:
using Test
using LinearAlgebra
using SparseArrays
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allapprox(v1, v2)
    return all(v1 .≈ v2)
end

function allapprox(v1, v2, tol)
    flag = length(v1) == length(v2)
    return flag && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:length(v1)])
end

function allequal(v1, v2)
    return all(v1 .== v2)
end

matrix = HDGElasticity.SystemMatrix()
@test length(matrix.rows) == 0
@test length(matrix.vals) == 0
@test length(matrix.cols) == 0

oprows = rand(1:20, 5)
opcols = rand(1:20, 5)
opvals = rand(5)
HDGElasticity.assemble!(matrix, oprows, opcols, opvals)
@test allapprox(matrix.rows, oprows)
@test allapprox(matrix.cols, opcols)
@test allapprox(matrix.vals, opvals)

rhs = HDGElasticity.SystemRHS()
@test length(rhs.rows) == 0
@test length(rhs.vals) == 0

rows = [1, 2, 5]
vals = rand(3)
HDGElasticity.assemble!(rhs, rows, vals)
@test allequal(rhs.rows, rows)
@test allapprox(rhs.vals, vals)

@test HDGElasticity.element_dof_start(5, 4) == 17
@test HDGElasticity.element_dof_start(6, 6) == 31

@test HDGElasticity.element_dof_stop(4, 8) == 32
@test HDGElasticity.element_dof_stop(6, 8) == 48

@test allequal(HDGElasticity.element_dofs(3, 5), 11:15)

rows = [1, 2, 3]
cols = [4, 5, 6]
r, c = HDGElasticity.operator_dofs(rows, cols)
testrows = [1, 2, 3, 1, 2, 3, 1, 2, 3]
testcols = [4, 4, 4, 5, 5, 5, 6, 6, 6]
@test allequal(r, testrows)
@test allequal(c, testcols)

matrix = HDGElasticity.SystemMatrix()
vals = rand(36)
HDGElasticity.assemble!(matrix, vals, 2, 3, 6)
rowtest = repeat(7:12, 6)
coltest = repeat(13:18, inner = 6)
@test allequal(rowtest, matrix.rows)
@test allequal(coltest, matrix.cols)
@test allapprox(vals, matrix.vals)

matrix = HDGElasticity.SystemMatrix()
vals = rand(108)
HDGElasticity.assemble!(matrix, vals, 3, [2, 3, 5], 6)
rowtest = repeat(13:18, 18)
coltest = vcat((
    repeat(7:12, inner = 6),
    repeat(13:18, inner = 6),
    repeat(25:30, inner = 6),
)...)
@test allequal(matrix.rows, rowtest)
@test allequal(matrix.cols, coltest)
@test allapprox(matrix.vals, vals)

matrix = HDGElasticity.SystemMatrix()
vals = rand(256)
HDGElasticity.assemble!(matrix, vals, [1, 2, 3, 4], [1, 2, 3, 4], 4)
rowtest = repeat(1:16, 16)
coltest = repeat(1:16, inner = 16)
@test allequal(matrix.rows, rowtest)
@test allequal(matrix.cols, coltest)
@test allapprox(matrix.vals, vals)

dofspernode = 2
basis = TensorProductBasis(1, 1)
NF = HDGElasticity.number_of_basis_functions(basis)
dofsperelement = dofspernode * NF
quad = tensor_product_quadrature(1, 2)
HH = HDGElasticity.mass_matrix(basis, quad, 1.0, 2)
rhsvals = repeat(quad.points, inner = (2, 1))
rhs = HDGElasticity.linear_form(rhsvals, basis, quad)
sysmatrix = HDGElasticity.SystemMatrix()
HDGElasticity.assemble_displacement_face!(sysmatrix, basis, quad, 1.0,
    1, dofsperelement)
K = dropzeros!(sparse(
    sysmatrix.rows,
    sysmatrix.cols,
    sysmatrix.vals,
    dofsperelement,
    dofsperelement,
))

sol = K\rhs
sol = reshape(sol,2,:)
testsol = [-1. 1.
           -1. 1.]
@test allapprox(sol,testsol)

vbasis = TensorProductBasis(2,1)
vquad = tensor_product_quadrature(2,2)
sbasis = TensorProductBasis(1,1)
NHF = HDGElasticity.number_of_basis_functions(sbasis)
dofsperelement = NHF*2
squad = tensor_product_quadrature(1,2)
facequads = repeat([squad],4)
facemaps = HDGElasticity.reference_cell_facemaps(2)
normals = HDGElasticity.reference_normals()
lambda = 1.
mu = 2.
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
stabilization = 0.1
cellmap = HDGElasticity.CellMap([1.,1.],[2.,3.])
facescale = HDGElasticity.face_determinant_jacobian(cellmap)

solvercomponents = HDGElasticity.LocalSolverComponents(vbasis,vquad,sbasis,
    facequads,facemaps,normals,D1,stabilization,cellmap)

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

HDGElasticity.assemble_mixed_face!(system_matrix,vbasis,sbasis,facequads[1],
    facemaps[1],normals[1],[0.,1.],[1.,0.],D1,stabilization,facescale[1],
    solvercomponents.iLLxLH,1,1:4,dofsperelement)

HDGElasticity.assemble_mixed_face!(system_matrix,vbasis,sbasis,facequads[4],
    facemaps[4],normals[4],[1.,0.],[0.,1.],D1,stabilization,facescale[4],
    solvercomponents.iLLxLH,4,1:4,dofsperelement)

HDGElasticity.assemble_traction_face!(system_matrix,sbasis,facequads[2],
    facescale[2],stabilization,solvercomponents.fLH[2]',
    solvercomponents.iLLxLH,2,[1,2,3,4],dofsperelement)

HDGElasticity.assemble_traction_face!(system_matrix,sbasis,facequads[3],
    facescale[3],stabilization,solvercomponents.fLH[3]',
    solvercomponents.iLLxLH,3,[1,2,3,4],dofsperelement)

rhsvals = HDGElasticity.linear_form(-[1.,0.],sbasis,squad)
HDGElasticity.assemble!(system_rhs,rhsvals,2,dofsperelement)

K = HDGElasticity.sparse(system_matrix,16)
rhs = HDGElasticity.rhs(system_rhs,16)

sol = K\rhs
L = solvercomponents.iLLxLH*sol
sigma = -D1*reshape(L[1:12],3,:)
displacement = reshape(L[13:20],2,:)

e11 = ((lambda+2mu) - lambda^2/(lambda+2mu))^-1
e22 = -lambda/(lambda+2mu)*e11
u1 = e11*1.0
u2 = e22*2
s11 = 1.0

testsigma = zeros(3,4)
testsigma[1,:] .= 1.0
testu = zeros(2,4)
testu[:,2] .= [0.,u2]
testu[:,3] .= [u1,0.]
testu[:,4] .= [u1,u2]

@test allapprox(testsigma,sigma,1e-10)
@test allapprox(testu,displacement,1e-10)

# ###############################################################################
# # Test a coherent interface
function distance_function(coords,xc)
    return coords[1,:] .- xc
end

polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
NF = HDGElasticity.number_of_basis_functions(levelset.basis)
mesh = UniformMesh([0.,0.],[2.,1.],[1,1])
coords = HDGElasticity.nodal_coordinates(mesh,levelset.basis)
interface_location = 0.5
levelsetcoeffs = reshape(distance_function(coords,interface_location),NF,1)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
NHF = HDGElasticity.number_of_basis_functions(ufs.sbasis)
dofsperelement = 2*NHF

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

lambda,mu = 1.,2.
stabilization = 1e-3
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D1,stabilization)

# Cell 1 phase 1
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,1,[0.,1.],[1.,0.],1,1:4,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,1,1,2,
    2,1:4,dofsperelement)
rhs = -HDGElasticity.linear_form(dgmesh.facescale[2]*[0.1,0.0],ufs.sbasis,ufs.fquads[1,1][2])
HDGElasticity.assemble!(system_rhs,rhs,2,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,3,3,1:4,dofsperelement)

# # Cell 1 phase 2
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,1,[0.,1.],[1.,0.],
    5,5:8,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    2,1,3,6,5:8,dofsperelement)
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,4,
    [1.,0.],[0.,1.],7,5:8,dofsperelement)

# # Interface condition
HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    1,4,1:4,8,5:8,dofsperelement)

K = dropzeros!(HDGElasticity.sparse(system_matrix,32))
R = HDGElasticity.rhs(system_rhs,32)

sol = K\R
H = reshape(sol,4,:)

dofspercell = 4*dofsperelement
H1 = sol[1:dofspercell]
L1 = cellsolvers[1,1].iLLxLH*H1
numstressdofs = NF*3
S1 = -D1*reshape(L1[1:numstressdofs],3,:)
numdispdofs = NF*2
U1 = reshape(L1[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)

H2 = sol[17:32]
L2 = cellsolvers[2,1].iLLxLH*H2
S2 = -D1*reshape(L2[1:12],3,:)
U2 = reshape(L2[13:20],2,:)

e11 = 0.1/((lambda+2mu) - lambda^2/(lambda+2mu))
e22 = -lambda/(lambda+2mu)*e11
u1 = e11*2
u2 = e22*1
ui = e11*interface_location

testH1 = [0.0 u1  u1  u1 0.0 u1 ui  ui
          0.0 0.0 0.0 u2 u2  u2 0.0 u2]
testH2 = [0.0 u1  0.0 u1 0.0 0.0 ui  ui
          0.0 0.0 u2  u2 0.0 u2  0.0 u2]
testU = copy(coords)
testU[1,:] .*= e11
testU[2,:] .*= e22
testS = zeros(3,NF)
testS[1,:] .= 0.1
@test allapprox(H1,testH1,1e-10)
@test allapprox(H2,testH2,1e-10)
@test allapprox(U1,testU,1e-10)
@test allapprox(U2,testU,1e-10)
@test allapprox(S1,testS,1e-10)
@test allapprox(S2,testS,1e-10)

###########################################################################
# Change to a sloping interface
function plane_distance_function(coords,normal,xc)
    return (coords .- xc)'*normal
end

polyorder = 1
numqp = 2
levelset = InterpolatingPolynomial(1,2,polyorder)
mesh = UniformMesh([0.,0.],[2.,1.],[1,1])
xc = [0.5,0.0]
normal = 1/sqrt(2)*[1,-1]
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->plane_distance_function(x,normal,xc),mesh,levelset.basis)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
NHF = HDGElasticity.number_of_basis_functions(ufs.sbasis)
dofsperelement = 2*NHF

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

lambda,mu = 1.,2.
stabilization = 1e-3
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D1,stabilization)

# Cell 1 phase 1
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,1,[0.,1.],[1.,0.],1,1:4,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,1,1,2,
    2,1:4,dofsperelement)
rhs = -HDGElasticity.linear_form(dgmesh.facescale[2]*[0.1,0.0],ufs.sbasis,ufs.fquads[1,1][2])
HDGElasticity.assemble!(system_rhs,rhs,2,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,3,3,1:4,dofsperelement)

# # Cell 1 phase 2
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,1,[0.,1.],[1.,0.],
    5,5:8,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    2,1,3,6,5:8,dofsperelement)
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,4,
    [1.,0.],[0.,1.],7,5:8,dofsperelement)

# # Interface condition
HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    1,4,1:4,8,5:8,dofsperelement)

K = dropzeros!(HDGElasticity.sparse(system_matrix,32))
R = HDGElasticity.rhs(system_rhs,32)

sol = K\R

dofspercell = 4*dofsperelement
H1 = sol[1:dofspercell]
L1 = cellsolvers[1,1].iLLxLH*H1
numstressdofs = NF*3
S1 = -D1*reshape(L1[1:numstressdofs],3,:)
numdispdofs = NF*2
U1 = reshape(L1[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)

H2 = sol[17:32]
L2 = cellsolvers[2,1].iLLxLH*H2
S2 = -D1*reshape(L2[1:12],3,:)
U2 = reshape(L2[13:20],2,:)

testui = [0.5 1.5
          0.0 1.0]
testui[1,:] .*= e11
testui[2,:] .*= e22
testH1[1:2,7:8] .= testui
testH2[1:2,7:8] .= testui

@test allapprox(H1,testH1,1e-10)
@test allapprox(H2,testH2,1e-10)
@test allapprox(U1,testU,1e-10)
@test allapprox(U2,testU,1e-10)
@test allapprox(S1,testS,1e-10)
@test allapprox(S2,testS,1e-10)

###########################################################################
# Change to a quadratic basis
polyorder = 2
numqp = 4
levelset = InterpolatingPolynomial(1,2,polyorder)
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->distance_function(x,0.5),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
    levelsetcoeffs,levelset)
NHF = HDGElasticity.number_of_basis_functions(ufs.sbasis)
NF = HDGElasticity.number_of_basis_functions(ufs.vbasis)
dofsperelement = 2*NHF
coords = HDGElasticity.nodal_coordinates(mesh,levelset.basis)

system_matrix = HDGElasticity.SystemMatrix()
system_rhs = HDGElasticity.SystemRHS()

lambda,mu = 1.,2.
stabilization = 1e-3
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D1,stabilization)

# Cell 1 phase 1
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,1,[0.,1.],[1.,0.],1,1:4,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,1,1,2,
    2,1:4,dofsperelement)
rhs = -HDGElasticity.linear_form(dgmesh.facescale[2]*[0.1,0.0],ufs.sbasis,ufs.fquads[1,1][2])
HDGElasticity.assemble!(system_rhs,rhs,2,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    1,1,3,3,1:4,dofsperelement)

# Cell 1 phase 2
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,1,[0.,1.],[1.,0.],
    5,5:8,dofsperelement)
HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
    2,1,3,6,5:8,dofsperelement)
HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,4,
    [1.,0.],[0.,1.],7,5:8,dofsperelement)

# Coherent interface condition
HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
    1,4,1:4,8,5:8,dofsperelement)

ndofs = dofsperelement*4*2
K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
R = HDGElasticity.rhs(system_rhs,ndofs)

sol = K\R

dofspercell = 4*dofsperelement
H1 = sol[1:dofspercell]
L1 = cellsolvers[1,1].iLLxLH*H1
numstressdofs = NF*3
S1 = -D1*reshape(L1[1:numstressdofs],3,:)
numdispdofs = NF*2
U1 = reshape(L1[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)

H2 = sol[(dofspercell+1):end]
L2 = cellsolvers[2,1].iLLxLH*H2
S2 = -D1*reshape(L2[1:numstressdofs],3,:)
U2 = reshape(L2[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)

testH1 = [0.0  0.5u1  u1   u1   u1    u1  0.0  0.5u1  u1  ui   ui    ui
          0.0  0.0    0.0  0.0  0.5u2 u2  u2   u2     u2  0.0  0.5u2 u2]
testH2 = [0.0  0.5u1  u1   0.0  0.5u1 u1  0.0  0.0    0.0 ui   ui    ui
          0.0  0.0    0.0  u2   u2    u2  0.0  0.5u2  u2  0.0  0.5u2 u2]

testS = zeros(3,NF)
testS[1,:] .= 0.1
testU = copy(coords)
testU[1,:] .*= e11
testU[2,:] .*= e22
@test allapprox(H1,testH1,1e-10)
@test allapprox(H2,testH2,1e-10)
@test allapprox(testS,S1,1e-10)
@test allapprox(testS,S2,1e-10)
@test allapprox(testU,U1,1e-10)
@test allapprox(testU,U2,1e-8)

###########################################################################
# Check hybrid element numbering

visited = [true,false,false,true]
isactiveface = [true,false,true,false]
facetohelid = zeros(Int,4)
helid = HDGElasticity.assign_cell_hybrid_element_ids!(facetohelid,visited,
    isactiveface,1,4)
@test helid == 2
@test allequal(facetohelid,[0,0,1,0])

mesh = UniformMesh([0.,0.],[2.,1.],[2,1])
levelset = InterpolatingPolynomial(1,2,1)
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->plane_distance_function(x,[1.,0.],[0.5,0.]),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
isactiveface = Matrix{Vector{Bool}}(undef,2,2)
isactiveface[1,1] = [true,true,true,false]
isactiveface[1,2] = [true,true,true,true]
isactiveface[2,1] = [true,false,true,true]
isactiveface[2,2] = zeros(Bool,4)
visited = [zeros(Bool,4) for i = 1:2, j = 1:2]
facetohelid = [zeros(Int,4) for i = 1:2, j = 1:2]
HDGElasticity.assign_cell_hybrid_element_ids!(facetohelid[1,1],
    visited[1,1],isactiveface[1,1],1,4)
@test allequal(facetohelid[1,1],[1,2,3,0])
@test allequal(facetohelid[2,1],[0,0,0,0])
@test allequal(facetohelid[1,2],[0,0,0,0])
@test allequal(facetohelid[2,2],[0,0,0,0])
@test allequal(visited[1,1],[1,1,1,0])
@test allequal(visited[2,1],[0,0,0,0])
@test allequal(visited[1,2],[0,0,0,0])
@test allequal(visited[2,2],[0,0,0,0])

HDGElasticity.assign_neighbor_cell_hybrid_ids!(facetohelid,
    visited,1,1,dgmesh.connectivity,4)
@test allequal(facetohelid[1,1],[1,2,3,0])
@test allequal(facetohelid[2,1],[0,0,0,0])
@test allequal(facetohelid[1,2],[0,0,0,2])
@test allequal(facetohelid[2,2],[0,0,0,0])
@test allequal(visited[1,1],[1,1,1,0])
@test allequal(visited[2,1],[0,0,0,0])
@test allequal(visited[1,2],[0,0,0,1])
@test allequal(visited[2,2],[0,0,0,0])

facetohelid,interfacehelid,helid = HDGElasticity.cell_hybrid_element_ids(dgmesh.cellsign,
    isactiveface,dgmesh.connectivity,4)
@test helid == 12
@test size(facetohelid) == (2,2)
@test allequal(facetohelid[1,1],[1,2,3,0])
@test allequal(facetohelid[2,1],[5,0,6,7])
@test allequal(facetohelid[1,2],[9,10,11,2])
@test allequal(facetohelid[2,2],[0,0,0,0])
@test size(interfacehelid) == (2,2)
@test allequal(interfacehelid,[4 0;8 0])

levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->plane_distance_function(x,1/sqrt(2)*[1.,-1.],[0.5,0.]),mesh,levelset.basis
)
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,4,levelsetcoeffs,levelset)
facetohelid,interfacehelid,helid = HDGElasticity.cell_hybrid_element_ids(dgmesh.cellsign,
    ufs.isactiveface,dgmesh.connectivity,4)
@test helid == 15
@test size(facetohelid) == (2,2)
@test allequal(facetohelid[1,1],[1,2,0,0])
@test allequal(facetohelid[2,1],[4,5,6,7])
@test allequal(facetohelid[1,2],[9,10,11,2])
@test allequal(facetohelid[2,2],[0,0,13,5])
@test size(interfacehelid) == (2,2)
@test allequal(interfacehelid,[3 12;8 14])

hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
@test hybrid_element_numbering.number_of_hybrid_elements == 14
facetohelid = hybrid_element_numbering.facetohelid
@test size(facetohelid) == (2,2)
@test allequal(facetohelid[1,1],[1,2,0,0])
@test allequal(facetohelid[2,1],[4,5,6,7])
@test allequal(facetohelid[1,2],[9,10,11,2])
@test allequal(facetohelid[2,2],[0,0,13,5])
interfacehelid = hybrid_element_numbering.interfacehelid
@test size(interfacehelid) == (2,2)
@test allequal(interfacehelid,[3 12;8 14])

###########################################################################
# Change to a curved interface
# function circle_distance_function(coords,xc,r)
#     dist = [r-norm(coords[:,i]-xc) for i = 1:size(coords)[2]]
# end

# polyorder = 2
# numqp = 6
# levelset = InterpolatingPolynomial(1,2,polyorder)
# circle_center = [2.5,0.5]
# radius = 1.5
# levelsetcoeffs = HDGElasticity.levelset_coefficients(
#     x->circle_distance_function(x,circle_center,radius),mesh,levelset.basis
# )
# dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
# ufs = HDGElasticity.UniformFunctionSpace(dgmesh,polyorder,numqp,
#     levelsetcoeffs,levelset)
# NHF = HDGElasticity.number_of_basis_functions(ufs.sbasis)
# NF = HDGElasticity.number_of_basis_functions(ufs.vbasis)
# dofsperelement = 2*NHF
# coords = HDGElasticity.nodal_coordinates(mesh,levelset.basis)
#
# system_matrix = HDGElasticity.SystemMatrix()
# system_rhs = HDGElasticity.SystemRHS()
#
# lambda,mu = 1.,2.
# stabilization = 1e-3
# D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(lambda,mu))
# cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D1,stabilization)
#
# # Cell 1 phase 1
# HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,
#     1,1,1,[0.,1.],[1.,0.],1,1:4,dofsperelement)
# HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,1,1,2,
#     2,1:4,dofsperelement)
# rhs = -HDGElasticity.linear_form(dgmesh.facescale[2]*[0.1,0.0],ufs.sbasis,ufs.fquads[1,1][2])
# HDGElasticity.assemble!(system_rhs,rhs,2,dofsperelement)
# HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
#     1,1,3,3,1:4,dofsperelement)
#
# # Cell 1 phase 2
# HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,1,[0.,1.],[1.,0.],
#     5,5:8,dofsperelement)
# HDGElasticity.assemble_traction_face!(system_matrix,dgmesh,ufs,cellsolvers,
#     2,1,3,6,5:8,dofsperelement)
# HDGElasticity.assemble_mixed_face!(system_matrix,dgmesh,ufs,cellsolvers,2,1,4,
#     [1.,0.],[0.,1.],7,5:8,dofsperelement)
#
# # Coherent interface condition
# HDGElasticity.assemble_coherent_interface!(system_matrix,dgmesh,ufs,cellsolvers,
#     1,4,1:4,8,5:8,dofsperelement)
#
# ndofs = dofsperelement*4*2
# K = dropzeros!(HDGElasticity.sparse(system_matrix,ndofs))
# R = HDGElasticity.rhs(system_rhs,ndofs)
#
# sol = K\R
#
# dofspercell = 4*dofsperelement
# H1 = sol[1:dofspercell]
# L1 = cellsolvers[1,1].iLLxLH*H1
# numstressdofs = NF*3
# S1 = -D1*reshape(L1[1:numstressdofs],3,:)
# numdispdofs = NF*2
# U1 = reshape(L1[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)
#
# H2 = sol[(dofspercell+1):end]
# L2 = cellsolvers[2,1].iLLxLH*H2
# S2 = -D1*reshape(L2[1:numstressdofs],3,:)
# U2 = reshape(L2[(numstressdofs+1):(numstressdofs+numdispdofs)],2,:)
#
# testH1 = [0.0  0.5u1  u1   u1   u1    u1  0.0  0.5u1  u1  ui   ui    ui
#           0.0  0.0    0.0  0.0  0.5u2 u2  u2   u2     u2  0.0  0.5u2 u2]
# testH2 = [0.0  0.5u1  u1   0.0  0.5u1 u1  0.0  0.0    0.0 ui   ui    ui
#           0.0  0.0    0.0  u2   u2    u2  0.0  0.5u2  u2  0.0  0.5u2 u2]


text:
export AbstractNDSparse, NDSparse, ndsparse

abstract type AbstractNDSparse end

mutable struct NDSparse{T, D<:Tuple, C<:Columns, V<:AbstractVector} <: AbstractNDSparse
    index::C
    data::V
    _table::NextTable

    index_buffer::C
    data_buffer::V
end

function NextTable(nds::NDSparse; kwargs...)
    convert(NextTable, nds.index, nds.data; kwargs...)
end

convert(::Type{NextTable}, nd::NDSparse) = NextTable(nd)

Base.@deprecate_binding IndexedTable NDSparse

# optional, non-exported name
Base.@deprecate_binding Table NDSparse

"""
`ndsparse(indices, data; agg, presorted, copy, chunks)`

Construct an NDSparse array with the given indices and data. Each vector in `indices` represents the index values for one dimension. On construction, the indices and data are sorted in lexicographic order of the indices.

# Arguments:

* `agg::Function`: If `indices` contains duplicate entries, the corresponding data items are reduced using this 2-argument function.
* `presorted::Bool`: If true, the indices are assumed to already be sorted and no sorting is done.
* `copy::Bool`: If true, the storage for the new array will not be shared with the passed indices and data. If false (the default), the passed arrays will be copied only if necessary for sorting. The only way to guarantee sharing of data is to pass `presorted=true`.
* `chunks::Integer`: distribute the table into `chunks` (Integer) chunks (a safe bet is nworkers()). Not distributed by default. See [Distributed](@distributed) docs.

# Examples:

1-dimensional NDSparse can be constructed with a single array as index.
```jldoctest ndsparse
julia> x = ndsparse(["a","b"],[3,4])
1-d NDSparse with 2 values (Int64):
1   │
────┼──
"a" │ 3
"b" │ 4

julia> keytype(x), eltype(x)
(Tuple{String}, Int64)

```

A dimension will be named if constructed with a named tuple of columns as index.
```jldoctest ndsparse
julia> x = ndsparse(@NT(date=Date.(2014:2017)), [4:7;])
1-d NDSparse with 4 values (Int64):
date       │
───────────┼──
2014-01-01 │ 4
2015-01-01 │ 5
2016-01-01 │ 6
2017-01-01 │ 7

```

```jldoctest ndsparse
julia> x[Date("2015-01-01")]
5

julia> keytype(x), eltype(x)
(Tuple{Date}, Int64)

```

Multi-dimensional `NDSparse` can be constructed by passing a tuple of index columns:

```jldoctest ndsparse
julia> x = ndsparse((["a","b"],[3,4]), [5,6])
2-d NDSparse with 2 values (Int64):
1    2 │
───────┼──
"a"  3 │ 5
"b"  4 │ 6

julia> keytype(x), eltype(x)
(Tuple{String,Int64}, Int64)

julia> x["a", 3]
5
```

The data itself can also contain tuples (these are stored in columnar format, just like in `table`.)

```jldoctest ndsparse
julia> x = ndsparse((["a","b"],[3,4]), ([5,6], [7.,8.]))
2-d NDSparse with 2 values (2-tuples):
1    2 │ 3  4
───────┼───────
"a"  3 │ 5  7.0
"b"  4 │ 6  8.0

julia> x = ndsparse(@NT(x=["a","a","b"],y=[3,4,4]),
                    @NT(p=[5,6,7], q=[8.,9.,10.]))
2-d NDSparse with 3 values (2 field named tuples):
x    y │ p  q
───────┼────────
"a"  3 │ 5  8.0
"a"  4 │ 6  9.0
"b"  4 │ 7  10.0

julia> keytype(x), eltype(x)
(Tuple{String,Int64}, NamedTuples._NT_p_q{Int64,Float64})

julia> x["a", :]
2-d NDSparse with 2 values (2 field named tuples):
x    y │ p  q
───────┼───────
"a"  3 │ 5  8.0
"a"  4 │ 6  9.0

```

Passing a `chunks` option to `ndsparse`, or constructing with a distributed array will cause the result to be distributed. Use `distribute` function to distribute an array.

```jldoctest ndsparse
julia> x = ndsparse(@NT(date=Date.(2014:2017)), [4:7.;], chunks=2)
1-d Distributed NDSparse with 4 values (Float64) in 2 chunks:
date       │
───────────┼────
2014-01-01 │ 4.0
2015-01-01 │ 5.0
2016-01-01 │ 6.0
2017-01-01 │ 7.0

julia> x = ndsparse(@NT(date=Date.(2014:2017)), distribute([4:7.0;], 2))
1-d Distributed NDSparse with 4 values (Float64) in 2 chunks:
date       │
───────────┼────
2014-01-01 │ 4.0
2015-01-01 │ 5.0
2016-01-01 │ 6.0
2017-01-01 │ 7.0
```

Distribution is done to match the first distributed column from left to right. Specify `chunks` to override this.
"""
function ndsparse end

function ndsparse(I::Tup, d::Union{Tup, AbstractVector};
                  chunks=nothing, kwargs...)
    if chunks !== nothing
        impl = Val{:distributed}()
    else
        impl = _impl(astuple(I)...)
        if impl === Val{:serial}()
            impl = isa(d, Tup) ?
                _impl(impl, astuple(d)...) : _impl(d)
        end
    end
    ndsparse(impl, I, d; chunks=chunks, kwargs...)
end

function ndsparse(::Val{:serial}, ks::Tup, vs::Union{Tup, AbstractVector};
                  agg=nothing, presorted=false,
                  chunks=nothing, copy=true)

    I = rows(ks)
    d = vs isa Tup ? Columns(vs) : vs

   #if !isempty(filter(x->!isa(x, Int),
   #                   intersect(colnames(I), colnames(d))))
   #    error("All column names, including index and data columns, must be distinct")
   #end
    length(I) == length(d) || error("index and data must have the same number of elements")

    if !presorted && !issorted(I)
        p = sortperm(I)
        I = I[p]
        d = d[p]
    elseif copy
        if agg !== nothing
            iter = GroupReduce(agg, I, d, Base.OneTo(length(I)))
            I, d = collect_columns(iter).columns
            agg = nothing
        else
            I = Base.copy(I)
            d = Base.copy(d)
        end
    end
    stripnames(x) = isa(x, Columns) ? rows(astuple(columns(x))) : rows((x,))
    _table = convert(NextTable, I, stripnames(d); presorted=true, copy=false)
    nd = NDSparse{eltype(d),astuple(eltype(I)),typeof(I),typeof(d)}(
        I, d, _table, similar(I,0), similar(d,0)
    )
    agg===nothing || aggregate!(agg, nd)
    return nd
end

function ndsparse(x::AbstractVector, y; kwargs...)
    ndsparse((x,), y; kwargs...)
end

function ndsparse(x::Tup, y::Columns; kwargs...)
    ndsparse(x, columns(y); kwargs...)
end

function ndsparse(x::Columns, y::AbstractVector; kwargs...)
    ndsparse(columns(x), y; kwargs...)
end

ndsparse(c::Columns{<:Pair}; kwargs...) =
    convert(NDSparse, c.columns.first, c.columns.second; kwargs...)

# backwards compat
NDSparse(idx::Columns, data; kwargs...) = ndsparse(idx, data; kwargs...)

# TableLike API
Base.@pure function colnames(t::NDSparse)
    dnames = colnames(t.data)
    if all(x->isa(x, Integer), dnames)
        dnames = map(x->x+ncols(t.index), dnames)
    end
    vcat(colnames(t.index), dnames)
end

columns(nd::NDSparse) = concat_tup(columns(nd.index), columns(nd.data))

# IndexedTableLike API

permcache(t::NDSparse) = permcache(t._table)
cacheperm!(t::NDSparse, p) = cacheperm!(t._table, p)

"""
    pkeynames(t::NDSparse)

Names of the primary key columns in `t`.

# Example

```jldoctest

julia> x = ndsparse([1,2],[3,4])
1-d NDSparse with 2 values (Int64):
1 │
──┼──
1 │ 3
2 │ 4

julia> pkeynames(x)
(1,)

```
"""
pkeynames(t::NDSparse) = (dimlabels(t)...)

# For an NDSparse, valuenames is either a tuple of fieldnames or a
# single name for scalar values
function valuenames(t::NDSparse)
    if isa(values(t), Columns)
        T = eltype(values(t))
        ((ndims(t) + (1:nfields(eltype(values(t)))))...)
    else
        ndims(t) + 1
    end
end

"""
`NDSparse(columns...; names=Symbol[...], kwargs...)`

Construct an NDSparse array from columns. The last argument is the data column, and the rest are index columns. The `names` keyword argument optionally specifies names for the index columns (dimensions).
"""
function NDSparse(columns...; names=nothing, rest...)
    keys, data = columns[1:end-1], columns[end]
    ndsparse(Columns(keys..., names=names), data; rest...)
end

similar(t::NDSparse) = NDSparse(similar(t.index, 0), similar(t.data, 0))

function copy(t::NDSparse)
    flush!(t)
    NDSparse(copy(t.index), copy(t.data), presorted=true)
end

function (==)(a::NDSparse, b::NDSparse)
    flush!(a); flush!(b)
    return a.index == b.index && a.data == b.data
end

function empty!(t::NDSparse)
    empty!(t.index)
    empty!(t.data)
    empty!(t.index_buffer)
    empty!(t.data_buffer)
    return t
end

_convert(::Type{<:Tuple}, tup::Tuple) = tup
_convert{T<:NamedTuple}(::Type{T}, tup::Tuple) = T(tup...)
convertkey(t::NDSparse{V,K,I}, tup::Tuple) where {V,K,I} = _convert(eltype(I), tup)

ndims(t::NDSparse) = length(t.index.columns)
length(t::NDSparse) = (flush!(t);length(t.index))
eltype{T,D,C,V}(::Type{NDSparse{T,D,C,V}}) = T
Base.keytype{T,D,C,V}(::Type{NDSparse{T,D,C,V}}) = D
Base.keytype(x::NDSparse) = keytype(typeof(x))
dimlabels{T,D,C,V}(::Type{NDSparse{T,D,C,V}}) = fieldnames(eltype(C))

# Generic ndsparse constructor that also works with distributed
# arrays in JuliaDB

Base.@deprecate itable(x, y) ndsparse(x, y)

# Keys and Values iterators

keys(t::NDSparse) = t.index
"""
`keys(x::NDSparse[, select::Selection])`

Get the keys of an `NDSparse` object. Same as [`rows`](@ref) but acts only on the index columns of the `NDSparse`.
"""
keys(t::NDSparse, which...) = rows(keys(t), which...)

# works for both NextTable and NDSparse
pkeys(t::NDSparse, which...) = keys(t, which...)

values(t::NDSparse) = t.data
"""
`values(x::NDSparse[, select::Selection])`

Get the values of an `NDSparse` object. Same as [`rows`](@ref) but acts only on the value columns of the `NDSparse`.
"""
function values(t::NDSparse, which...)
    if values(t) isa Columns
        rows(values(t), which...)
    else
        if which[1] != 1
            error("column $which not found")
        end
        values(t)
    end
end

## Some array-like API

"""
`dimlabels(t::NDSparse)`

Returns an array of integers or symbols giving the labels for the dimensions of `t`.
`ndims(t) == length(dimlabels(t))`.
"""
dimlabels(t::NDSparse) = dimlabels(typeof(t))

start(a::NDSparse) = start(a.data)
next(a::NDSparse, st) = next(a.data, st)
done(a::NDSparse, st) = done(a.data, st)

function permutedims(t::NDSparse, p::AbstractVector)
    if !(length(p) == ndims(t) && isperm(p))
        throw(ArgumentError("argument to permutedims must be a valid permutation"))
    end
    flush!(t)
    NDSparse(Columns(t.index.columns[p]), t.data, copy=true)
end

# showing

import Base.show
function show(io::IO, t::NDSparse{T,D}) where {T,D}
    flush!(t)
    if !(values(t) isa Columns)
        cnames = colnames(keys(t))
        eltypeheader = "$(eltype(t))"
    else
        cnames = colnames(t)
        nf = nfields(eltype(t))
        if eltype(t) <: NamedTuple
            eltypeheader = "$(nf) field named tuples"
        else
            eltypeheader = "$(nf)-tuples"
        end
    end
    header = "$(ndims(t))-d NDSparse with $(length(t)) values (" * eltypeheader * "):"
    showtable(io, t; header=header,
              cnames=cnames, divider=length(columns(keys(t))))
end

import Base: @md_str

function showmeta(io, t::NDSparse, cnames)
    nc = length(columns(t))
    nidx = length(columns(keys(t)))
    nkeys = length(columns(values(t)))

    print(io,"    ")
    with_output_format(:underline, println, io, "Dimensions")
    metat = Columns(([1:nidx;], [Text(get(cnames, i, "<noname>")) for i in 1:nidx],
                     eltype.([columns(keys(t))...])))
    showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true)
    print(io,"\n    ")
    with_output_format(:underline, println, io, "Values")
    if isa(values(t), Columns)
        metat = Columns(([nidx+1:nkeys+nidx;], [Text(get(cnames, i, "<noname>")) for i in nidx+1:nkeys+nidx],
                         eltype.(Any[columns(values(t))...])))
        showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true)
    else
        show(io, eltype(values(t)))
    end
end

abstract type SerializedNDSparse end

function serialize(s::AbstractSerializer, x::NDSparse)
    flush!(x)
    Base.Serializer.serialize_type(s, SerializedNDSparse)
    serialize(s, x.index)
    serialize(s, x.data)
end

function deserialize(s::AbstractSerializer, ::Type{SerializedNDSparse})
    I = deserialize(s)
    d = deserialize(s)
    NDSparse(I, d, presorted=true)
end

@noinline convert(::Type{NDSparse}, ks::ANY, vs::ANY; kwargs...) = ndsparse(ks, vs; kwargs...)
@noinline convert(T::Type{NDSparse}, c::Columns{<:Pair}; kwargs...) = convert(T, c.columns.first, c.columns.second; kwargs...)

# map and convert

"""
    map(f, x::NDSparse; select)

Apply `f` to every data value in `x`. `select` selects fields
passed to `f`. By default, the data values are selected.

If the return value of `f` is a tuple or named tuple the result
will contain many data columns.

# Examples

```jldoctest
julia> x = ndsparse(@NT(t=[0.01, 0.05]), @NT(x=[1,2], y=[3,4]))
1-d NDSparse with 2 values (2 field named tuples):
t    │ x  y
─────┼─────
0.01 │ 1  3
0.05 │ 2  4

julia> manh = map(row->row.x + row.y, x)
1-d NDSparse with 2 values (Int64):
t    │
─────┼──
0.01 │ 4
0.05 │ 6

julia> vx = map(row->row.x/row.t, x, select=(:t,:x))
1-d NDSparse with 2 values (Float64):
t    │
─────┼──────
0.01 │ 100.0
0.05 │ 40.0

julia> polar = map(p->@NT(r=hypot(p.x + p.y), θ=atan2(p.y, p.x)), x)
1-d NDSparse with 2 values (2 field named tuples):
t    │ r    θ
─────┼─────────────
0.01 │ 4.0  1.24905
0.05 │ 6.0  1.10715

julia> map(sin, polar, select=:θ)
1-d NDSparse with 2 values (Float64):
t    │
─────┼─────────
0.01 │ 0.948683
0.05 │ 0.894427

```
"""
function map(f, x::NDSparse; select=x.data)
    ndsparse(copy(x.index), map_rows(f, rows(x, select)),
             presorted=true, copy=false)
end

# """
# `columns(x::NDSparse, names...)`
#
# Given an NDSparse array with multiple data columns (its data vector is a `Columns` object), return a
# new array with the specified subset of data columns. Data is shared with the original array.
# """
# columns(x::NDSparse, which...) = NDSparse(x.index, Columns(x.data.columns[[which...]]), presorted=true)

#columns(x::NDSparse, which) = NDSparse(x.index, x.data.columns[which], presorted=true)

#column(x::NDSparse, which) = columns(x, which)

# NDSparse uses lex order, Base arrays use colex order, so we need to
# reorder the data. transpose and permutedims are used for this.
convert(::Type{NDSparse}, m::SparseMatrixCSC) = NDSparse(findnz(m.')[[2,1,3]]..., presorted=true)

# special method to allow selection on
# ndsparse with repeating names in keys and values
function column(x::NDSparse, which::Integer)
    @assert which > 0
    if which <= ndims(x)
        keys(x, which)
    else
        values(x, which-ndims(x))
    end
end

function convert{T}(::Type{NDSparse}, a::AbstractArray{T})
    n = length(a)
    nd = ndims(a)
    a = permutedims(a, [nd:-1:1;])
    data = reshape(a, (n,))
    idxs = [ Vector{Int}(n) for i = 1:nd ]
    i = 1
    for I in CartesianRange(size(a))
        for j = 1:nd
            idxs[j][i] = I[j]
        end
        i += 1
    end
    NDSparse(Columns(reverse(idxs)...), data, presorted=true)
end

# aggregation

"""
`aggregate!(f::Function, arr::NDSparse)`

Combine adjacent rows with equal indices using the given 2-argument reduction function,
in place.
"""
function aggregate!(f, x::NDSparse)
    idxs, data = x.index, x.data
    n = length(idxs)
    newlen = 0
    i1 = 1
    while i1 <= n
        val = data[i1]
        i = i1+1
        while i <= n && roweq(idxs, i, i1)
            val = f(val, data[i])
            i += 1
        end
        newlen += 1
        if newlen != i1
            copyrow!(idxs, newlen, i1)
        end
        data[newlen] = val
        i1 = i
    end
    resize!(idxs, newlen)
    resize!(data, newlen)
    x
end

function subtable(x::NDSparse, idx; presorted=true)
    ndsparse(keys(x)[idx], values(x)[idx], presorted=presorted, copy=false)
end


text:
@doc raw"""
    SigmoidKernel([a=1 [,c=1]])

The Sigmoid Kernel is given by:
```
    κ(x,y) = tanh(a⋅xᵀy + c)
```

# Examples

```jldoctest; setup = :(using MLKernels)
julia> SigmoidKernel()
SigmoidKernel{Float64}(1.0,1.0)

julia> SigmoidKernel(0.5f0)
SigmoidKernel{Float32}(0.5,1.0)

julia> SigmoidKernel(0.5f0, 0.5)
SigmoidKernel{Float64}(0.5,0.5)
```
"""
struct SigmoidKernel{T<:AbstractFloat} <: Kernel{T}
    a::T
    c::T
    function SigmoidKernel{T}(a::Real=T(1), c::Real=T(1)) where {T<:AbstractFloat}
        @check_args(SigmoidKernel, a, a >  zero(T), "a > 0")
        @check_args(SigmoidKernel, c, c >= zero(T), "c ≧ 0")
        return new{T}(a, c)
    end
end
function SigmoidKernel(a::T₁=1.0, c::T₂=T₁(1)) where {T₁<:Real,T₂<:Real}
    SigmoidKernel{promote_float(T₁,T₂)}(a,c)
end

@inline basefunction(::SigmoidKernel) = ScalarProduct()

@inline kappa(κ::SigmoidKernel{T}, xᵀy::T) where {T} = tanh(κ.a*xᵀy + κ.c)

function convert(::Type{K}, κ::SigmoidKernel) where {K>:SigmoidKernel{T}} where T
    return SigmoidKernel{T}(κ.a, κ.c)
end