text:
"""
Parent Centric Crossover (PCX).

See
    Deb, K., Anand, A., and Joshi, D., "A Computationally Efficient Evolutionary Algorithm for Real-Parameter Optimization," Evolutionary Computation, vol. 10, no. 4, pp. 371-395, 2002.
"""
struct ParentCentricCrossover{NP} <: CrossoverOperator{NP,1}
    ζ::Float64 # sd for the orthogonal directions defined by 2nd,3rd etc parents
    η::Float64 # sd for the direction of the 1st parent

    function ParentCentricCrossover{NP}(ζ::Number, η::Number) where NP
        ζ > 0 || throw(ArgumentError("ζ must be positive"))
        η > 0 || throw(ArgumentError("η must be non-negative"))
        new{NP}(ζ, η)
    end
    ParentCentricCrossover{NP}(params::Parameters) where NP =
        new{NP}(params[:PCX_ζ], params[:PCX_η])
end

const PCX_DefaultOptions = ParamsDict(
    :PCX_ζ => 0.1,
    :PCX_η => 0.5
)

function apply!(xover::ParentCentricCrossover,
                target::AbstractIndividual, targetIndex::Int,
                pop, parentIndices)
    @assert length(parentIndices) == numparents(xover)

    center = mean(viewer(pop, parentIndices), dims=2)
    # project other parents vectors to the subspace
    # orthogonal to the 1st (central) parent
    parent_1 = viewer(pop, parentIndices[1]) .- dropdims(center, dims=2)
    other_parents_centered = viewer(pop, parentIndices[2:end]) .- center
    if (parent_1_sqnorm = sum(abs2, parent_1)) > 0.0
        dotprods = transpose(parent_1) * other_parents_centered
        dotprods ./= parent_1_sqnorm
        other_parents_centered .-= dotprods .* parent_1
    end
    sd = mean(sqrt.(sum(abs2, other_parents_centered, dims=2)))
    if sd > 1E-8
        svdf = svd(other_parents_centered, full=false)
        svals_norm = norm(svdf.S)
        svals = svdf.S * sd * xover.ζ / svals_norm
        svals .*= randn(length(svals))
        mul!(target, svdf.U, svals)
    else # degenerated
        fill!(target, zero(eltype(target)))
    end

    selParent = viewer(pop, parentIndices[1])
    selScale = randn() * xover.η
    @inbounds target .+= selParent .+ selScale .* parent_1
    return target
end


text:
module Utils

# Dependancies:
using CSV
using DataFrames

"""
    ReadLaborTracker(fName::String)

Function to read ACTUAL_PLANNED hours from SAP.
Report should be saved in CSV format.
Rows with missing names of resources are ignored.

# Arguments
- `fName::String`: filename of the report

# Returns
- `df::DataFrame`: DataFrame of the report

# Example:

```julia
    df = ReadLaborTracker("C:\\Users\\james\\Desktop\\labor_tracker.csv")
```

"""
function ReadLaborTracker(fName::String)

    df = CSV.read(fName, DataFrame)
    df = dropmissing(df, :"Employee Name");

    return df;
end

"""
    ReadCostTracker(fName::String)

Function to read cost tracker report from SAP.
On SAP Report, you would input all project numbers which you may be interested. The report will pull them all.

# Arguments
- `fName::String`: filename of the report

# Returns
- `df::DataFrame`: DataFrame of the report

"""
function ReadCostTracker(fName::String)
    # This function takes report from SAP and saved as CSV file
    # Report from SAP to generate report: Cost Tracker - Sub Product Line

    df = CSV.read(fName, DataFrame)
    df = dropmissing(df, :"Sub - product line");

    if (df.Quoted[1] |> typeof |> supertype == InlineString)
        sum(occursin.("," , df.Quoted)) > 0 ? df.Quoted = replace.(df.Quoted, "," => "") : ();
        df.Quoted = parse.(Float64, df.Quoted);
    end

    if (df.Projected[1] |> typeof |> supertype == InlineString)
        sum(occursin.("," , df.Projected)) > 0 ? df.Projected = replace.(df.Projected, "," => "") : ();
        df.Projected = parse.(Float64, df.Projected);
    end

    if (df.Actual[1] |> typeof |> supertype == InlineString)
        sum(occursin.("," , df.Actual)) > 0 ? df.Actual = replace.(df.Actual, "," => "") : ();
        df.Actual = parse.(Float64, df.Actual);
    end

    if (df.:"Antic."[1] |> typeof |> supertype == InlineString)
        sum(occursin.("," , df.:"Antic.")) > 0 ? df.:"Antic." = replace.(df.:"Antic.", "," => "") : ();
        df.:"Antic." = parse.(Float64, df.:"Antic.");
    end

    if (df.Var[1] |> typeof |> supertype == InlineString)
        sum(occursin.("," , df.Var)) > 0 ? df.Var = replace.(df.Var, "," => "") : ();
        df.Var = parse.(Float64, df.Var);
    end
    # convert df.:"Project Definition"  Column from Int to String type
    df.:"Project Definition" = string.(df.:"Project Definition");

    return df;

end

"""
    ReadAvailHours(fName::String)

Function to read Available hours on a given month from SAP.
Report should be saved in CSV format.
YLF_UTIL_DEPT is used to pull available hours from.
Since SAP report only goes fwd 12 months, spreadhsheet can be
edited to include avail hours beyond 12 months.

# Arguments
- `fName::String`: filename of the report

# Returns
- `df::DataFrame`: DataFrame of the report

# Example:

```julia
    df = ReadAvailHours("C:\\Users\\james\\Desktop\\labor_tracker.csv")
```

"""
function ReadAvailHours(fName::String)

    df = CSV.read(fName, DataFrame)

    return df;
end

"""
    vec_to_dic(v_vals::Vector, v_keys::Vector)

Function to convert vector pair to Dictionary format

# Arguments
- `v_vals::Vector`: vector of values
- `v_keys::Vector`: vector of keys

# Returns
- `dic::Dictionary`: Dictionary

"""
function vec_to_dic(v_keys::Vector, v_vals::Vector)

    @assert length(v_vals) == length(v_keys) "Vector lengths must match!"

    dic = Dict();
    for i = 1:length(v_vals)
        dic[v_keys[i]] = v_vals[i];
    end

    return dic;
end

"""
    dic_to_vec(dic::Dict)
Function to return Dictionary key:value pair  as two vectors: keys and values

# Arguments
- `dic::Dict`: Dictionary

# Returns
- `v_keys::Vector`: vector of keys
- `v_vals::Vector`: vector of values
"""
function dic_to_vec(dic::Dict)

    v_keys = [];
    v_vals = [];
    for k in keys(dic), v in values(dic)
        push!(v_keys, k);
        push!(v_vals, v);
    end

    return v_keys, v_vals;
end

end

ext: jl


text:
@doc raw"""
    exp(M::AbstractManifold, p, X)
    exp(M::AbstractManifold, p, X, t::Real = 1)

Compute the exponential map of tangent vector `X`, optionally scaled by `t`,  at point `p`
from the manifold [`AbstractManifold`](@ref) `M`, i.e.

```math
\exp_p X = γ_{p,X}(1),
```
where ``γ_{p,X}`` is the unique geodesic starting in ``γ(0)=p`` such that ``\dot γ(0) = X``.

See also [`shortest_geodesic`](@ref), [`retract`](@ref).
"""
function exp(M::AbstractManifold, p, X)
    q = allocate_result(M, exp, p, X)
    exp!(M, q, p, X)
    return q
end
exp(M::AbstractManifold, p, X, t::Real) = exp(M, p, t * X)

"""
    exp!(M::AbstractManifold, q, p, X)
    exp!(M::AbstractManifold, q, p, X, t::Real = 1)

Compute the exponential map of tangent vector `X`, optionally scaled by `t`,  at point `p`
from the manifold [`AbstractManifold`](@ref) `M`.
The result is saved to `q`.

See also [`exp`](@ref).
"""
function exp!(M::AbstractManifold, q, p, X)
    return error(manifold_function_not_implemented_message(M, exp!, q, p, X))
end
exp!(M::AbstractManifold, q, p, X, t::Real) = exp!(M, q, p, t * X)

@doc raw"""
    geodesic(M::AbstractManifold, p, X) -> Function

Get the geodesic with initial point `p` and velocity `X` on the [`AbstractManifold`](@ref) `M`.
A geodesic is a curve of zero acceleration. That is for the curve ``γ_{p,X}: I → \mathcal M``,
with ``γ_{p,X}(0) = p`` and ``\dot γ_{p,X}(0) = X`` a geodesic further fulfills

```math
∇_{\dot γ_{p,X}(t)} \dot γ_{p,X}(t) = 0,
```

i.e. the curve is acceleration free with respect to the Riemannian metric.
This yields, that the curve has constant velocity that is locally distance-minimizing.

This function returns a function of (time) `t`.

    geodesic(M::AbstractManifold, p, X, t::Real)
    geodesic(M::AbstractManifold, p, X, T::AbstractVector) -> AbstractVector

Return the point at time `t` or points at times `t` in `T` along the geodesic.
"""
geodesic(M::AbstractManifold, p, X) = t -> exp(M, p, X, t)
geodesic(M::AbstractManifold, p, X, t::Real) = exp(M, p, X, t)
geodesic(M::AbstractManifold, p, X, T::AbstractVector) = map(t -> exp(M, p, X, t), T)

"""
    log(M::AbstractManifold, p, q)

Compute the logarithmic map of point `q` at base point `p` on the [`AbstractManifold`](@ref) `M`.
The logarithmic map is the inverse of the [`exp`](@ref)onential map.
Note that the logarithmic map might not be globally defined.

See also [`inverse_retract`](@ref).
"""
function log(M::AbstractManifold, p, q)
    X = allocate_result(M, log, p, q)
    log!(M, X, p, q)
    return X
end

"""
    log!(M::AbstractManifold, X, p, q)

Compute the logarithmic map of point `q` at base point `p` on the [`AbstractManifold`](@ref) `M`.
The result is saved to `X`.
The logarithmic map is the inverse of the [`exp!`](@ref)onential map.
Note that the logarithmic map might not be globally defined.

see also [`log`](@ref) and [`inverse_retract!`](@ref),
"""
function log!(M::AbstractManifold, X, p, q)
    return error(manifold_function_not_implemented_message(M, log!, X, p, q))
end

@doc raw"""
    shortest_geodesic(M::AbstractManifold, p, q) -> Function

Get a [`geodesic`](@ref) $γ_{p,q}(t)$ whose length is the shortest path between the
points `p`and `q`, where $γ_{p,q}(0)=p$ and $γ_{p,q}(1)=q$. When there are
multiple shortest geodesics, a deterministic choice will be returned.

This function returns a function of time, which may be a `Real` or an `AbstractVector`.

    shortest_geodesic(M::AabstractManifold, p, q, t::Real)
    shortest_geodesic(M::AbstractManifold, p, q, T::AbstractVector) -> AbstractVector

Return the point at time `t` or points at times `t` in `T` along the shortest [`geodesic`](@ref).
"""
shortest_geodesic(M::AbstractManifold, p, q) = geodesic(M, p, log(M, p, q))
shortest_geodesic(M::AbstractManifold, p, q, t::Real) = geodesic(M, p, log(M, p, q), t)
function shortest_geodesic(M::AbstractManifold, p, q, T::AbstractVector)
    return geodesic(M, p, log(M, p, q), T)
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using DirectTrajectoryOptimization
using LinearAlgebra
using Plots

# ## horizon
T = 101

# ## acrobot
nx = 4
nu = 1
nw = 0

function acrobot(x, u, w)
    mass1 = 1.0
    inertia1 = 0.33
    length1 = 1.0
    lengthcom1 = 0.5

    mass2 = 1.0
    inertia2 = 0.33
    length2 = 1.0
    lengthcom2 = 0.5

    gravity = 9.81
    friction1 = 0.1
    friction2 = 0.1

    function M(x, w)
        a = (inertia1 + inertia2 + mass2 * length1 * length1
            + 2.0 * mass2 * length1 * lengthcom2 * cos(x[2]))

        b = inertia2 + mass2 * length1 * lengthcom2 * cos(x[2])

        c = inertia2

       return [a b; b c]
    end

    function τ(x, w)
        a = (-1.0 * mass1 * gravity * lengthcom1 * sin(x[1])
            - mass2 * gravity * (length1 * sin(x[1])
            + lengthcom2 * sin(x[1] + x[2])))

        b = -1.0 * mass2 * gravity * lengthcom2 * sin(x[1] + x[2])

        return [a; b]
    end

    function C(x, w)
        a = -2.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
        b = -1.0 * mass2 * length1 * lengthcom2 * sin(x[2]) * x[4]
        c = mass2 * length1 * lengthcom2 * sin(x[2]) * x[3]
        d = 0.0

        return [a b; c d]
    end

    function B(x, w)
        [0.0; 1.0]
    end

    q = view(x, 1:2)
    v = view(x, 3:4)

    qdd = M(q, w) \ (-1.0 * C(x, w) * v
            + τ(q, w) + B(q, w) * u[1] - [friction1; friction2] .* v)

    return [x[3]; x[4]; qdd[1]; qdd[2]]
end

function midpoint_implicit(y, x, u, w)
    h = 0.05 # timestep
    y - (x + h * acrobot(0.5 * (x + y), u, w))
end

# ## model
dt = Dynamics(midpoint_implicit, nx, nx, nu, nw=nw)
dyn = [dt for t = 1:T-1]
model = DynamicsModel(dyn)

# ## initialization
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [0.0; π; 0.0; 0.0]

# ## objective
ot = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4]) + 0.1 * dot(u, u)
oT = (x, u, w) -> 0.1 * dot(x[3:4], x[3:4])
ct = Cost(ot, nx, nu, nw, [t for t = 1:T-1])
cT = Cost(oT, nx, 0, nw, [T])
obj = [ct, cT]

# ## constraints
x_init = Bound(nx, nu, [1], xl=x1, xu=x1)
x_goal = Bound(nx, 0, [T], xl=xT, xu=xT)
cons = ConstraintSet([x_init, x_goal], [StageConstraint()])

# ## problem
trajopt = TrajectoryOptimizationProblem(obj, model, cons)
s = Solver(trajopt)

# ## initialize
x_interpolation = linear_interpolation(x1, xT, T)
u_guess = [1.0 * randn(nu) for t = 1:T-1]
z0 = zeros(s.p.num_var)
for (t, idx) in enumerate(s.p.trajopt.model.idx.x)
    z0[idx] = x_interpolation[t]
end
for (t, idx) in enumerate(s.p.trajopt.model.idx.u)
    z0[idx] = u_guess[t]
end
initialize!(s, z0)

# ## solve
@time solve!(s)

# ## solution
@show trajopt.x[1]
@show trajopt.x[T]

# ## state
plot(hcat(trajopt.x...)')

# ## control
plot(hcat(trajopt.u[1:end-1]..., trajopt.u[end-1])', linetype = :steppost)


text:
using RecyclePermanent
using RNGPool
using JLD2
using BenchmarkTools

include("../models/gkModel.jl")
include("compare.jl")

@load "data/gkData.jld" ys θ
ys = ys[1:100]
model = makeGKModel(θ, ys)

N = 80*model.n

trials = 1000

setRNGs(12345)

n = 100
N = 50*n

rpio = RPIO{model.Tx}(N, n, 1)

simpleOutput = testEstimator(() -> simpleEstimate(model, rpio), trials)
biasedOutput = testEstimator(() -> biasedEstimate(model, rpio), trials)
rpOutput = testEstimator(() -> rpEstimate(model, rpio), trials)

kmio = KMIO{model.Tx}(N, n)
# ϵ = 1e-2
ϵ = 1e-5
kmEstimate(model, kmio, ϵ)
kmOutput = testEstimator(() -> kmEstimate(model, kmio, ϵ), trials)

@btime rpEstimate($model, $rpio)


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames
using Suppressor

@testset "mixtures EM algorithm tests" begin
    num_mix = 3
    data = DataFrame(BitArray([1 0 1 0 1 0 1 0 1 0;
                                1 1 1 1 1 1 1 1 1 1;
                                0 0 0 0 0 0 0 0 0 0;
                                0 1 1 0 1 0 0 1 0 1]), :auto)
    num_v = num_features(data)

    # init
    @test_throws Exception initial_weights(data, num_mix; alg="")
    w1 = initial_weights(data, num_mix; alg="cluster")
    w2 = initial_weights(data, num_mix; alg="random")
    @test sum(w1) ≈ 1.0
    @test sum(w2) ≈ 1.0
    @test length(w1) == length(w2) == num_mix

    # em
    pc = fully_factorized_circuit(ProbCircuit, num_v)
    estimate_parameters(pc, data; pseudocount=1.0)
    spc = compile(SharedProbCircuit, pc, num_mix)
    values, flows, node2id = satisfies_flows(spc, data)
    ll0 = log_likelihood_per_instance_per_component(spc, data, values, flows, node2id)
    ll = log_likelihood_per_instance(pc, data)
    @test all(ll .≈ ll0)
    w = [0.1 0.2 0.7;]
    ll1, w = one_step_em(spc, data, values, flows, node2id, w;pseudocount=1.0)
    @test sum(ll1) ≈ sum(ll)
    ll2, w = one_step_em(spc, data, values, flows, node2id, w;pseudocount=1.0)
    ll3, w = one_step_em(spc, data, values, flows, node2id, w;pseudocount=1.0)
    @test sum(ll3) > sum(ll2)
    foreach(spc) do n
        if n isa SharedSumNode
            @test all(sum(exp.(n.log_probs), dims=1) .≈ 1.0)
        end
    end
    @test_nowarn learn_circuit_mixture(pc, data; num_mix=1, em_maxiter=2, verbose = false)
    @test_nowarn learn_circuit_mixture(pc, data; num_mix=num_mix, em_maxiter=2, verbose = false)

    @test_nowarn learn_strudel(data; num_mix = num_mix, init_maxiter = 2, em_maxiter = 2, verbose = false)
end


text:
function gaussconditional(μ::Array{Float64,1}, Σ::Array{Float64,2}, a_idx::Array{Int64,1}, b_idx::Array{Int64,1})

    # make sure the given indices do not overlap
    @assert(isempty(intersect(a_idx, b_idx)))

    Σ_aa = Σ[a_idx, a_idx]

    Σ_ab = Σ[a_idx, b_idx]

    Σ_ba = Σ[b_idx, a_idx]

    Σ_bb = Σ[b_idx, b_idx]

    μ_a  = @view μ[a_idx]

    μ_b  = @view μ[a_idx]

    # L = cholesky(Σ_bb).L
    #
    # Q = L \ Σ_ba
    #
    # μ_a_given_b(x) = μ_a + Σ_ab*(L' \ (L \ (x .- μ_b)))
    #
    # Σ_a_given_b    = PDMat(Σ_aa - Q'*Q)

    μ_a_given_b(x) = μ_a + Σ_ab*(Σ_bb \ (x .- μ_b))

     Σ_a_given_b    = PDMat(Σ_aa - Σ_ab*(Σ_bb\Σ_ba))

    return μ_a_given_b, Σ_a_given_b

end

function MvNormalConditional(g::AbstractMvNormal, a_idx::Array{Int64,1}, b_idx::Array{Int64,1})

    μcond, Σcond = gaussconditional(g.μ, g.Σ.mat, a_idx, b_idx)

    x -> MvNormal(μcond(x), Σcond)

end

function MvNormalPartition(g::AbstractMvNormal, keep_idx::Array{Int64,1})

    μnew = g.μ[keep_idx]

    Σnew = PDMat(g.Σ.mat[keep_idx, keep_idx])

    MvNormal(μnew, Σnew)

end

function testMvNormalConditional()

    # let's do a comparison on an example calculated in two different ways

end


text:
@defstruct RandomNormalLayer Layer (
  name :: AbstractString = "random-normal",
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (output_dims :: Vector{Int} = Int[], length(output_dims) >0),
  (eltype :: DataType = Float32),
  (batch_sizes :: Vector{Int} = Int[], (length(batch_sizes)==length(tops))),
)

@characterize_layer(RandomNormalLayer,
  is_source => true
)

type RandomNormalLayerState <: LayerState
    layer :: RandomNormalLayer
    blobs :: Vector{Blob}
    etc        :: Vector{Any}

    RandomNormalLayerState(backend::Backend, layer::RandomNormalLayer) = begin
      blobs = Array{Blob}(length(layer.tops))
      for i = 1:length(blobs)
        dims = tuple(layer.output_dims..., layer.batch_sizes[i])
        blobs[i] = make_blob(backend, layer.eltype, dims...)
    end
    new(layer, blobs, Any[])
  end
end

function setup_etc(backend::CPUBackend, layer::RandomNormalLayer)
  return Any[]
end

function destroy_etc(backend::CPUBackend, state::RandomNormalLayerState)
  # do nothing
end

function setup(backend::Backend, layer::RandomNormalLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  @assert length(inputs) == 0
  state = RandomNormalLayerState(backend, layer)
  state.etc = setup_etc(backend, layer)
  return state
end

function shutdown(backend::Backend, state::RandomNormalLayerState)
  map(destroy, state.blobs)
  destroy_etc(backend, state)
end

function forward(backend::CPUBackend, state::RandomNormalLayerState, inputs::Vector{Blob})
  for i = 1:length(state.blobs)
      randn!(state.blobs[i])
  end
end

function backward(backend::Backend, state::RandomNormalLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
end


text:
using Documenter
using NNFEM
using Markdown
makedocs(sitename="NNFEM", modules=[NNFEM],
pages = Any[
    "index.md",
    "Examples: Inverse Problems"=>["verify.md", "verify_linear.md", "hyperelasticity.md", "verify_function.md"],
    "Examples: Forward Computation"=>["ex_simulation.md", "ex_ad.md"],
    "Manual"=>["representation.md", "solvers.md", "gmsh.md"],
    "api.md"
],
authors = "Kailai Xu and Daniel (Zhengyu) Huang")

deploydocs(
    repo = "github.com/kailaix/NNFEM.jl.git",
)


text:
using Documenter, SeisNetwork

makedocs(;
    modules=[SeisNetwork],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/kura-okubo/SeisNetwork.jl/blob/{commit}{path}#L{line}",
    sitename="SeisNetwork.jl",
    authors="kurama",
    assets=String[],
)

deploydocs(;
    repo="github.com/kura-okubo/SeisNetwork.jl",
)


text:
@testset "n-Queens" begin

    N = 8
    P = 100
    generatePositions = ()->collect(1:N)[randperm(N)]

    # Vector of N cols filled with numbers from 1:N specifying row position
    function nqueens(queens::Vector{Int})
        n = length(queens)
        fitness = 0
        for i=1:(n-1)
            for j=(i+1):n
                k = abs(queens[i] - queens[j])
                if (j-i) == k || k == 0
                    fitness += 1
                end
                # println("$(i),$(queens[i]) <=> $(j),$(queens[j]) : $(fitness)")
            end
        end
        return fitness
    end
    @test nqueens([2,4,1,3]) == 0
    @test nqueens([3,1,2]) == 1

    # Testing: GA solution with various mutations
    for muts in [inversion, insertion, swap2, scramble, shifting]
        result = Evolutionary.optimize(nqueens,
            generatePositions,
            GA(
                populationSize = P,
                selection = susinv,
                crossover = PMX,
                mutation = muts
            ));
        # show(result)
        println("GA:PMX:$(string(muts))(N=$(N), P=$(P)) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
        @test nqueens(Evolutionary.minimizer(result)) == 0
        @test Evolutionary.iterations(result) < 100 # Test early stopping
    end

    # Testing: ES
    for muts in [inversion, insertion, swap2, scramble, shifting]
        for sel in [:plus, :comma]
            μ = 20
            result = Evolutionary.optimize(nqueens, generatePositions,
                ES(mutation = mutationwrapper(muts), μ=μ, ρ=1, λ=P, selection=sel),
                Evolutionary.Options(show_trace=false))
            println("($μ$(sel == :plus ? "+" : ",")$(P))-ES:$(string(muts)) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
            # show(result)
            @test nqueens(Evolutionary.minimizer(result)) == 0
        end
    end

end


text:
function crossoverone(individual1, individual2)
    s = size(individual1, 1)
    crosspoint = rand(1:s)
    individual1[1:crosspoint], individual2[1:crosspoint] = individual2[1:crosspoint], individual1[1:crosspoint]
    return individual1, individual2
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
"Compute a transmission coefficient for a Gaussian absorption line with given depth.  Optionally, truncate absorption at limit_line_effect."
function std_gaussian_line(x::Number, depth::Number; limit_line_effect = Inf)
   return abs(x)>limit_line_effect ? one(x) : one(x)-depth*exp(-0.5*x^2)
end

"Compute spectrum for an Gaussian absorption line give it's location, width and depth"
function absorption_line(x::Number, location::Number, width::Number, depth::Number; limit_line_effect = Inf)
    std_gaussian_line((x-location)/width, depth, limit_line_effect=limit_line_effect)
end

"Compute product of absorption lines at one wavelength given a list of locations, widths and depths"
function absorption_lines(x::T1, locations::A2, widths::A3, depths::A4; limit_line_effect = Inf) where { T1<:Number, T2<:Number, T3<:Number, T4<:Number, A2<:AbstractVector{T2}, A3<:AbstractVector{T3}, A4<:AbstractVector{T4} }
    @assert length(locations) == length(widths) == length(depths) >= 1
    @inbounds trans = absorption_line(x, locations[1], widths[1], depths[1], limit_line_effect = limit_line_effect )
    for i in 2:length(locations)
        @inbounds trans = trans * absorption_line(x, locations[i],widths[i], depths[i], limit_line_effect = limit_line_effect )
    end
    return trans
end

function absorption_lines(x::A1, locations::A2, widths::A3, depths::A4; limit_line_effect = Inf) where { T1<:Number, T2<:Number, T3<:Number, T4<:Number, A1<:AbstractArray{T1}, A2<:AbstractVector{T2}, A3<:AbstractVector{T3}, A4<:AbstractVector{T4} }
    @assert length(locations) == length(widths) == length(depths) >= 1
    trans = absorption_line.(x, locations[1], widths[1], depths[1], limit_line_effect = limit_line_effect )
    @inbounds for i in 2:length(locations)
        trans = trans .* absorption_line.(x, locations[i],widths[i], depths[i], limit_line_effect = limit_line_effect )
    end
    return trans
end


text:
using BayesTests
using Base.Test

binomBF = BinomTest(1, 1)
samples = posterior(binomBF; iter = 10000)
credible = hpd(samples)

@test binomBF.BF == 1.0
@test_approx_eq_eps credible[1] .167 .1
@test_approx_eq_eps credible[2] .986 .1

update(binomBF, 3, 3)

corBF = CorrelationTest(1, 1; a = 1)
corBF2 = CorrelationTest(3, 1; a = 1)
@test isnan(corBF.BF)
@test !isnan(corBF2.BF)

x1 = collect(1:10)
x2 = collect(11:20)

ttestBF = TTest(x1, x2)
update(ttestBF, x1, x2)

ttestBF = TTest(x1, x2; typ = :twosample)
update(ttestBF, x1, x2)


text:
"""
    Base.read!(io::IO, model::FileFormats.MOF.Model)

Read `io` in the MathOptFormat file format and store the result in `model`.
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read model from file as destination model is not empty.")
    end
    options = get_options(model)
    if options.validate
        validate(io)
    end
    object = JSON.parse(io; dicttype=Object)
    file_version = _parse_mof_version(object["version"])
    if file_version > VERSION
        error("Sorry, the file $(filename) can't be read because this library" *
              " supports v$(VERSION) of MathOptFormat, but the file you are " *
              "trying to read is v$(file_version).")
    end
    name_map = read_variables(model, object)
    read_objective(model, object, name_map)
    read_constraints(model, object, name_map)
    return
end

function read_variables(model::Model, object::Object)
    indices = MOI.add_variables(model, length(object["variables"]))
    name_map = Dict{String, MOI.VariableIndex}()
    for (index, variable) in zip(indices, object["variables"])
        if !haskey(variable, "name")
            error("Variable is missing a `name` field.")
        end
        name = variable["name"]
        if name == ""
            error("Variable name is `\"\"`. MathOptFormat variable names must" *
                  " be unique and non-empty.")
        end
        MOI.set(model, MOI.VariableName(), index, name)
        name_map[name] = index
    end
    return name_map
end

function read_objective(
    model::Model, object::Object, name_map::Dict{String, MOI.VariableIndex}
)
    obj = object["objective"]
    sense = read_objective_sense(obj["sense"])
    MOI.set(model, MOI.ObjectiveSense(), sense)
    if sense == MOI.FEASIBILITY_SENSE
        return
    end
    func = function_to_moi(obj["function"], model, name_map)
    MOI.set(model, MOI.ObjectiveFunction{typeof(func)}(), func)
    return
end

function read_constraints(model::Model, object::Object,
                          name_map::Dict{String, MOI.VariableIndex})
    for constraint in object["constraints"]
        foo = function_to_moi(constraint["function"], model, name_map)
        set = set_to_moi(constraint["set"], model, name_map)
        index = MOI.add_constraint(model, foo, set)
        if haskey(constraint, "name")
            MOI.set(model, MOI.ConstraintName(), index, constraint["name"])
        end
    end
end

function read_objective_sense(sense::String)
    if sense == "min"
        return MOI.MIN_SENSE
    elseif sense == "max"
        return MOI.MAX_SENSE
    else
        @assert sense == "feasibility"
        return MOI.FEASIBILITY_SENSE
    end
end

"""
    function_to_moi(x::OrderedDict, model::Model)

Convert `x` from an OrderedDict representation into a MOI representation.
"""
function function_to_moi(x::Object, args...)
    val_type = Val{Symbol(x["head"])}()
    return function_to_moi(val_type, x, args...)
end

function function_to_moi(::Val{FunctionSymbol}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex}) where FunctionSymbol
    error("Version $(VERSION) of MathOptFormat does not support the function:" *
          " $(FunctionSymbol).")
end

# ========== Non-typed scalar functions ==========

function function_to_moi(::Val{:SingleVariable}, object::Object, model::Model,
                         name_map::Dict{String, MOI.VariableIndex})
    return MOI.SingleVariable(name_map[object["variable"]])
end

# ========== Typed scalar functions ==========

# Here, we deal with a special case: ScalarAffineTerm, ScalarQuadraticTerm,
# VectorAffineTerm, and VectorQuadraticTerm do not contain a "head" field
# (because it is unnecessary at the JSON level).

function parse_scalar_affine_term(
        object::Object, model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.ScalarAffineTerm(
        object["coefficient"],
        name_map[object["variable"]]
    )
end

function parse_scalar_affine_terms(
        terms, model::Model, name_map::Dict{String, MOI.VariableIndex})
    out_terms = MOI.ScalarAffineTerm{Float64}[]
    for term in terms
        push!(out_terms, parse_scalar_affine_term(term, model, name_map))
    end
    return out_terms
end

function function_to_moi(::Val{:ScalarAffineFunction}, object::Object,
                         model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.ScalarAffineFunction(
        parse_scalar_affine_terms(object["terms"], model, name_map),
        object["constant"]
    )
end

function parse_scalar_quadratic_term(
        object::Object, model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.ScalarQuadraticTerm(
        object["coefficient"],
        name_map[object["variable_1"]],
        name_map[object["variable_2"]]
    )
end

function parse_scalar_quadratic_terms(
        terms, model::Model, name_map::Dict{String, MOI.VariableIndex})
    out_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for term in terms
        push!(out_terms, parse_scalar_quadratic_term(term, model, name_map))
    end
    return out_terms
end

function function_to_moi(::Val{:ScalarQuadraticFunction}, object::Object,
                         model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.ScalarQuadraticFunction(
        parse_scalar_affine_terms(object["affine_terms"], model, name_map),
        parse_scalar_quadratic_terms(object["quadratic_terms"], model, name_map),
        object["constant"]
    )
end

# ========== Non-typed vector functions ==========

function function_to_moi(::Val{:VectorOfVariables}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex})
    return MOI.VectorOfVariables(
        [name_map[variable] for variable in object["variables"]]
    )
end

# ========== Typed vector functions ==========

function parse_vector_affine_term(
        object::Object, model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.VectorAffineTerm(
        object["output_index"],
        parse_scalar_affine_term(object["scalar_term"], model, name_map)
    )
end

function parse_vector_affine_terms(
        terms, model::Model, name_map::Dict{String, MOI.VariableIndex})
    out_terms = MOI.VectorAffineTerm{Float64}[]
    for term in terms
        push!(out_terms, parse_vector_affine_term(term, model, name_map))
    end
    return out_terms
end

function function_to_moi(::Val{:VectorAffineFunction}, object::Object, model::Model,
                         name_map::Dict{String, MOI.VariableIndex})
    return MOI.VectorAffineFunction(
        parse_vector_affine_terms(object["terms"], model, name_map),
        Float64.(object["constants"])
    )
end

function parse_vector_quadratic_term(
        object::Object, model::Model, name_map::Dict{String, MOI.VariableIndex})
    return MOI.VectorQuadraticTerm(
        object["output_index"],
        parse_scalar_quadratic_term(object["scalar_term"], model, name_map)
    )
end

function parse_vector_quadratic_terms(
        terms, model::Model, name_map::Dict{String, MOI.VariableIndex})
    out_terms = MOI.VectorQuadraticTerm{Float64}[]
    for term in terms
        push!(out_terms, parse_vector_quadratic_term(term, model, name_map))
    end
    return out_terms
end

function function_to_moi(::Val{:VectorQuadraticFunction}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex})
    return MOI.VectorQuadraticFunction(
        parse_vector_affine_terms(object["affine_terms"], model, name_map),
        parse_vector_quadratic_terms(object["quadratic_terms"], model, name_map),
        Float64.(object["constants"])
    )
end

# ========== Default fallback ==========
"""
    set_to_moi(x::OrderedDict, model::Model, name_map::Dict{String, MOI.VariableIndex})

Convert `x` from an OrderedDict representation into a MOI representation.
"""
function set_to_moi(x::Object, args...)
    val_type = Val{Symbol(x["head"])}()
    return set_to_moi(val_type, x, args...)
end

"""
    set_info(::Type{Val{HeadName}}) where HeadName

Return a tuple of the corresponding MOI set and an ordered list of fieldnames.

`HeadName` is a symbol of the string returned by `head_name(set)`.

    HeadName = Symbol(head_name(set))
    typeof(set_info(Val{HeadName})[1]) == typeof(set)
"""
function set_info(::Type{Val{SetSymbol}}) where SetSymbol
    error("Version $(VERSION) of MathOptFormat does not support the set: " *
          "$(SetSymbol).")
end

function set_to_moi(::Val{SetSymbol}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex}) where SetSymbol
    args = set_info(Val{SetSymbol})
    SetType = args[1]
    if length(args) > 1
        return SetType([object[key] for key in args[2:end]]...)
    else
        return SetType()
    end
end

function set_to_moi(::Val{:SOS1}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex})
    return MOI.SOS1(Float64.(object["weights"]))
end
function set_to_moi(::Val{:SOS2}, object::Object, model::Model,
                       name_map::Dict{String, MOI.VariableIndex})
    return MOI.SOS2(Float64.(object["weights"]))
end

function set_to_moi(
    ::Val{:IndicatorSet}, object::Object, model::Model,
    name_map::Dict{String, MOI.VariableIndex}
)
    set = set_to_moi(object["set"], model, name_map)
    indicator = if object["activate_on"] == "one"
        MOI.ACTIVATE_ON_ONE
    else
        @assert object["activate_on"] == "zero"
        MOI.ACTIVATE_ON_ZERO
    end
    return MOI.IndicatorSet{indicator}(set)
end

# ========== Non-typed scalar sets ==========
set_info(::Type{Val{:ZeroOne}}) = (MOI.ZeroOne,)
set_info(::Type{Val{:Integer}}) = (MOI.Integer,)

# ========== Typed scalar sets ==========
set_info(::Type{Val{:LessThan}}) = (MOI.LessThan, "upper")
set_info(::Type{Val{:GreaterThan}}) = (MOI.GreaterThan, "lower")
set_info(::Type{Val{:EqualTo}}) = (MOI.EqualTo, "value")
set_info(::Type{Val{:Interval}}) = (MOI.Interval, "lower", "upper")
set_info(::Type{Val{:Semiinteger}}) = (MOI.Semiinteger, "lower", "upper")
set_info(::Type{Val{:Semicontinuous}}) = (MOI.Semicontinuous, "lower", "upper")

# ========== Non-typed vector sets ==========
set_info(::Type{Val{:Zeros}}) = (MOI.Zeros, "dimension")
set_info(::Type{Val{:Reals}}) = (MOI.Reals, "dimension")
set_info(::Type{Val{:Nonnegatives}}) = (MOI.Nonnegatives, "dimension")
set_info(::Type{Val{:Nonpositives}}) = (MOI.Nonpositives, "dimension")
set_info(::Type{Val{:SecondOrderCone}}) = (MOI.SecondOrderCone, "dimension")
function set_info(::Type{Val{:RotatedSecondOrderCone}})
    return MOI.RotatedSecondOrderCone, "dimension"
end
set_info(::Type{Val{:GeometricMeanCone}}) = (MOI.GeometricMeanCone, "dimension")
set_info(::Type{Val{:NormOneCone}}) = (MOI.NormOneCone, "dimension")
set_info(::Type{Val{:NormInfinityCone}}) = (MOI.NormInfinityCone, "dimension")
set_info(::Type{Val{:RelativeEntropyCone}}) = (MOI.RelativeEntropyCone, "dimension")
set_info(::Type{Val{:NormSpectralCone}}) = (MOI.NormSpectralCone, "row_dim", "column_dim")
set_info(::Type{Val{:NormNuclearCone}}) = (MOI.NormNuclearCone, "row_dim", "column_dim")
function set_info(::Type{Val{:RootDetConeTriangle}})
    return MOI.RootDetConeTriangle, "side_dimension"
end
function set_info(::Type{Val{:RootDetConeSquare}})
    return MOI.RootDetConeSquare, "side_dimension"
end
function set_info(::Type{Val{:LogDetConeTriangle}})
    return MOI.LogDetConeTriangle, "side_dimension"
end
function set_info(::Type{Val{:LogDetConeSquare}})
    return MOI.LogDetConeSquare, "side_dimension"
end
function set_info(::Type{Val{:PositiveSemidefiniteConeTriangle}})
    return MOI.PositiveSemidefiniteConeTriangle, "side_dimension"
end
function set_info(::Type{Val{:PositiveSemidefiniteConeSquare}})
    return MOI.PositiveSemidefiniteConeSquare, "side_dimension"
end
set_info(::Type{Val{:ExponentialCone}}) = (MOI.ExponentialCone, )
set_info(::Type{Val{:DualExponentialCone}}) = (MOI.DualExponentialCone, )

# ========== Typed vector sets ==========
set_info(::Type{Val{:PowerCone}}) = (MOI.PowerCone, "exponent")
set_info(::Type{Val{:DualPowerCone}}) = (MOI.DualPowerCone, "exponent")


text:
# # Rossby Parameter
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/example_beta.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/example_beta.ipynb)
#
# [from wikipedia:](https://en.wikipedia.org/wiki/Rossby_parameter)
#
# The Rossby parameter `β` arises due to the variation of the Coriolis
# force `f(φ)` as a function of latitude `φ`. It is important
# in the generation of Rossby waves and given by
#
# $$\beta(\varphi) = 2\Omega / R * cos(\varphi)$$

using PhysicalOceanography
isapprox(β(45),1.61687441686e-11)


text:
module QuantumOperators

using Printf: @printf, @sprintf
using ...Prerequisites: atol, rtol, decimaltostr
using ...Prerequisites.Traits: efficientoperations, getcontent, contentorder
using ...Prerequisites.Traits: rawtype, fulltype, parametertype, parameterpairs, reparameter, promoteparameters

import LaTeXStrings: latexstring
import LinearAlgebra: ishermitian
import ..Essentials: dtype
import ...Interfaces: id, value, rank, add!, sub!, mul!, div!, ⊗, ⋅, permute
import ...Prerequisites.Traits: contentnames, dissolve, isparameterbound, parameternames

export QuantumOperator, OperatorUnit, ID, OperatorPack, OperatorProd, OperatorSum, Operator, Operators, LaTeX
export Transformation, Identity, Numericalization, MatrixRepresentation, Permutation, AbstractSubstitution, AbstractUnitSubstitution, UnitSubstitution, RankFilter
export ishermitian, idtype, optype, sequence, latexname, latexformat, superscript, subscript, script, matrix, matrix!

# Generic quantum operator
"""
    QuantumOperator

The abstract type of any quantum operator.
"""
abstract type QuantumOperator end
@inline Base.:(==)(m₁::QuantumOperator, m₂::QuantumOperator) = ==(efficientoperations, m₁, m₂)
@inline Base.isequal(m₁::QuantumOperator, m₂::QuantumOperator) = isequal(efficientoperations, m₁, m₂)

"""
    replace(m::QuantumOperator; kwargs...) -> typeof(m)

Return a copy of a concrete `QuantumOperator` with some of the field values replaced by the keyword arguments.
"""
@inline Base.replace(m::QuantumOperator; kwargs...) = replace(efficientoperations, m; kwargs...)

# Operator unit
"""
    OperatorUnit <: QuantumOperator

An operator unit is the irreducible symbolic unit to represent a quantum operator.

It plays the role of the symbols as in usual computer algebras while it can host internal structures, which is convenient for quantum operators in representative of the internal degrees of freedom.
"""
abstract type OperatorUnit <: QuantumOperator end
@inline Base.show(io::IO, u::OperatorUnit) = @printf io "%s(%s)" nameof(typeof(u)) join(map(repr, ntuple(i->getfield(u, i), Val(fieldcount(typeof(u))))), ", ")

# ID of a composite quantum operator
"""
    ID{U<:OperatorUnit, N}

The id of a composite quantum operator, which is an ordered set of operator units.

Type alias for `NTuple{N, U} where {U<:OperatorUnit}`.
"""
const ID{U<:OperatorUnit, N} = NTuple{N, U}
@inline Base.show(io::IO, id::Tuple{OperatorUnit, Vararg{OperatorUnit}}) = @printf io "ID(%s)" join(id, ", ")
@inline Base.promote_rule(::Type{Tuple{}}, I::Type{<:Tuple{OperatorUnit, Vararg{OperatorUnit}}}) = ID{I|>eltype}
@inline Base.promote_rule(I::Type{<:Tuple{OperatorUnit, Vararg{OperatorUnit}}}, ::Type{Tuple{}}) = ID{I|>eltype}

"""
    ID(id::OperatorUnit...)
    ID(u::OperatorUnit, id::ID{OperatorUnit})
    ID(id::ID{OperatorUnit}, u::OperatorUnit)
    ID(id₁::ID{OperatorUnit}, id₂::ID{OperatorUnit})

Get the id from operator units/ids.
"""
@inline ID(id::OperatorUnit...) = id
@inline ID(u::OperatorUnit, id::ID{OperatorUnit}) = ID(u, id...)
@inline ID(id::ID{OperatorUnit}, u::OperatorUnit) = ID(id..., u)
@inline ID(id₁::ID{OperatorUnit}, id₂::ID{OperatorUnit}) = ID(id₁..., id₂...)

"""
    ID(::Type{U}, attrs::Vararg{NTuple{N}, M}) where {U<:OperatorUnit, N, M}

Get the composite id from the components of singular ids.
"""
@inline @generated function ID(::Type{U}, attrs::Vararg{NTuple{N, Any}, M}) where {U<:OperatorUnit, N, M}
    exprs = []
    for i = 1:N
        args = [:(attrs[$j][$i]) for j = 1:M]
        push!(exprs, :(U($(args...))))
    end
    return :(ID($(exprs...)))
end

"""
    propertynames(::Type{I}) where I<:ID{OperatorUnit} -> Tuple{Vararg{Symbol}}

Get the property names of a composite id.
"""
@inline @generated function Base.propertynames(I::ID{OperatorUnit})
    exprs = [QuoteNode(Symbol(name, 's')) for name in I|>eltype|>fieldnames]
    return Expr(:tuple, exprs...)
end

"""
    getproperty(id::ID{OperatorUnit}, name::Symbol)

Get the property of a composite id.
"""
@inline Base.getproperty(id::ID{OperatorUnit}, name::Symbol) = idgetproperty(id, Val(name), Val(id|>propertynames))
@inline @generated function idgetproperty(id::ID{OperatorUnit}, ::Val{name}, ::Val{names}) where {name, names}
    index = findfirst(isequal(name), names)::Int
    exprs = [:(getfield(id[$i], $index)) for i = 1:fieldcount(id)]
    return Expr(:tuple, exprs...)
end

"""
    rank(id::ID{OperatorUnit}) -> Int
    rank(::Type{<:ID{OperatorUnit}}) -> Any
    rank(::Type{<:ID{OperatorUnit, N}}) where N -> Int

Get the rank of an id.
"""
@inline rank(id::ID{OperatorUnit}) = rank(typeof(id))
@inline rank(::Type{<:ID{OperatorUnit}}) = Any
@inline rank(::Type{<:ID{OperatorUnit, N}}) where N = N

"""
    adjoint(id::ID{OperatorUnit}) -> ID

Get the adjoint of an id.
"""
@inline Base.adjoint(id::ID{OperatorUnit}) = map(adjoint, reverse(id))

"""
    ishermitian(id::ID{OperatorUnit}) -> Bool

Judge whether an id is Hermitian.
"""
function ishermitian(id::ID{OperatorUnit})
    for i = 1:((rank(id)+1)÷2)
        id[i]'==id[rank(id)+1-i] || return false
    end
    return true
end

# Operator pack
"""
    OperatorPack{V, I<:ID{OperatorUnit}} <: QuantumOperator

The entity that represent the pack of a number and several quantum units.

Basically, a concrete subtype should contain two predefined contents:
- `value::V`: the coefficient of the pack
- `id::I`: the total id of the pack
"""
abstract type OperatorPack{V, I<:ID{OperatorUnit}} <: QuantumOperator end
@inline contentnames(::Type{<:OperatorPack}) = (:value, :id)
@inline parameternames(::Type{<:OperatorPack}) = (:value, :id)
@inline isparameterbound(::Type{<:OperatorPack}, ::Val{:value}, ::Type{V}) where V = false
@inline isparameterbound(::Type{<:OperatorPack}, ::Val{:id}, ::Type{I}) where {I<:ID{OperatorUnit}} = !isconcretetype(I)
@inline function Base.promote_rule(::Type{M₁}, ::Type{M₂}) where {M₁<:OperatorPack, M₂<:OperatorPack}
    M₁<:M₂ && return M₂
    M₂<:M₁ && return M₁
    r₁, r₂ = rank(M₁), rank(M₂)
    M = r₂==0 ? rawtype(M₁) : r₁==0 ? rawtype(M₂) : typejoin(rawtype(M₁), rawtype(M₂))
    return fulltype(M, promoteparameters(parameterpairs(M₁), parameterpairs(M₂)))
end
@inline Base.promote_rule(::Type{M}, ::Type{N}) where {M<:OperatorPack, N<:Number} = reparameter(M, :value, promote_type(valtype(M), N))

"""
    rank(m::OperatorPack) -> Int
    rank(::Type{M}) where {M<:OperatorPack} -> Int

Get the rank of an `OperatorPack`.
"""
@inline rank(m::OperatorPack) = rank(typeof(m))
@inline rank(::Type{M}) where {M<:OperatorPack} = rank(idtype(M))

"""
    valtype(m::OperatorPack)
    valtype(::Type{T}) where {T<:OperatorPack}

Get the type of the value of an `OperatorPack`.
"""
@inline Base.valtype(m::OperatorPack) = valtype(typeof(m))
@inline @generated Base.valtype(::Type{T}) where {T<:OperatorPack} = parametertype(supertype(T, :OperatorPack), :value)

"""
    idtype(m::OperatorPack)
    idtype(::Type{T}) where {T<:OperatorPack}

The type of the id of an `OperatorPack`.
"""
@inline idtype(m::OperatorPack) = idtype(typeof(m))
@inline @generated idtype(::Type{T}) where {T<:OperatorPack} = parametertype(supertype(T, :OperatorPack), :id)

"""
    dtype(m::OperatorPack)
    dtype(::Type{T}) where {T<:OperatorPack}

The data type of the coefficient of an `OperatorPack`.
"""
@inline dtype(m::OperatorPack) = dtype(typeof(m))
@inline dtype(::Type{T}) where {T<:OperatorPack} = valtype(T)

"""
    value(m::OperatorPack) -> valtype(m)

Get the value of an `OperatorPack`.
"""
@inline value(m::OperatorPack) = getcontent(m, :value)

"""
    id(m::OperatorPack) -> idtype(m)

Get the id of an `OperatorPack`.
"""
@inline id(m::OperatorPack) = getcontent(m, :id)

@inline newvalue(m::OperatorPack, v) = v
"""
    replace(m::OperatorPack, v) -> OperatorPack

Replace the value of an `OperatorPack`.
"""
@inline Base.replace(m::OperatorPack, v) = rawtype(typeof(m))(dissolve(m, newvalue, (v,))...)
@inline dissolve(m::OperatorPack, ::Val{:value}, ::typeof(newvalue), args::Tuple, kwargs::NamedTuple) = newvalue(m, args...; kwargs...)

"""
    isapprox(m₁::OperatorPack, m₂::OperatorPack; atol::Real=atol, rtol::Real=rtol) -> Bool

Compare two `OperatorPack`s and judge whether they are approximate to each other.
"""
@inline function Base.isapprox(m₁::OperatorPack, m₂::OperatorPack; atol::Real=atol, rtol::Real=rtol)
    isapprox(efficientoperations, contentorder(typeof(m₁), :value)|>Val, dissolve(m₁), dissolve(m₂); atol=atol, rtol=rtol)::Bool
end

"""
    one(::Type{M}) where M<:OperatorPack
    one(m::OperatorPack)

Get the identity quantum operator.
"""
@inline function Base.one(::Type{M}) where M<:OperatorPack
    rtype = rawtype(M)
    vtype = isconcretetype(valtype(M)) ? valtype(M) : Int
    @assert fieldnames(rtype) == (:value, :id) "one error: not supported type($(nameof(rtype)))."
    return rtype(one(vtype), ())
end
@inline Base.one(m::OperatorPack) = rawtype(typeof(m))(dissolve(m, one)...)
@inline dissolve(m::OperatorPack, ::Val{:value}, ::typeof(one), ::Tuple, ::NamedTuple) = one(valtype(m))
@inline dissolve(::OperatorPack, ::Val{:id}, ::typeof(one), ::Tuple, ::NamedTuple) = ID()

"""
    convert(::Type{M}, m::Number) where {M<:OperatorPack}
    convert(::Type{M}, m::OperatorPack) where {M<:OperatorPack}

1) Convert a number to a quantum operator.
2) Convert a quantum operator from one type to another.
"""
@inline function Base.convert(::Type{M}, m::Number) where {M<:OperatorPack}
    @assert isa(one(M), M) "convert error: not convertible."
    return one(M)*convert(dtype(M), m)
end
@inline function Base.convert(::Type{M}, m::OperatorPack) where {M<:OperatorPack}
    (typeof(m) <: M) && return m
    @assert convertible(M, m) "convert error: $(typeof(m)) cannot be converted to $M."
    return replace(m, convert(valtype(M), value(m)))
end
function convertible(::Type{M}, m::OperatorPack) where {M<:OperatorPack}
    !(rawtype(typeof(m)) <: rawtype(M)) && return false
    S = supertype(typeof(m), nameof(M))
    for (i, name) in enumerate(parameternames(M))
        (name != :value) && !(parametertype(S, i) <: parametertype(M, i)) && return false
    end
    return true
end

# Operator prod
"""
    OperatorProd{V, I<:ID{OperatorUnit}} <: OperatorPack{V, I}

A special kind of `OperatorPack`, where the relation between the coefficient and quantum units could be viewed as product.
"""
abstract type OperatorProd{V, I<:ID{OperatorUnit}} <: OperatorPack{V, I} end
@inline Base.eltype(m::OperatorProd) = eltype(typeof(m))
@inline Base.eltype(::Type{M}) where {M<:OperatorProd} = eltype(idtype(M))
@inline Base.iterate(m::OperatorProd) = iterate(id(m))
@inline Base.iterate(m::OperatorProd, state) = iterate(id(m), state)

"""
    length(m::OperatorProd) -> Int

Get the length of an `OperatorProd`.
"""
@inline Base.length(m::OperatorProd) = rank(m)
@inline Base.firstindex(m::OperatorProd) = 1
@inline Base.lastindex(m::OperatorProd) = rank(m)

"""
    getindex(m::OperatorProd, i::Integer) -> OperatorUnit
    getindex(m::OperatorProd, slice) -> OperatorProd

Overloaded `[]`.
"""
@inline Base.getindex(m::OperatorProd, i::Integer) = id(m)[i]
@inline Base.getindex(m::OperatorProd, slice) = rawtype(typeof(m))(dissolve(m, getindex, (slice,))...)
@inline dissolve(m::OperatorProd, ::Val{:value}, ::typeof(getindex), ::Tuple{Any}, ::NamedTuple) = one(valtype(m))
@inline dissolve(m::OperatorProd, ::Val{:id}, ::typeof(getindex), slice::Tuple{Any}, ::NamedTuple) = id(m)[first(slice)]

"""
    split(m::OperatorProd) -> Tuple{valtype(m), Vararg{OperatorUnit}}

Split an `OperatorProd` into the coefficient and a sequence of `OperatorUnit`s.
"""
@inline Base.split(m::OperatorProd) = (value(m), id(m)...)

"""
    sequence(m::OperatorProd, table) -> NTuple{rank(m), Int}

Get the sequence of the id of a quantum operator according to a table.
"""
@inline sequence(m::OperatorProd, table) = map(u->table[u], id(m))

# Operator
"""
    Operator{V<:Number, I<:ID{OperatorUnit}} <: OperatorProd{V, I}

Operator.
"""
struct Operator{V<:Number, I<:ID{OperatorUnit}} <: OperatorProd{V, I}
    value::V
    id::I
end
@inline Operator(value::Number, id::OperatorUnit...) = Operator(value, id)
function Base.show(io::IO, m::Operator)
    @printf io "%s(%s%s%s)" nameof(typeof(m)) decimaltostr(value(m)) (rank(m)>0 ? ", " : "") join(id(m), ", ")
end

"""
    adjoint(m::Operator) -> Operator

Get the adjoint of an operator.
"""
@inline Base.adjoint(m::Operator) = rawtype(typeof(m))(value(m)', id(m)')

"""
    ishermitian(m::Operator) -> Bool

Judge whether an operator is Hermitian.
"""
@inline ishermitian(m::Operator) = isa(value(m), Real) && ishermitian(id(m))

"""
    convert(::Type{M}, u::OperatorUnit) where {M<:Operator{<:Number, <:ID{OperatorUnit}}}

Convert an operator unit to an operator.
"""
@inline function Base.convert(::Type{M}, u::OperatorUnit) where {M<:Operator{<:Number, <:ID{OperatorUnit}}}
    @assert Tuple{typeof(u)} <: idtype(M) "convert error: not convertible."
    return Operator(one(valtype(M)), ID(u))
end

# Operator sum
"""
    OperatorSum{M<:OperatorPack, I<:ID{OperatorUnit}} <: QuantumOperator

The sum of `OperatorPack`s.

Similar items are automatically merged with the aid of the id system.
"""
struct OperatorSum{M<:OperatorPack, I<:ID{OperatorUnit}} <: QuantumOperator
    contents::Dict{I, M}
    OperatorSum(contents::Dict{<:ID{OperatorUnit}, <:OperatorPack}) = new{valtype(contents), keytype(contents)}(contents)
end
@inline Base.eltype(ms::OperatorSum) = eltype(typeof(ms))
@inline Base.eltype(::Type{<:OperatorSum{M}}) where {M<:OperatorPack} = M
@inline Base.iterate(ms::OperatorSum) = iterate(values(ms.contents))
@inline Base.iterate(ms::OperatorSum, state) = iterate(values(ms.contents), state)
@inline Base.length(ms::OperatorSum) = length(ms.contents)
function Base.show(io::IO, ms::OperatorSum)
    @printf io "%s with %s %s\n" summary(ms) length(ms) nameof(eltype(ms))
    for m in ms
        @printf io "  %s\n" m
    end
end
@inline Base.haskey(ms::OperatorSum, id::ID{OperatorUnit}) = haskey(ms.contents, id)
@inline Base.getindex(ms::OperatorSum, id::ID{OperatorUnit}) = ms.contents[id]
@inline Base.setindex!(ms::OperatorSum, m::OperatorPack, id::ID{OperatorUnit}) = (ms.contents[id] = m; m)
@inline Base.empty(ms::OperatorSum) = OperatorSum(empty(ms.contents))
@inline Base.empty!(ms::OperatorSum) = (empty!(ms.contents); ms)
@inline Base.merge(ms::OperatorSum, another::OperatorSum) = OperatorSum(merge(ms.contents, another.contents))
@inline Base.merge!(ms::OperatorSum, another::OperatorSum) = (merge!(ms.contents, another.contents); ms)
@inline function Base.promote_rule(::Type{MS₁}, ::Type{MS₂}) where {MS₁<:OperatorSum, MS₂<:OperatorSum}
    M = promote_type(eltype(MS₁), eltype(MS₂))
    return OperatorSum{M, idtype(M)}
end
@inline function Base.convert(::Type{MS}, ms::OperatorSum) where {MS<:OperatorSum}
    @assert eltype(MS)>:eltype(ms) "convert error: cannot convert an object of $(typeof(ms)) to an object of $(MS)."
    return add!(zero(MS), ms)
end

"""
    OperatorSum(ms)
    OperatorSum(ms::QuantumOperator...)
    OperatorSum{M}(ms) where {M<:OperatorPack}
    OperatorSum{M}(ms::QuantumOperator...) where {M<:OperatorPack}

Get the sum of `OperatorPack`s.
"""
@inline OperatorSum(ms::QuantumOperator...) = OperatorSum{eltype(ms)}(ms)
@inline OperatorSum(ms) = OperatorSum{eltype(ms)}(ms)
@inline OperatorSum{M}(ms::QuantumOperator...) where {M<:OperatorPack} = OperatorSum{M}(ms)
function OperatorSum{M}(ms) where {M<:OperatorPack}
    result = OperatorSum(Dict{idtype(M), M}())
    for m in ms
        add!(result, m)
    end
    return result
end

"""
    isapprox(ms₁::OperatorSum, ms₂::OperatorSum; atol::Real=atol, rtol::Real=rtol) -> Bool

Compare two `OperatorSum`s and judge whether they are approximate to each other.
"""
@inline function Base.isapprox(ms₁::OperatorSum, ms₂::OperatorSum; atol::Real=atol, rtol::Real=rtol)
    for m in ms₁
        isapprox(value(m), 0, atol=atol, rtol=rtol) && continue
        k = id(m)
        haskey(ms₂, k) || return false
        isapprox(m, ms₂[k]) || return false
    end
    for m in ms₂
        isapprox(value(m), 0,  atol=atol, rtol=rtol) && continue
        k = id(m)
        haskey(ms₁, k) || return false
        isapprox(m, ms₁[k]) || return false
    end
    return true
end

"""
    zero(ms::OperatorSum) -> OperatorSum
    zero(::Type{MS}) where {MS<:OperatorSum} -> OperatorSum

Get the zero sum.
"""
@inline Base.zero(ms::OperatorSum) = zero(typeof(ms))
@inline Base.zero(::Type{MS}) where {MS<:OperatorSum} = OperatorSum{eltype(MS)}()

"""
    add!(ms::OperatorSum) -> typeof(ms)
    add!(ms::OperatorSum, m::Union{Number, OperatorUnit, OperatorPack}) -> typeof(ms)
    add!(ms::OperatorSum, mms::OperatorSum) -> typeof(ms)

Get the in-place addition of quantum operators.
"""
@inline add!(ms::OperatorSum) = ms
@inline function add!(ms::OperatorSum, m::Union{Number, OperatorUnit, OperatorPack})
    isa(m, Number) && m==0 && return ms
    m = convert(eltype(ms), m)
    old = get(ms.contents, id(m), nothing)
    new = isnothing(old) ? m : replace(old, value(old)+value(m))
    value(new)==0 ? delete!(ms.contents, id(m)) : (ms[id(m)] = new)
    return ms
end
@inline function add!(ms::OperatorSum, mms::OperatorSum)
    for m in mms
        add!(ms, m)
    end
    return ms
end

"""
    sub!(ms::OperatorSum) -> typeof(ms)
    sub!(ms::OperatorSum, m::Union{Number, OperatorUnit, OperatorPack}) -> typeof(ms)
    sub!(ms::OperatorSum, mms::OperatorSum) -> typeof(ms)

Get the in-place subtraction of quantum operators.
"""
@inline sub!(ms::OperatorSum) = ms
@inline function sub!(ms::OperatorSum, m::Union{Number, OperatorUnit, OperatorPack})
    isa(m, Number) && abs(m)==0 && return ms
    m = convert(eltype(ms), m)
    old = get(ms.contents, id(m), nothing)
    new = isnothing(old) ? -m : replace(old, value(old)-value(m))
    value(new)==0 ? delete!(ms.contents, id(m)) : (ms[id(m)] = new)
    return ms
end
@inline function sub!(ms::OperatorSum, mms::OperatorSum)
    for m in mms
        sub!(ms, m)
    end
    return ms
end

"""
    mul!(ms::OperatorSum, factor::Number) -> OperatorSum

Get the in-place multiplication of an `OperatorSum` with a number.
"""
function mul!(ms::OperatorSum, factor::Number)
    abs(factor)==0 && return empty!(ms)
    for m in ms
        new = replace(m, value(m)*factor)
        value(new)==0 ? delete!(ms.contents, id(m)) : (ms[id(m)] = new)
    end
    return ms
end

"""
    div!(ms::OperatorSum, factor::Number) -> OperatorSum

Get the in-place division of an `OperatorSum` with a number.
"""
@inline div!(ms::OperatorSum, factor::Number) = mul!(ms, one(dtype(eltype(ms)))/factor)

"""
    optype(m::QuantumOperator)
    optype(::Type{<:QuantumOperator})

Get the corresponding `OperatorPack` type of a generic quantum operator.
"""
@inline optype(m::QuantumOperator) = optype(typeof(m))
@inline optype(::Type{M}) where {M<:OperatorUnit} = fulltype(Operator, NamedTuple{(:value, :id), Tuple{Int, Tuple{M}}})
@inline optype(::Type{M}) where {M<:OperatorPack} = M
@inline optype(::Type{M}) where {M<:OperatorSum} = eltype(M)

"""
    +(m::QuantumOperator) -> typeof(m)
    +(m₁::QuantumOperator, m₂::QuantumOperator) -> OperatorSum
    +(factor::Number, m::QuantumOperator) -> OperatorSum
    +(m::QuantumOperator, factor::Number) -> OperatorSum

Overloaded `+` between quantum operators.
"""
@inline Base.:+(m::QuantumOperator) = m
@inline function Base.:+(m₁::QuantumOperator, m₂::QuantumOperator)
    M = promote_type(optype(m₁), optype(m₂))
    result = OperatorSum{M}()
    add!(result, m₁)
    add!(result, m₂)
    return result
end
@inline Base.:+(factor::Number, m::QuantumOperator) = one(optype(m))*factor + m
@inline Base.:+(m::QuantumOperator, factor::Number) = m + one(optype(m))*factor

"""
    -(m::QuantumOperator) -> QuantumOperator
    -(m₁::QuantumOperator, m₂::QuantumOperator) -> OperatorSum
    -(factor::Number, m::QuantumOperator) -> OperatorSum
    -(m::QuantumOperator, factor::Number) -> OperatorSum

Overloaded `-` between quantum operators.
"""
@inline Base.:-(m::QuantumOperator) = m*(-1)
@inline function Base.:-(m₁::QuantumOperator, m₂::QuantumOperator)
    M = promote_type(optype(m₁), optype(m₂))
    result = OperatorSum{M}()
    add!(result, m₁)
    sub!(result, m₂)
    return result
end
@inline Base.:-(factor::Number, m::QuantumOperator) = one(optype(m))*factor - m
@inline Base.:-(m::QuantumOperator, factor::Number) = m - one(optype(m))*factor

"""
    *(factor::Number, m::OperatorUnit) -> Operator
    *(m::OperatorUnit, factor::Number) -> Operator
    *(m₁::OperatorUnit, m₂::OperatorUnit) -> Operator
    *(factor::Number, m::OperatorPack) -> OperatorPack
    *(m::OperatorPack, factor::Number) -> OperatorPack
    *(m₁::OperatorPack, m₂::OperatorPack) -> OperatorPack
    *(factor::Number, ms::OperatorSum) -> OperatorSum
    *(ms::OperatorSum, factor::Number) -> OperatorSum
    *(m::OperatorPack, ms::OperatorSum) -> OperatorSum
    *(ms::OperatorSum, m::OperatorPack) -> OperatorSum
    *(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `*` between quantum operators or a quantum operator and a number.
"""
@inline Base.:*(factor::Number, m::OperatorUnit) = Operator(factor, ID(m))
@inline Base.:*(m::OperatorUnit, factor::Number) = Operator(factor, ID(m))
@inline Base.:*(m₁::OperatorUnit, m₂::OperatorUnit) = Operator(1, ID(m₁, m₂))
@inline Base.:*(factor::Number, m::OperatorPack) = replace(m, factor*value(m))
@inline Base.:*(m::OperatorPack, factor::Number) = replace(m, value(m)*factor)
@inline function Base.:*(m₁::OperatorPack, m₂::OperatorPack)
    M₁, M₂ = typeof(m₁), typeof(m₂)
    @assert nameof(M₁)==nameof(M₂) && contentnames(M₁)==(:value, :id)==contentnames(M₂) "\"*\" error: not implemented between $(nameof(M₁)) and $(nameof(M₂))."
    return rawtype(M₁)(value(m₁)*value(m₂), ID(id(m₁), id(m₂)))
end
@inline Base.:*(factor::Number, ms::OperatorSum) = ms * factor
function Base.:*(ms::OperatorSum, factor::Number)
    abs(factor)==0 && return zero(ms)
    result = OperatorSum{promote_type(eltype(ms), typeof(factor))}()
    for m in ms
        add!(result, m*factor)
    end
    return result
end
@inline Base.:*(m::OperatorPack, ms::OperatorSum) = OperatorSum(collect(m*mm for mm in ms))
@inline Base.:*(ms::OperatorSum, m::OperatorPack) = OperatorSum(collect(mm*m for mm in ms))
@inline Base.:*(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum(collect(m₁*m₂ for m₁ in ms₁ for m₂ in ms₂))

"""
    /(m::QuantumOperator, factor::Number) -> QuantumOperator

Overloaded `/` between a quantum operator and a number.
"""
@inline Base.:/(m::QuantumOperator, factor::Number) = m * (one(dtype(optype(m)))/factor)

"""
    //(m::QuantumOperator, factor::Integer) -> QuantumOperator

Overloaded `//` between a quantum operator and an integer.
"""
@inline Base.://(m::QuantumOperator, factor::Integer) = m * (1//factor)

"""
    ^(m::QuantumOperator, n::Integer) -> QuantumOperator

Overloaded `^` between a quantum operator and an integer.
"""
@inline Base.:^(m::QuantumOperator, n::Integer) = (@assert n>0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->m, Val(n)), init=1))

"""
    ⊗(m::Union{OperatorUnit, OperatorPack}, ms::OperatorSum) -> OperatorSum
    ⊗(ms::OperatorSum, m::Union{OperatorUnit, OperatorPack}) -> OperatorSum
    ⊗(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `⊗` between quantum operators.
"""
@inline ⊗(m::Union{OperatorUnit, OperatorPack}, ms::OperatorSum) = OperatorSum(collect(m ⊗ mm for mm in ms))
@inline ⊗(ms::OperatorSum, m::Union{OperatorUnit, OperatorPack}) = OperatorSum(collect(mm ⊗ m for mm in ms))
@inline ⊗(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum(collect(m₁ ⊗ m₂ for m₁ in ms₁ for m₂ in ms₂))

"""
    ⋅(m::Union{OperatorUnit, OperatorPack}, ms::OperatorSum) -> OperatorSum
    ⋅(ms::OperatorSum, m::Union{OperatorUnit, OperatorPack}) -> OperatorSum
    ⋅(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `⋅` between quantum operators.
"""
@inline ⋅(m::Union{OperatorUnit, OperatorPack}, ms::OperatorSum) = OperatorSum(collect(m ⋅ mm for mm in ms))
@inline ⋅(ms::OperatorSum, m::Union{OperatorUnit, OperatorPack}) = OperatorSum(collect(mm ⋅ m for mm in ms))
@inline ⋅(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum(collect(m₁ ⋅ m₂ for m₁ in ms₁ for m₂ in ms₂))

# Operators
"""
    Operators{O<:Operator, I<:ID{OperatorUnit}}

A set of operators.

Type alias for `OperatorSum{O<:Operator, I<:ID{OperatorUnit}}`.
"""
const Operators{O<:Operator, I<:ID{OperatorUnit}} = OperatorSum{O, I}
@inline Base.summary(io::IO, opts::Operators) = @printf io "Operators"

"""
    Operators(opts::Operator...)
    Operators{M}(opts::Operator...)

Get a set of operators.
"""
@inline Operators(opts::Operator...) = OperatorSum(opts)
@inline Operators{M}(opts::Operator...) where {M<:Operator} = OperatorSum{M}(opts)

"""
    adjoint(opts::Operators) -> Operators

Get the adjoint of a set of operators.
"""
function Base.adjoint(opts::Operators)
    result = zero(opts)
    for opt in opts
        add!(result, opt')
    end
    return result
end

"""
    ishermitian(opts::Operators) -> Bool

Judge whether a set of operators as a whole is Hermitian.
"""
@inline ishermitian(opts::Operators) = opts == opts'

# LaTeX format of quantum operators
"""
    LaTeX{SP, SB}(body, spdelimiter::String=", ", sbdelimiter::String=", "; options...) where {SP, SB}

LaTeX string representation of quantum operators.
"""
struct LaTeX{SP, SB, B, O}
    body::B
    spdelimiter::String
    sbdelimiter::String
    options::O
    function LaTeX{SP, SB}(body, spdelimiter::String=", ", sbdelimiter::String=", "; options...) where {SP, SB}
        @assert isa(SP, Tuple{Vararg{Symbol}}) && isa(SB, Tuple{Vararg{Symbol}}) "LaTeX error: SP and SB must be tuple of symbols."
        new{SP, SB, typeof(body), typeof(options)}(body, spdelimiter, sbdelimiter, options)
    end
end
@inline superscript(::Type{<:LaTeX{SP}}) where SP = SP
@inline subscript(::Type{<:LaTeX{SP, SB} where SP}) where SB = SB

"""
    latexname(T::Type{<:OperatorUnit}) -> Symbol

Get the name of a type of `OperatorUnit` in the latex format lookups.
"""
@inline latexname(T::Type{<:OperatorUnit}) = nameof(T)

const latexformats = Dict{Symbol, LaTeX}()
"""
    latexformat(T::Type{<:OperatorUnit}) -> LaTeX
    latexformat(T::Type{<:OperatorUnit}, l::LaTeX) -> LaTeX

Get/Set the LaTeX format for a subtype of `OperatorUnit`.
"""
@inline latexformat(T::Type{<:OperatorUnit}) = latexformats[latexname(T)]
@inline latexformat(T::Type{<:OperatorUnit}, l::LaTeX) = latexformats[latexname(T)] = l

"""
    script(::Val{:BD}, u::OperatorUnit, l::LaTeX) -> Any
    script(::Val{:SP}, u::OperatorUnit, l::LaTeX) -> Tuple
    script(::Val{:SB}, u::OperatorUnit, l::LaTeX) -> Tuple

Get the body/superscript/subscript of the LaTeX string representation of an operator unit.
"""
@inline script(::Val{:BD}, u::OperatorUnit, l::LaTeX) = l.body
@inline @generated script(::Val{:SP}, u::OperatorUnit, l::LaTeX) = Expr(:tuple, [:(script(Val($sup), u; l.options...)) for sup in QuoteNode.(l|>superscript)]...)
@inline @generated script(::Val{:SB}, u::OperatorUnit, l::LaTeX) = Expr(:tuple, [:(script(Val($sub), u; l.options...)) for sub in QuoteNode.(l|>subscript)]...)

"""
    latexstring(u::OperatorUnit) -> String

LaTeX string representation of an operator unit.
"""
@inline function latexstring(u::OperatorUnit)
    l = latexformat(typeof(u))
    return @sprintf "%s^{%s}_{%s}" script(Val(:BD), u, l) join(script(Val(:SP), u, l), l.spdelimiter) join(script(Val(:SB), u, l), l.sbdelimiter)
end

"""
    latexstring(opt::Operator) -> String

Get the string representation of an operator in the LaTeX format.
"""
function latexstring(opt::Operator)
    rank(opt)==0 && return replace(valuetolatextext(value(opt)), " "=>"")
    poses = Int[]
    push!(poses, 1)
    for i = 2:rank(opt)
        id(opt)[i]≠id(opt)[i-1] && push!(poses, i)
    end
    push!(poses, rank(opt)+1)
    result = valuetostr(value(opt))
    for i = 1:length(poses)-1
        order = poses[i+1] - poses[i]
        if order == 1
            result = @sprintf "%s%s" result latexstring(id(opt)[poses[i]])
        else
            result = @sprintf "%s(%s)^%s" result latexstring(id(opt)[poses[i]]) order
        end
    end
    return result
end
function valuetostr(v)
    v==+1 && return ""
    v==-1 && return "-"
    result = valuetolatextext(v)
    if occursin(" ", result) || (isa(v, Complex) && real(v)≠0 && imag(v)≠0)
        bra, ket = occursin("(", result) ? ("[", "]") : ("(", ")")
        result = @sprintf "%s%s%s" bra replace(result, " "=>"") ket
    end
    return result
end
@inline valuetolatextext(value::Union{Real, Complex}) = decimaltostr(value)
function valuetolatextext(value)
    io = IOBuffer()
    if showable(MIME"text/latex"(), value)
        show(IOContext(io, :limit=>false), MIME"text/latex"(), value)
    else
        show(IOContext(io, :limit=>false), MIME"text/plain"(), value)
    end
    return replace(replace(replace(replace(String(take!(io)), "\\begin{equation*}" => ""), "\\end{equation*}" => ""), "\$" => ""), "\n" => "")
end

"""
    latexstring(opts::Operators) -> String

Get the string representation of a set of operators in the LaTeX format.
"""
function latexstring(opts::Operators)
    result = String[]
    for (i, opt) in enumerate(values(opts))
        rep = latexstring(opt)
        i>1 && rep[1]≠'-' && push!(result, "+")
        push!(result, rep)
    end
    return join(result, "")
end

"""
    show(io::IO, ::MIME"text/latex", m::QuantumOperator)

Show a quantum operator.
"""
Base.show(io::IO, ::MIME"text/latex", m::QuantumOperator) = show(io, MIME"text/latex"(), latexstring(latexstring(m)))

# Transformations
"""
    add!(destination, transformation::Function, op::QuantumOperator; kwargs...) -> typeof(destination)

Add the result of the transformation on a quantum operator to the destination.
"""
@inline function add!(destination, transformation::Function, op::QuantumOperator; kwargs...)
    add!(destination, transformation(op; kwargs...))
end

"""
    map!(transformation::Function, ms::OperatorSum; kwargs...) -> typeof(ms)

In place map of an `OperatorSum` by a function.
"""
function Base.map!(transformation::Function, ms::OperatorSum; kwargs...)
    for m in ms
        ms[id(m)] = transformation(m; kwargs...)
    end
    return ms
end

"""
    map!(transformation::Function, destination, ms::OperatorSum; kwargs...) -> typeof(destination)

In place map of an `OperatorSum` by a function.
"""
function Base.map!(transformation::Function, destination, ms::OperatorSum; kwargs...)
    for m in ms
        add!(destination, transformation, m; kwargs...)
    end
    return destination
end

"""
    Transformation <: Function

Abstract transformation that could transform the quantum operators.
"""
abstract type Transformation <: Function end
@inline Base.:(==)(transformation₁::Transformation, transformation₂::Transformation) = ==(efficientoperations, transformation₁, transformation₂)
@inline Base.isequal(transformation₁::Transformation, transformation₂::Transformation) = isequal(efficientoperations, transformation₁, transformation₂)
@inline Base.valtype(transformation::Transformation, m::QuantumOperator) = valtype(typeof(transformation), typeof(m))
@inline Base.zero(transformation::Transformation, m::QuantumOperator) = zero(valtype(transformation, m))

"""
    (transformation::Transformation)(ms::OperatorSum; kwargs...) -> OperatorSum

Get the transformed quantum operators.
"""
function (transformation::Transformation)(ms::OperatorSum; kwargs...)
    result = zero(transformation, ms)
    for m in ms
        add!(result, transformation, m; kwargs...)
    end
    return result
end

"""
    Identity <: Transformation

The identity transformation.
"""
struct Identity <: Transformation end
@inline Base.valtype(::Type{Identity}, M::Type{<:QuantumOperator}) = M
@inline (i::Identity)(m::Union{OperatorUnit, OperatorPack}; kwargs...) = m

"""
    Numericalization{T<:Number} <: Transformation

The numericalization transformation, which converts the value of a quantum operator to a number of type `T`.
"""
struct Numericalization{T<:Number} <: Transformation end
@inline Base.valtype(num::Numericalization) = valtype(typeof(num))
@inline Base.valtype(::Type{<:Numericalization{T}}) where {T<:Number} = T
@inline Base.valtype(T::Type{<:Numericalization}, M::Type{<:OperatorPack}) = reparameter(M, :value, valtype(T))
@inline function Base.valtype(T::Type{<:Numericalization}, M::Type{<:OperatorSum})
    V = valtype(T, eltype(M))
    return OperatorSum{V, idtype(V)}
end
@inline (n::Numericalization)(m::OperatorPack; kwargs...) = convert(valtype(n, m), m)

"""
    MatrixRepresentation <: Transformation

The matrix representation.
"""
abstract type MatrixRepresentation <: Transformation end

"""
    matrix

Generic matrix representation.
"""
function matrix end

"""
    matrix!

Generic matrix representation.
"""
function matrix! end

"""
    Permutation{T} <: Transformation

The permutation transformation.
"""
struct Permutation{T} <: Transformation
    table::T
end
@inline function Base.valtype(::Type{<:Permutation}, M::Type{<:OperatorProd})
    M = reparameter(M, :id, ID{eltype(M)})
    return OperatorSum{M, idtype(M)}
end
@inline Base.valtype(P::Type{<:Permutation}, M::Type{<:OperatorSum}) = valtype(P, eltype(M))

"""
    (permutation::Permutation)(m::OperatorProd; kwargs...) -> OperatorSum

Permute the operator units of an `OperatorProd` to the descending order according to the table contained in `permutation`.

!!! note
    To use this function, the user must implement a method of `permute`, which computes the result of the permutation of two operator units:
    ```julia
    permute(u₁::OperatorUnit, u₂::OperatorUnit) -> Union{OperatorProd, OperatorSum}
    ```
    Here, `u₁` and `u₂` are two arbitrary operator units contained in `id(m)`.
"""
@inline function (permutation::Permutation)(m::OperatorProd; kwargs...)
    result = zero(permutation, m)
    cache = eltype(result)[m]
    while length(cache) > 0
        current = pop!(cache)
        pos = operatorprodcommuteposition(sequence(current, permutation.table))
        if isnothing(pos)
            add!(result, current)
        else
            left = current[1:pos-1] * value(m)
            right = current[pos+2:end]
            for middle in permute(id(current)[pos], id(current)[pos+1])
                temp = left * middle * right
                isnothing(temp) || push!(cache, temp)
            end
        end
    end
    return result
end
function operatorprodcommuteposition(sequences)
    pos = 1
    while pos < length(sequences)
        (sequences[pos] < sequences[pos+1]) && return pos
        pos += 1
    end
    return nothing
end

"""
    AbstractSubstitution <: Transformation

The substitution transformation.
"""
abstract type AbstractSubstitution <: Transformation end

"""
    AbstractUnitSubstitution{U<:OperatorUnit, S<:OperatorSum} <: AbstractSubstitution

The "unit substitution" transformation, which substitutes evergy `OperatorUnit` in the old quantum operators to a new expression represented by an `OperatorSum`.
"""
abstract type AbstractUnitSubstitution{U<:OperatorUnit, S<:OperatorSum} <: AbstractSubstitution end
@inline function Base.valtype(::Type{<:AbstractUnitSubstitution{U, S}}, M::Type{<:OperatorProd}) where {U<:OperatorUnit, S<:OperatorSum}
    @assert U<:eltype(eltype(M)) "valtype error: mismatched unit transformation."
    V = fulltype(eltype(S), NamedTuple{(:value, :id), Tuple{promote_type(valtype(eltype(S)), valtype(M)), ID{eltype(eltype(S))}}})
    return OperatorSum{V, idtype(V)}
end
@inline Base.valtype(P::Type{<:AbstractUnitSubstitution}, M::Type{<:OperatorSum}) = valtype(P, eltype(M))

"""
    (unitsubstitution::AbstractUnitSubstitution)(m::OperatorProd; kwargs...) -> OperatorSum

Substitute every `OperatorUnit` in an `OperatorProd` with a new `OperatorSum`.
"""
function (unitsubstitution::AbstractUnitSubstitution)(m::OperatorProd; kwargs...)
    return prod(ntuple(i->unitsubstitution(m[i]; kwargs...), Val(rank(m))), init=value(m))
end

"""
    UnitSubstitution{U<:OperatorUnit, S<:OperatorSum, T<:AbstractDict{U, S}} <: AbstractUnitSubstitution{U, S}

A concrete "unit substitution" transformation, which stores every substitution of the old `OperatorUnit`s in its table as a dictionary.
"""
struct UnitSubstitution{U<:OperatorUnit, S<:OperatorSum, T<:AbstractDict{U, S}} <: AbstractUnitSubstitution{U, S}
    table::T
    function UnitSubstitution(table::AbstractDict{<:OperatorUnit, <:OperatorSum})
        new{keytype(table), valtype(table), typeof(table)}(table)
    end
end
(unitsubstitution::UnitSubstitution)(m::OperatorUnit; kwargs...) = unitsubstitution.table[m]

"""
    RankFilter{R} <: Transformation

Rank filter, which filters out the `OperatorPack` with a given rank `R`.
"""
struct RankFilter{R} <: Transformation
    function RankFilter(rank::Int)
        @assert rank>=0 "RankFilter error: the wanted rank must be non-negative."
        new{rank}()
    end
end
@inline Base.valtype(::Type{RankFilter{R}}, M::Type{<:OperatorPack}) where R = reparameter(M, :id, ID{eltype(M), R})
@inline function Base.valtype(R::Type{<:RankFilter}, M::Type{<:OperatorSum})
    V = valtype(R, eltype(M))
    return OperatorSum{V, idtype(V)}
end
@inline rank(rf::RankFilter) = rank(typeof(rf))
@inline rank(::Type{RankFilter{R}}) where R = R
@inline @generated (rf::RankFilter)(m::OperatorPack; kwargs...) = rank(m)==rank(rf) ? :(m) : 0

end #module


text:
using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using DataFrames
using FreqTables
using Distributions

#(l)
function q1()
    #(a)
    rng = Random.seed!(1234)
    A = rand(rng,-5:10,(10,7))
    d = Normal(-2,15)
    B = rand(rng,d,(10,7))
    C = hcat(A[1:5,1:5],B[1:5,6:end])
    D = copy(A)
    D[D .> 0] .= 0
    #(b)
    println(length(A))
    #(c)
    println(length(unique(D)))
    #(d)
    E = reshape(B,(1,:))
    #(e)
    F =cat(A, B, dims=3)
    #(f)
    F = permutedims(F, (3,1,2))
    #(g)
    G = kron(B, C)
    #if I use kron(C, F), there is a dimentional mismatch
    #(h)
    @save "matrixpractice.jld"  A B C D E F G
    #(i)
    @save "firstmatrix.jld"  A B C D
    #(j)
    C = DataFrame(C,:auto)
    CSV.write( "Cmatrix.csv",  C)
    #(k)
    D = DataFrame(D,:auto)
    CSV.write( "Dmatrix.dat",  D)
    return A,B,C,D
end

A,B,C,D = q1()

#(2)
#(f)
function q2(A,B,C)
    #(a)
    AB = zeros(size(A))
    for j = 1:size(AB,2)
        for i = 1:size(AB,1)
            AB[i,j]= A[i,j]*B[i,j]
        end
    end
    AB2 = A.*B
    #{b}
    Cprime = vec(Array(C))
    Cprime = [k for k in Cprime if k <5 && k>=-5]
    Cprime2 = vec(Array(C))
    Cprime2 = filter(t ->  - 5 <= t <= 5, Cprime2)
    #(c)
    X = zeros(15,6,5)
    for t in 1:size(X,3)
        for k in 1:size(X,2)
            if k == 1
                X[:,k,t] = ones(size(X,1))
            elseif k ==2
                p =0.75*(6-t)/5
                d = Binomial(1,p)
                N = size(X,1)
                X[:,k,t] = rand(d,N)
            elseif k ==3
                m = 15 +t -1
                dev = 5*(t-1)
                N = size(X,1)
                d = Normal(m,dev)
                X[:,k,t] = rand(d,N)
            elseif k==4
                m = pi*(6-t)/3
                dev = 1/ℯ
                N = size(X,1)
                d = Normal(m,dev)
                X[:,k,t] = rand(d,N)
            elseif k ==5
                d = Binomial(20,0.6)
                N = size(X,1)
                X[:,k,t] = rand(d,N)
            else
                d = Binomial(20,0.5)
                N = size(X,1)
                X[:,k,t] = rand(d,N)
            end
        end
    end
    #(d)
    β = zeros(6,5)
    for k in 1:size(β,1)
        for t in 1:size(β,2)
            if k == 1
                β[k,t] = 0.75 + 0.25 *t
            elseif k == 2
                β[k,t] = log(t)
            elseif k == 3
                β[k,t] = - sqrt(t)
            elseif k == 4
                β[k,t] = ℯ^t - ℯ^(t+1)
            elseif k == 5
                β[k,t] = t
            else
                β[k,t] = t/3
            end
        end
    end
    #(e)
    Y = zeros(15,5)
    for t in 1:size(X,3)
        d = Normal(0,0.36)
        N =size(X,1)
        ϵ =  rand(d,N)
        Y[:,t]= X[:,:,t]*β[:,t] + ϵ
    end
    return
end
q2(A,B,C)

#(3)
#(g)
function q3()
    df = CSV.read("nlsw88.csv",DataFrame)
    @save "nlsw88.jld"  df
    #(b)
    println("Never married rate: ",size(df[(df.never_married.==1),:],1)/size(df,1))
    println("college graduates rate: ",size(df[(df.collgrad.==1),:],1)/size(df,1))
    #(c)
    freqtable(df.race)
    #(d)
    summarystats = describe(df,:mean, :std, :min, :q25, :median, :q75, :max, :nunique, :nmissing)
    println("Number of missing values of grade is ",summarystats[summarystats.variable.==:grade,:].nmissing)
    summarystats = Array(summarystats)
    #(e)
    cv = freqtable(df.industry,df.occupation)
    #(f)
    df1 = select(df, [:industry,:occupation, :wage])
    gdf = groupby(df1, [:industry,:occupation])
    df_mean = combine(gdf,:wage => mean)
    return
end

q3()

#(4)
#(h)
function q4()
    #(a)
    d = load("firstmatrix.jld")
    A = d["A"]
    B = d["B"]
    C = d["C"]
    D = d["D"]
    #(b)
    function matrixops(A,B)
        #=
        The function does the following calculations:
        (i) AB2: the element-by-element product of the inputs, (ii) AB3: the product A0B, and (iii) AB4: the sum of all the elements of A+B.
        =#
        if size(A)!=size(B)
            throw(error("inputs must have the same size."))
            return
        end
        AB2 = A.*B
        AB3 = A'B
        AB4 = sum(A+B)
        return AB2,AB3, AB4
    end
    #(d)
    println(matrixops(A,B))
    #(f)
    #println(matrixops(C,D))
    #(g)
    df = CSV.read("nlsw88.csv",DataFrame)
    ttl_exp = convert(Array,df.ttl_exp)
    wage = convert(Array,df.wage)
    println(matrixops(ttl_exp,wage))
    return
end

q4()


text:
# Running all testcases and checking their performance

using ODE
using DASSL
const S = Solvers

# Test that test-cases produce the same result.  A bit fiddly as
# failures can come from both this package as well as the used ODE
# solver.

atols = [1e-6]
rtols = atols

typealias TestDict Dict{Symbol, Dict{Solver, Float64}}
tests = TestDict()
# need at least one entry for each IVPTestSuite.tc_all
tests[:hires]     = Dict(S.allsolvers[ODE.ode23s] => 3.3206089)
tests[:vdpol]     = Dict(S.allsolvers[ODE.ode23s] =>  4.610889958728395)
tests[:threebody] = Dict(S.allsolvers[ODE.ode45_dp]  => 2.3352306868067383) # there was a regression at some point, used to be:  2.6330632171040236)
tests[:rober]     = Dict(S.allsolvers[ODE.ode23s] => 0.8920654209711261) # there was a regression at some point, used to be: 1.4917671146318976)
tests[:bruss1d]   = Dict(S.allsolvers[ODE.ode23s] => 4.068350511153728)
tests[:chemakzo]  = Dict(S.allsolvers[DASSL.dasslSolve]  => 4.352199985825764)
tests[:plei]  = Dict(S.allsolvers[ODE.ode45_dp] => 3.4750436027559948)

#@test length(tests)==length(IVPTestSuite.tc_all)

for (name,tc) in IVPTestSuite.tc_all
    println(name)
    for (solver,scd) in tests[name]
        println(solver)
        suite = TestSuite(tc, solver, atols, rtols, [NaN])
        tr = collect(suite)[1]
        re = IVPTestSuite.run_ode_test_throwerror(tr)
        @test abs(re.scd-scd)<1e-3 # pass if less than 0.1% change
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 22cec81f-c503-43de-8297-ec0591e4e6c5
# ╠═╡ show_logs = false
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add("Plots")
		Pkg.add("BenchmarkTools")
		Pkg.add(url="https://github.com/Dale-Black/DistanceTransforms.jl")
	end

	using PlutoUI
	using Plots
	using BenchmarkTools
	using DistanceTransforms
end

# ╔═╡ ed7fea36-60ca-47e0-9046-bf2565fdd41d
md"""
## Import packages
"""

# ╔═╡ afeb0802-60c1-40e0-a9e5-b13ad6af27bc
TableOfContents()

# ╔═╡ 4f2e365c-6e11-490e-933f-c9036ae4e018
md"""
## GPU enabled distance transforms
One of the advantages of a Julia-based distance transform library is how accessible something like GPU programming is for complex algorithms. The `SquaredEuclidean` distance transform is highly parallelizable and is set up to take advantage of the [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl) library and GPU hardware.

Unfortunately, [glassnotebook.io](https://glassnotebook.io) does not yet support GPU compatible static exports, so this documentation, as it stands currently, only shows how to set up a GPU-based distance transform, without any execution.

The code is copy-pastable though for anyone interested in testing it out on their own GPU-enabled computer.
"""

# ╔═╡ f94cf7fa-91de-4fbe-ae52-2063bd01f2ba
md"""
```julia
array = CUDA.CuArray(
	boolean_indicator(
		[
			0 1 1 1 0
			1 1 1 1 1
			1 0 0 0 1
			1 0 0 0 1
			1 0 0 0 1
			1 1 1 1 1
			0 1 1 1 0
		]
	)
)

dt = CuArray{Float32}(undef, size(x))
v = CUDA.ones(Int64, size(x))
z = CUDA.zeros(Float32, size(x) .+ 1)
tfm = DistanceTransforms.SquaredEuclidean(x, dt, v, z)

transform!(x, tfm)
```

This would output a CuArray, like so:
```julia
CUDA.CuArray(
	[
		1.0 0.0 0.0 0.0 1.0
		0.0 0.0 0.0 0.0 0.0
		0.0 1.0 1.0 1.0 0.0
		0.0 4.0 4.0 4.0 0.0
		0.0 1.0 1.0 1.0 0.0
		0.0 0.0 0.0 0.0 0.0
		1.0 0.0 0.0 0.0 1.0
	]
)
```
"""

# ╔═╡ 83901aae-6c82-4d4e-b85d-cfb8fbefa23c
md"""
## GPU enabled distance transform use cases

Distance transforms are ubiquitous in the fields of computer vision and image processing. From object recognition and path planning to deep learning and image segmentation, distance transforms are increasingly useful. DistanceTransforms.jl was created to give developers a simple, consistent API for implementing various distance transforms and to give end-users a seamless way to utilize distance transforms for various tasks, especially GPU-related tasks.

One such example can be seen in the next tutorial on using distance transforms within various [loss functions](link...)
"""

# ╔═╡ Cell order:
# ╟─ed7fea36-60ca-47e0-9046-bf2565fdd41d
# ╠═22cec81f-c503-43de-8297-ec0591e4e6c5
# ╠═afeb0802-60c1-40e0-a9e5-b13ad6af27bc
# ╟─4f2e365c-6e11-490e-933f-c9036ae4e018
# ╟─f94cf7fa-91de-4fbe-ae52-2063bd01f2ba
# ╟─83901aae-6c82-4d4e-b85d-cfb8fbefa23c

max_stars_repo_name: Dale-Black/DistanceTransformsDocs

max_issues_repo_name: Dale-Black/DistanceTransformsDocs

max_forks_repo_name: Dale-Black/DistanceTransformsDocs


text:
@testset "KO Filtering tests:" begin

    # 1D case

    # N=3
    dim = 20
    delta_    = zeros(dim)
    delta_[3] = 1.0

    order = 3
    sigma = 16.0
    ko_filter = Jecco.KO_Filter(order, sigma, dim)

    ko_filter(delta_)

    @test delta_[1:5] ≈ ko_filter.kernel

    # N=5
    dim = 20
    delta_    = zeros(dim)
    delta_[4] = 1.0

    order = 5
    sigma = 64.0
    ko_filter = Jecco.KO_Filter(order, sigma, dim)

    ko_filter(delta_)

    @test delta_[1:7] ≈ ko_filter.kernel

    # 3D case

    dimx = 12
    dimy = 16
    dimz = 10
    delta_         = zeros(dimx, dimy, dimz)
    delta_[:,4,:] .= 1

    order = 5
    sigma = 64.0
    ko_filter = Jecco.KO_Filter{2}(order, sigma, dimx, dimy, dimz)

    ko_filter(delta_)

    myarr = []
    for k in 1:dimz
        for i in 1:dimx
            push!(myarr, delta_[i,1:7,k] ≈ ko_filter.kernel)
        end
    end

    @test all(myarr)
end

@testset "Exponential filtering tests:" begin
    # kernel test
    dim = 49
    γ   = 8.0
    exp_filter   = Jecco.Exp_Filter(γ, dim)
    kernel_bench = [ 1.0, 9.99999999999889e-01, 9.99999999377987e-01, 9.99997121269861e-01,
                     9.88950030743732e-01, 2.22034255471165e-16]
    @test exp_filter.kernel[end-5:end] ≈ kernel_bench

    # 1D case
    dim = 16
    delta_    = zeros(dim)
    delta_[1] = 1.0

    exp_filter1D = Jecco.Exp_Filter{1}(γ, dim)

    f = copy(delta_)

    exp_filter1D(f)

    @test exp_filter1D.fft_plan * f ≈ exp_filter1D.kernel

    # 3D case
    dimx = 12
    dimy = 16
    dimz = 8
    delta3D_ = zeros(dimx, dimy, dimz)
    delta3D_[:,1,:] .= 1

    exp_filter3D = Jecco.Exp_Filter{2}(γ, dimx, dimy, dimz)

    f = copy(delta3D_)

    exp_filter3D(f)

    tmp = exp_filter3D.fft_plan * f
    myarr = []
    for k in 1:dimz
        for i in 1:dimx
            push!(myarr, tmp[i,:,k] ≈ exp_filter3D.kernel)
        end
    end

    @test all(myarr)
end