text:
@mix @columns struct AdditiveGeneticVariance{VA}
    # Field                      | Default  | Flat  | Bounds      | Description
    additiveGeneticVariance::VA  | 2.0      | true  | (0.0, 1e9)  | "Additive genetic variance"
end

"""
Extends CellRule for rules of selection gradient effect

For best performance these should be chained with other
CellRule or following an NeighborhoodRule.
"""
abstract type SelectionGradientRule{R,W} <: CellRule{R,W} end

"""
Extends SelectionGradientRule for selction rules using an heterogenous
exposure layer.
"""
abstract type SelectionGradientMapRule{R,W} <: SelectionGradientRule{R,W} end

DynamicGrids.precalcrules(rule::SelectionGradientMapRule, data) = begin
    precalclayer(layer(rule, data), rule, data)
end

DynamicGrids.precalcrules(rule::SelectionGradientRule, data) = rule

@Exposure @HillCoefficient @AdditiveGeneticVariance struct SelectionGradientSurv{R,W} <: SelectionGradientRule{R,W} end

@inline applyrule(data, rule::SelectionGradientSurv, LC50, args...) = begin
    LC50 > zero(LC50) || return zero(LC50)
    @fastmath LC50 + rule.additiveGeneticVariance * rule.hillcoefficient*(rule.exposure/LC50)^rule.hillcoefficient / ( LC50 *((rule.exposure/LC50)^rule.hillcoefficient +1))
end

@Layers @HillCoefficient @AdditiveGeneticVariance struct SelectionGradientSurvMap{R,W} <: SelectionGradientMapRule{R,W} end

@inline applyrule(data, rule::SelectionGradientSurvMap, LC50, index, args...) = begin
    LC50 > zero(LC50) || return zero(LC50)
    exposure = layer(rule, data, index)
    @fastmath LC50 + rule.additiveGeneticVariance * rule.hillcoefficient*(exposure/LC50)^rule.hillcoefficient / ( LC50 *((exposure/LC50)^rule.hillcoefficient +1))
end

# ######### ----- Phenotype handling single locus genotype
struct LandeVariable{AF,PH}
    alleleFrequency::AF
    phenotype::PH
end
# const LV = LandeVariable

Base.:*(lv::LandeVariable, x::Number) = LandeVariable(x * lv.alleleFrequency, x * lv.phenotype)
Base.:*(x::Number, lv::LandeVariable) = LandeVariable(x * lv.alleleFrequency, x * lv.phenotype)
Base.:+(lv::LandeVariable, x::Number) = LandeVariable(x + lv.alleleFrequency, x + lv.phenotype)
Base.:+(x::Number, lv::LandeVariable) = LandeVariable(x + lv.alleleFrequency, x + lv.phenotype)
Base.:+(lv1::LandeVariable, lv2::LandeVariable) = LandeVariable(lv1.alleleFrequency + lv2.alleleFrequency, lv1.phenotype + lv2.phenotype)
Base.:-(lv::LandeVariable, x::Number) = LandeVariable(x - lv.alleleFrequency, x - lv.phenotype)
Base.:-(x::Number, lv::LandeVariable) = LandeVariable(x - lv.alleleFrequency, x - lv.phenotype)
Base.:-(lv1::LandeVariable, lv2::LandeVariable) = LandeVariable(lv1.alleleFrequency - lv2.alleleFrequency, lv1.phenotype - lv2.phenotype)

Base.zero(::Type{<:LandeVariable{T1,T2}}) where {T1,T2} = LandeVariable(zero(T1), zero(T2))

@Exposure @HillCoefficient @DominanceDegree @DeviationPhenotype struct SelectionGradient1locusSurv{R,W} <: SelectionGradientRule{R,W} end

@inline applyrule(data, rule::SelectionGradient1locusSurv, landeVariable, args...) = begin
    @fastmath LandeVariable(
    landeVariable.alleleFrequency +
    landeVariable.alleleFrequency * (1-landeVariable.alleleFrequency)*(rule.deviationPhenotype + rule.dominanceDegree*(1-2*landeVariable.alleleFrequency))*
    rule.hillcoefficient*(rule.exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((rule.exposure/landeVariable.phenotype)^rule.hillcoefficient +1))
    ,

    landeVariable.phenotype +
        2*landeVariable.alleleFrequency*(1-landeVariable.alleleFrequency) * (rule.deviationPhenotype + rule.dominanceDegree*(1-2*landeVariable.alleleFrequency))^2*
      rule.hillcoefficient*(rule.exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((rule.exposure/landeVariable.phenotype)^rule.hillcoefficient +1)) *
      (1-rule.dominanceDegree*rule.hillcoefficient*(rule.exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((rule.exposure/landeVariable.phenotype)^rule.hillcoefficient +1)))
    )
end

@Layers @HillCoefficient @DominanceDegree @DeviationPhenotype struct SelectionGradient1locusSurvMap{R,W} <: SelectionGradientMapRule{R,W} end

@inline applyrule(data, rule::SelectionGradient1locusSurvMap, landeVariable, index, args...) = begin
    exposure = layer(rule, data, index)
    @fastmath LandeVariable(
        landeVariable.alleleFrequency +
    landeVariable.alleleFrequency * (1-landeVariable.alleleFrequency)*(rule.deviationPhenotype + rule.dominanceDegree*(1-2*landeVariable.alleleFrequency))*
    rule.hillcoefficient*(exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((exposure/landeVariable.phenotype)^rule.hillcoefficient +1)),
    landeVariable.phenotype +
        2*landeVariable.alleleFrequency*(1-landeVariable.alleleFrequency) * (rule.deviationPhenotype + rule.dominanceDegree*(1-2*landeVariable.alleleFrequency))^2*
      rule.hillcoefficient*(exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((exposure/landeVariable.phenotype)^rule.hillcoefficient +1)) *
      (1-rule.dominanceDegree*rule.hillcoefficient*(exposure/landeVariable.phenotype)^rule.hillcoefficient / ( landeVariable.phenotype *((exposure/landeVariable.phenotype)^rule.hillcoefficient +1)))
    )
end


text:
function datetime_intervals end
function sequence_and_random_date_shift end

ext: jl


text:
using JuMP, BilevelJuMP, Cbc

#=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=PMAX1  |                      |
    Csta=0, Cprop=1     |                      |
      P1                |----------------------|
                        |        L12           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      D1                |                      | D2
=#

model = BilevelModel(Cbc.Optimizer)
#=
BilevelJuMP.AbstractBoundedMode{T}
BilevelJuMP.ComplementMode{T}
BilevelJuMP.IndicatorMode{T}
BilevelJuMP.NoMode{T}
BilevelJuMP.ProductMode{T}
BilevelJuMP.SOS1Mode{T}
BilevelJuMP.StrongDualityMode{T}
=#
BilevelJuMP.set_mode(model, BilevelJuMP.IndicatorMode())

PMAX1 = 100.

# equivalent testcase : no_problem
# works as expected
P1 = 40.
D1 = 10.
D2 = 30.
L12 = 35.

# # equivalent testcase : EOD_problem_needs_capping
# # works as expected
# P1 = 100.
# D1 = 10.
# D2 = 30.
# L12 = 35.

# # equivalent testcase : EOD_problem_needs_loss_of_load
# # FAIL : cuts all prod and conso for a cost of 130*e17+40
# # while a solution of cost 90*1e7 exists
# # adding @constraint(Lower(model), e_st == 0), gives the expected result
# P1 = 40.
# D1 = 100.
# D2 = 30.
# L12 = 35.

# # equivalent testcase : RSO_problem
# # FAIL : cuts all prod and conso for a cost of 50*e17+50
# # while a solution of cost 5*1e7+5 exists
# # adding @constraint(Lower(model), e_st == 5), gives the expected result
# P1 = 50.
# D1 = 10.
# D2 = 40.
# L12 = 35.

#upper
@variable(Upper(model), lolmin_st, lower_bound=0., upper_bound=D1+D2)
@variable(Upper(model), lol_n1st, lower_bound=0., upper_bound=D1)
@variable(Upper(model), lol_n2st, lower_bound=0., upper_bound=D2)

@variable(Upper(model), emin_st, lower_bound=0., upper_bound=P1)
@variable(Upper(model), e_n1st, lower_bound=0., upper_bound=P1)
@variable(Upper(model), penr_n1st, lower_bound=0., upper_bound=P1)#modif
@variable(Upper(model), plim_n1st, lower_bound=0., upper_bound=PMAX1)#modif
@variable(Upper(model), islimited_n1st, binary=true)
@variable(Upper(model), plimxislim_n1st, lower_bound=0., upper_bound=PMAX1)#modif

#lower
@variable(Lower(model), lol_st, lower_bound=0., upper_bound=D1+D2)

@variable(Lower(model), e_st, lower_bound=0., upper_bound=P1)
@variable(Lower(model), penr_st, lower_bound=0., upper_bound=P1)

# Upper
@objective(Upper(model), Min, 1e7*lolmin_st+emin_st)
@constraints(Upper(model), begin
    -35. <= 0.5*penr_n1st +0.5*lol_n1st -0.5*D1 -0.5*lol_n2st +0.5*D2 <= 35. #RSO
    lol_st == lol_n2st + lol_n1st
    e_n1st == P1 - penr_n1st
    #plimxislim_n1st = plim_n1st * islimited_n1st
    plimxislim_n1st <= plim_n1st
    plimxislim_n1st <= PMAX1*islimited_n1st
    plim_n1st <= PMAX1*(1-islimited_n1st) + plimxislim_n1st
    # penr_n1st min(P1, plim_n1st)
    #penr_n1st <= P1 #as bound
    penr_n1st <= plim_n1st
    penr_n1st == plimxislim_n1st + P1*(1 - islimited_n1st)
    penr_st == penr_n1st
    e_st == e_n1st
end)

# Lower
@objective(Lower(model), Min, 1e7*lol_st+e_st)
@constraints(Lower(model), begin
    P1-e_st == D1 + D2 - lol_st
    emin_st <= e_st
    lolmin_st <= lol_st
end)

optimize!(model)
# optimize!(model, lower_prob="low", upper_prob="up", bilevel_prob="bil", solver_prob="solv")

println("obj lower : ", objective_value(Lower(model)))
println("obj upper : ", objective_value(model))

println("lolmin : ", value(lolmin_st))
println("emin : ", value(emin_st))
println("lol : ", value(lol_st))
println("e : ", value(e_st))

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using IntervalArithmetic
using Base.Test

@test IntervalArithmetic.round_expr(:(a + b), RoundDown) == :($(Expr(:escape, :a)) + $(Expr(:escape, :b)) + $(RoundDown))

@test IntervalArithmetic.round_expr(:(sin(a)), RoundUp) ==
    :( sin( $(Expr(:escape, :a)), $(RoundUp) ) )


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools
using StaticArrays
pdb_filepath = "1ake.pdb"

comp = System(pdb_filepath)
dfs = DataFrameSystem(comp)

comp2 = System(pdb_filepath)
for x in AtomIterator(comp2)
    setPosition(x, SVector{3,Float64}(x.position_[1]/sqrt(3),x.position_[2]/sqrt(3),x.position_[3]/sqrt(3)))
end

rmsd(comp, comp2,false)

elapsed = @elapsed rmsd(comp,comp2,false)

println(elapsed)
@show @benchmark rmsd(comp,comp2,false)


text:
const PrecisionPiece{T} = Tuple{Vector{Int64}, Vector{Int64}, Vector{T}}

abstract type AbstractVecchiaConfig{H,D,F} end

# TODO (cg 2021/04/25 13:06): should these fields chunksize and blockrank be in
# here? Arguably the are redundant and encoded in the data/pts/condix values.
# And having them sort of provides a dangerously easy option to not check and
# make sure what those sizes really need to be.
struct VecchiaConfig{H,D,F} <: AbstractVecchiaConfig{H,D,F}
  chunksize::Int64
  blockrank::Int64
  kernel::F
  data::Vector{Matrix{H}}
  pts::Vector{Vector{SVector{D, Float64}}}
  condix::Vector{Vector{Int64}}
end

function Base.display(V::VecchiaConfig)
  println("Vecchia configuration with:")
  println("chunksize:  $(V.chunksize)")
  println("block rank: $(V.blockrank)")
  println("data size:  $(sum(x->size(x,1), V.data))")
  println("nsamples:   $(size(V.data[1], 2))")
end

# TODO (cg 2021/04/25 13:06): should these fields chunksize and blockrank be in
# here? Arguably the are redundant and encoded in the data/pts/condix values.
# And having them sort of provides a dangerously easy option to not check and
# make sure what those sizes really need to be.
struct ScalarVecchiaConfig{H,D,F} <: AbstractVecchiaConfig{H,D,F}
  chunksize::Int64
  blockrank::Int64
  kernel::F
  data::Vector{Matrix{H}}
  pts::Vector{Vector{Float64}}
  condix::Vector{Vector{Int64}}
end

function Base.display(V::ScalarVecchiaConfig)
  println("Scalarized Vecchia configuration with:")
  println("chunksize:  $(V.chunksize)")
  println("block rank: $(V.blockrank)")
  println("data size:  $(sum(x->size(x,1), V.data))")
  println("nsamples:   $(size(V.data[1], 2))")
end

struct RCholesky{T}
  diagonals::Vector{UpperTriangular{T,Matrix{T}}}
  odiagonals::Vector{Matrix{T}}
  condix::Vector{Vector{Int64}}
  idxs::Vector{UnitRange{Int64}}
  is_instantiated::Bool
end

# TODO (cg 2022/05/30 12:10): make this more information.
Base.display(U::RCholesky{T}) where{T} = println("RCholesky{$T}")
Base.display(Ut::Adjoint{T,RCholesky{T}}) where{T} = println("Adjoint{$T, RCholesky{$T}}")

struct MemoizedKernel{F,T}
  cache::Dict{Tuple{UInt64, UInt64},T}
  phash::UInt64
  kernel::F
end

# Not good code or anything. Just a quick and dirty way to make a Vecchia object
# with a KD-tree to choose the conditioning points.
function kdtreeconfig(data, pts, chunksize, blockrank, kfun)
  (data isa Vector) && return kdtreeconfig(hcat(data), pts, chunksize, blockrank, kfun)
  # Make a KDTree of the points with a certain leaf size
  tree    = KDTree(pts, leafsize=chunksize)
  # re-order the data accordingly:
  _d      = Dict(zip(pts, eachrow(data)))
  data_p  = [_d[xj] for xj in tree.data]
  # Get the data in chunks of size $chunksize.
  pts_out = map(copy, Iterators.partition(tree.data, chunksize))
  # TODO (cg 2022/04/21 16:17): The inner function here (vec_of_...), really
  # could use some optimization and is kind of just a stand-in at the moment. So
  # long as this constructor isn't in the hot loop of an optimizer it's probably
  # fine as-is, but just a note that there is some perf on the table here.
  dat_out = map(vec_of_vecs_to_matrows, Iterators.partition(data_p, chunksize))
  # re-order the centroids of THOSE leaves.
  centroids = map(x->sum(x)/length(x), pts_out)
  c_ix_dict = Dict(zip(centroids, eachindex(centroids)))
  c_tree    = KDTree(centroids)
  perm      = [c_ix_dict[x] for x in c_tree.data]
  @assert sort(perm) == collect(eachindex(centroids)) "debug test."
  pts_out   = pts_out[perm]
  dat_out   = dat_out[perm]
  # Create the conditioning meta-indices for the chunks.
  condix  = map(j->cond_ixs(j,blockrank), eachindex(pts_out))
  (H,D,F) = (eltype(data), length(first(pts)), typeof(kfun))
  VecchiaConfig{H,D,F}(min(chunksize, length(first(pts_out))),
                       min(blockrank, length(pts_out)),
                       kfun, dat_out, pts_out, condix)
end

# Even less good code, but here's a simple connection between Vecchia and the
# H-matrix with global Nystrom that I've worked on earlier: add a single
# collection of "global" points to the conditioning set of all subsequent
# points. Pretty easy to do with a couple extra tree constructions.
function nystrom_kdtreeconfig(data, pts, chunksize, blockrank, kfun, nys_size)
  (data isa Vector) && return nystrom_kdtreeconfig(hcat(data), pts, chunksize,
                                                   blockrank, kfun, nys_size)
  # Make a KDTree of the points with a certain leaf size
  tree   = KDTree(pts, leafsize=chunksize)
  _d     = Dict(zip(pts, eachrow(data)))
  data_p = reshape(reduce(vcat, [_d[xj] for xj in tree.data]), :, size(data,2))
  # Pick the "nystrom points" in the simplest way:
  nys_inds = 1:div(length(pts), nys_size):length(pts)
  nys_pts  = tree.data[nys_inds]
  # now for the rest of the points, do the ordering above:
  _cfg     = kdtreeconfig(data_p, tree.data[setdiff(1:length(pts), nys_inds)],
            chunksize, blockrank, kfun)
  # extract and modify the fields we actually want:
  pts_out  = vcat([nys_pts], _cfg.pts)
  nys_dat  = reshape(reduce(vcat, data_p[nys_inds]), :, size(data, 2))
  dat_out  = pushfirst!(_cfg.data, nys_dat)
  condix   = map(j->cond_ixs(j,blockrank), eachindex(pts_out))
  for j in 2:length(condix)
    cixj = condix[j]
    if !in(1, cixj)
      pushfirst!(cixj, 1)
    end
  end
  (H,D,F) = (eltype(data), length(first(pts)), typeof(kfun))
  VecchiaConfig{H, D,F}(min(chunksize+nys_size, length(first(pts_out))),
                        min(blockrank+1, length(pts_out)),
                        kfun, dat_out, pts_out, condix)
end

function scalarize(v::VecchiaConfig{H,D,F}, scalarized_kernel::G) where{H,D,F,G}
  scalarized_pts = map(x->reduce(vcat, x), v.pts)
  ScalarVecchiaConfig{H,D,G}(v.chunksize,
                             v.blockrank,
                             scalarized_kernel,
                             v.data,
                             scalarized_pts,
                             v.condix)
end

function MemoizedKernel(f, sample_x, p)
  f11 = f(sample_x, sample_x, p)
  ph  = hash(p)
  cache = Dict{Tuple{UInt64, UInt64},typeof(f11)}()
  MemoizedKernel(cache, ph, f)
end

function (k::MemoizedKernel{F,T})(x, y, p) where{F,T}
  _key = (hash(x), hash(y))
  if haskey(k.cache, _key)
    return k.cache[_key]
  else
    _val = k.kernel(x,y,p)
    k.cache[_key] = _val
    return _val
  end
end


text:
using UnsteadyFlowSolvers

alphadef = EldUpDef(45. *pi/180, 0.4, 0.8)

hdef = ConstDef(0.)

udef = ConstDef(1.)

full_kinem = KinemDef(alphadef, hdef, udef)

pvt = 0.0

geometry = "FlatPlate"

lespcrit = [0.11;]

surf = TwoDSurf(geometry, pvt, full_kinem, lespcrit)

curfield = TwoDFlowField()

dtstar = find_tstep(alphadef)

t_tot = 9.

nsteps = Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/18.

#delvort = delSpalart(500, 12, 1e-5)
delvort = delNone()

mat, surf, curfield = ldvm(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

makeVortPlots2D()

makeForcePlots2D()

cleanWrite()

max_stars_repo_name: MaximeUofG/UNSflowTutorials

max_issues_repo_name: MaximeUofG/UNSflowTutorials

max_forks_repo_name: MaximeUofG/UNSflowTutorials


text:
include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))
include(joinpath(LOOPVECBENCHDIR, "loadsharedlibs.jl"))

using BenchmarkTools, SharedArrays
struct SizedResults{V <: AbstractVector} <: AbstractMatrix{String}
    results::Matrix{Float64}
    sizes::V
end
function Base.size(sr::SizedResults)
    M, N = size(sr.results)
    N, M+1
end
struct BenchmarkResult{V}
    tests::Vector{String}
    sizedresults::SizedResults{V}
end
function BenchmarkResult(results, tests, sizes)
    ntests = length(tests); nsizes = length(sizes)
    BenchmarkResult(
        append!(["Size"], tests),
        SizedResults(results, sizes)
    )
end
function Base.getindex(br::SizedResults, row, col)
    col == 1 ? string(br.sizes[row]) : string(br.results[col - 1, row])
end
Base.setindex!(br::BenchmarkResult, v, i...) = br.sizedresults.results[i...] = v
function Base.vcat(br1::BenchmarkResult, br2::BenchmarkResult)
    BenchmarkResult(
        br1.tests,
        SizedResults(
            hcat(br1.sizedresults.results, br2.sizedresults.results),
            vcat(br1.sizedresults.sizes, br2.sizedresults.sizes)
        )
    )
end

tothreetuple(i::Int) = (i,i,i)
tothreetuple(i::NTuple{3,Int}) = i
function matmul_bench!(br, C, A, B, i)
    M, N = size(C); K = size(B,1);
    n_gflop = M*K*N*2e-9
    Cblas = A * B
    br[1,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
    @assert C ≈ Cblas "LoopVec gemm wrong?"; fill!(C, NaN)
    br[2,i] = n_gflop / @belapsed jgemm!($C, $A, $B)
    @assert C ≈ Cblas "Julia gemm wrong?"; fill!(C, NaN)
    br[3,i] = n_gflop / @belapsed cgemm!($C, $A, $B)
    @assert C ≈ Cblas "Clang gemm wrong?"; fill!(C, NaN)
    br[4,i] = n_gflop / @belapsed fgemm!($C, $A, $B)
    @assert C ≈ Cblas "Fort gemm wrong?"; fill!(C, NaN)
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icgemm!($C, $A, $B)
        @assert C ≈ Cblas "icc gemm wrong?"; fill!(C, NaN)
        br[6,i] = n_gflop / @belapsed ifgemm!($C, $A, $B)
        @assert C ≈ Cblas "ifort gemm wrong?"; fill!(C, NaN)
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed egemm!($C, $A, $B)
    @assert C ≈ Cblas "eigen gemm wrong?"; fill!(C, NaN)
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed iegemm!($C, $A, $B)
    @assert C ≈ Cblas "i-eigen gemm wrong?"; fill!(C, NaN)

    br[7+2INTEL_BENCH,i] = n_gflop / @belapsed fgemm_builtin!($C, $A, $B)
    @assert C ≈ Cblas "Fort builtin gemm wrong?"; fill!(C, NaN)
    if INTEL_BENCH
        br[8+2INTEL_BENCH,i] = n_gflop / @belapsed ifgemm_builtin!($C, $A, $B)
        @assert C ≈ Cblas "ifort builtin gemm wrong?"; fill!(C, NaN)
    end
    br[8+3INTEL_BENCH,i] = n_gflop / @belapsed gemmopenblas!($C, $A, $B);
    @assert C ≈ Cblas "OpenBLAS gemm wrong?"
    if MKL_BENCH
        br[9+3INTEL_BENCH,i] = n_gflop / @belapsed gemmmkl!($C, $A, $B)
        @assert C ≈ Cblas "MKL gemm wrong?"
    end
    # br[12,i] = n_gflop / @belapsed gemmavx!($C, $A, $B)
end
function A_mul_B_bench!(br, s, i)
    M, K, N = tothreetuple(s)
    C = Matrix{Float64}(undef, M, N)
    A = rand(M, K)
    B = rand(K, N)
    matmul_bench!(br, C, A, B, i)
end
function A_mul_B_bench!(br, s, i)
    M, K, N = tothreetuple(s)
    C = Matrix{Float64}(undef, M, N)
    A = rand(M, K)
    B = rand(K, N)
    matmul_bench!(br, C, A, B, i)
end
function A_mul_Bt_bench!(br, s, i)
    M, K, N = tothreetuple(s)
    C = Matrix{Float64}(undef, M, N)
    A = rand(M, K)
    B = rand(N, K)'
    matmul_bench!(br, C, A, B, i)
end
function At_mul_B_bench!(br, s, i)
    M, K, N = tothreetuple(s)
    C = Matrix{Float64}(undef, M, N)
    A = rand(K, M)'
    B = rand(K, N)
    matmul_bench!(br, C, A, B, i)
end
function At_mul_Bt_bench!(br, s, i)
    M, K, N = tothreetuple(s)
    C = Matrix{Float64}(undef, M, N)
    A = rand(K, M)'
    B = rand(N, K)'
    matmul_bench!(br, C, A, B, i)
end

function dot_bench!(br, s, i)
    a = rand(s); b = rand(s);
    dotblas = dot(a, b)
    n_gflop = s * 2e-9
    br[1,i] = n_gflop / @belapsed jdotavx($a, $b)
    @assert jdotavx(a,b) ≈ dotblas "LoopVec dot wrong?"
    br[2,i] = n_gflop / @belapsed jdot($a, $b)
    @assert jdot(a,b) ≈ dotblas "Julia dot wrong?"
    br[3,i] = n_gflop / @belapsed cdot($a, $b)
    @assert cdot(a,b) ≈ dotblas "Clang dot wrong?"
    br[4,i] = n_gflop / @belapsed fdot($a, $b)
    @assert fdot(a,b) ≈ dotblas "Fort dot wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icdot($a, $b)
        @assert icdot(a,b) ≈ dotblas "icc dot wrong?"
        br[6,i] = n_gflop / @belapsed ifdot($a, $b)
        @assert ifdot(a,b) ≈ dotblas "ifort dot wrong?"
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed edot($a, $b)
    @assert edot(a,b) ≈ dotblas "eigen dot wrong?"
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed iedot($a, $b)
    @assert iedot(a,b) ≈ dotblas "i-eigen dot wrong?"
    # br[9,i] = n_gflop / @belapsed dot($a, $b)
end
function selfdot_bench!(br, s, i)
    a = rand(s); b = rand(s);
    dotblas = dot(a, a)
    n_gflop = s * 2e-9
    br[1,i] = n_gflop / @belapsed jselfdotavx($a)
    @assert jselfdotavx(a) ≈ dotblas "LoopVec dot wrong?"
    br[2,i] = n_gflop / @belapsed jselfdot($a)
    @assert jselfdot(a) ≈ dotblas "Julia dot wrong?"
    br[3,i] = n_gflop / @belapsed cselfdot($a)
    @assert cselfdot(a) ≈ dotblas "Clang dot wrong?"
    br[4,i] = n_gflop / @belapsed fselfdot($a)
    @assert fselfdot(a) ≈ dotblas "Fort dot wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icselfdot($a)
        @assert cselfdot(a) ≈ dotblas "icc dot wrong?"
        br[6,i] = n_gflop / @belapsed ifselfdot($a)
        @assert fselfdot(a) ≈ dotblas "ifort dot wrong?"
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed eselfdot($a)
    @assert eselfdot(a) ≈ dotblas "eigen dot wrong?"
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed ieselfdot($a)
    @assert ieselfdot(a) ≈ dotblas "i-eigen dot wrong?"
    # br[9,i] = n_gflop / @belapsed dot($a, $a)
end

totwotuple(i::Int) = (i,i)
totwotuple(i::Tuple{Int,Int}) = i
function gemv_bench!(br, x, A, y, i)
    M, N = size(A)
    n_gflop = M*N * 2e-9
    xblas = A * y
    br[1,i] = n_gflop / @belapsed jgemvavx!($x, $A, $y)
    @assert x ≈ xblas "LoopVec wrong?"; fill!(x, NaN);
    br[2,i] = n_gflop / @belapsed jgemv!($x, $A, $y)
    @assert x ≈ xblas "Julia wrong?"; fill!(x, NaN);
    br[3,i] = n_gflop / @belapsed cgemv!($x, $A, $y)
    @assert x ≈ xblas "Clang wrong?"; fill!(x, NaN);
    br[4,i] = n_gflop / @belapsed fgemv!($x, $A, $y)
    @assert x ≈ xblas "Fort wrong?"; fill!(x, NaN);
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icgemv!($x, $A, $y)
        @assert x ≈ xblas "icc wrong?"; fill!(x, NaN);
        br[6,i] = n_gflop / @belapsed ifgemv!($x, $A, $y)
        @assert x ≈ xblas "ifort wrong?"; fill!(x, NaN);
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed egemv!($x, $A, $y)
    @assert x ≈ xblas "eigen wrong?"; fill!(x, NaN);
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed iegemv!($x, $A, $y)
    @assert x ≈ xblas "i-eigen wrong?"; fill!(x, NaN);
    br[7+2INTEL_BENCH,i] = n_gflop / @belapsed fgemv_builtin!($x, $A, $y)
    @assert x ≈ xblas "Fort wrong?"; fill!(x, NaN);
    br[8+2INTEL_BENCH,i] = n_gflop / @belapsed ifgemv_builtin!($x, $A, $y)
    @assert x ≈ xblas "ifort wrong?"; fill!(x, NaN);
    br[9+2INTEL_BENCH,i] = n_gflop / @belapsed dgemvopenblas!($x, $A, $y)
    @assert x ≈ xblas "gemvopenblas wrong?"; fill!(x, NaN);
    br[10+2INTEL_BENCH,i] = n_gflop / @belapsed dgemvmkl!($x, $A, $y)
    @assert x ≈ xblas "gemvmkl wrong?"; fill!(x, NaN);
end
function A_mul_vb_bench!(br, s, i)
    M, N = totwotuple(s)
    x = Vector{Float64}(undef, M);
    A = rand(M, N);
    y = rand(N);
    gemv_bench!(br, x, A, y, i)
end
function At_mul_vb_bench!(br, s, i)
    M, N = totwotuple(s)
    x = Vector{Float64}(undef, M);
    A = rand(N, M)';
    y = rand(N);
    gemv_bench!(br, x, A, y, i)
end

function dot3_bench!(br, s, i)
    M, N = totwotuple(s)
    x = rand(M); A = rand(M, N); y = rand(N);
    dotblas = dot(x, A, y)
    n_gflop = M*N * 3e-9
    br[1,i] = n_gflop / @belapsed jdot3avx($x, $A, $y)
    @assert jdot3avx(x, A, y) ≈ dotblas "LoopVec dot wrong?"
    br[2,i] = n_gflop / @belapsed jdot3($x, $A, $y)
    @assert jdot3(x, A, y) ≈ dotblas "Julia dot wrong?"
    br[3,i] = n_gflop / @belapsed cdot3($x, $A, $y)
    @assert cdot3(x, A, y) ≈ dotblas "Clang dot wrong?"
    br[4,i] = n_gflop / @belapsed fdot3($x, $A, $y)
    @assert fdot3(x, A, y) ≈ dotblas "Fort dot wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icdot3($x, $A, $y)
        @assert icdot3(x, A, y) ≈ dotblas "icc dot wrong?"
        br[6,i] = n_gflop / @belapsed ifdot3($x, $A, $y)
        @assert ifdot3(x, A, y) ≈ dotblas "ifort dot wrong?"
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed edot3($x, $A, $y)
    @assert edot3(x, A, y) ≈ dotblas "eigen dot wrong?"
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed iedot3($x, $A, $y)
    @assert iedot3(x, A, y) ≈ dotblas "c-eigen dot wrong?"
    br[7+2INTEL_BENCH,i] = n_gflop / @belapsed dot($x, $A, $y)
end
# BLAS.set_num_threads(1)
function sse!(Xβ, y, X, β)
    dgemvmkl!(copyto!(Xβ, y), X, β, 1.0, -1.0)
    jdot(Xβ, Xβ)
end
sse_totwotuple(s::NTuple{2}) = s
sse_totwotuple(s::Integer) = ((3s) >> 1, s >> 1)

function sse_bench!(br, s, i)
    N, P = sse_totwotuple(s)
    y = rand(N); β = rand(P);
    X = randn(N, P);
    Xβ = similar(y);
    lpblas = sse!(Xβ, y, X, β)
    n_gflop = 2e-9*(P*N + 2N)
    br[1,i] = n_gflop / @belapsed jOLSlp_avx($y, $X, $β)
    @assert jOLSlp_avx(y, X, β) ≈ lpblas "LoopVec wrong?"
    br[2,i] = n_gflop / @belapsed jOLSlp($y, $X, $β)
    @assert jOLSlp(y, X, β) ≈ lpblas "Julia wrong?"
    br[3,i] = n_gflop / @belapsed cOLSlp($y, $X, $β)
    @assert cOLSlp(y, X, β) ≈ lpblas "Clang wrong?"
    br[4,i] = n_gflop / @belapsed fOLSlp($y, $X, $β)
    @assert fOLSlp(y, X, β) ≈ lpblas "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icOLSlp($y, $X, $β)
        @assert icOLSlp(y, X, β) ≈ lpblas "icc wrong?"
        br[6,i] = n_gflop / @belapsed ifOLSlp($y, $X, $β)
        @assert ifOLSlp(y, X, β) ≈ lpblas "ifort wrong?"
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed eOLSlp($y, $X, $β)
    @assert eOLSlp(y, X, β) ≈ lpblas "eigen wrong?"
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed ieOLSlp($y, $X, $β)
    @assert ieOLSlp(y, X, β) ≈ lpblas "i-eigen wrong?"
    if MKL_BENCH
        br[7+2INTEL_BENCH,i] = n_gflop / @belapsed sse!($Xβ, $y, $X, $β)
        @assert sse!(Xβ, y, X, β) ≈ lpblas "MKL wrong?"
    end
end

function exp_bench!(br, s, i)
    a = rand(s); b = similar(a)
    n_gflop = 1e-9*s # not really gflops
    br[1,i] = n_gflop / @belapsed @turbo @. $b = exp($a)
    baseb = copy(b)
    br[2,i] = n_gflop / @belapsed @. $b = exp($a)
    @assert b ≈ baseb "LoopVec wrong?"
    br[3,i] = n_gflop / @belapsed cvexp!($b, $a)
    @assert b ≈ baseb "Clang wrong?"
    br[4,i] = n_gflop / @belapsed fvexp!($b, $a)
    @assert b ≈ baseb "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icvexp!($b, $a)
        @assert b ≈ baseb "icc wrong?"
        br[6,i] = n_gflop / @belapsed ifvexp!($b, $a)
        @assert b ≈ baseb "ifort wrong?"
    end
end

function aplusBc_bench!(br, s, i)
    M, N = totwotuple(s)
    a = rand(M); B = rand(M,N); c = rand(N);
    c′ = c'; D = similar(B)
    n_gflop = 2e-9 * M*N
    br[1,i] = n_gflop / @belapsed @turbo @. $D = $a + $B * $c′
    Dcopy = copy(D); fill!(D, NaN);
    br[2,i] = n_gflop / @belapsed @. $D = $a + $B * $c′
    @assert D ≈ Dcopy "LoopVec wrong?"
    br[3,i] = n_gflop / @belapsed caplusBc!($D, $a, $B, $c)
    @assert D ≈ Dcopy "Clang wrong?"; fill!(D, NaN);
    br[4,i] = n_gflop / @belapsed faplusBc!($D, $a, $B, $c)
    @assert D ≈ Dcopy "Fort wrong?"; fill!(D, NaN);
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icaplusBc!($D, $a, $B, $c)
        @assert D ≈ Dcopy "icc wrong?"; fill!(D, NaN);
        br[6,i] = n_gflop / @belapsed ifaplusBc!($D, $a, $B, $c)
        @assert D ≈ Dcopy "ifort wrong?"; fill!(D, NaN);
    end
    br[7,i] = n_gflop / @belapsed eaplusBc!($D, $a, $B, $c)
    @assert D ≈ Dcopy "eigen wrong?"; fill!(D, NaN);
    br[8,i] = n_gflop / @belapsed ieaplusBc!($D, $a, $B, $c)
    @assert D ≈ Dcopy "i-eigen wrong?"; fill!(D, NaN);
end

function AplusAt_bench!(br, s, i)
    A = rand(s,s); B = similar(A)
    n_gflop = 1e-9*s^2
    br[1,i] = n_gflop / @belapsed @turbo @. $B = $A + $A'
    baseB = copy(B); fill!(B, NaN);
    br[2,i] = n_gflop / @belapsed @. $B = $A + $A'
    @assert B ≈ baseB "LoopVec wrong?"
    br[3,i] = n_gflop / @belapsed cAplusAt!($B, $A)
    @assert B ≈ baseB "Clang wrong?"; fill!(B, NaN);
    br[4,i] = n_gflop / @belapsed fAplusAt!($B, $A)
    @assert B ≈ baseB "Fort wrong?"; fill!(B, NaN);
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icAplusAt!($B, $A)
        @assert B ≈ baseB "icc wrong?"; fill!(B, NaN);
        br[6,i] = n_gflop / @belapsed ifAplusAt!($B, $A)
        @assert B ≈ baseB "ifort wrong?"; fill!(B, NaN);
    end
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed eAplusAt!($B, $A)
    @assert B ≈ baseB "eigen wrong?"; fill!(B, NaN);
    br[6+2INTEL_BENCH,i] = n_gflop / @belapsed ieAplusAt!($B, $A)
    @assert B ≈ baseB "i-eigen wrong?"; fill!(B, NaN);
    br[7+2INTEL_BENCH,i] = n_gflop / @belapsed fAplusAt_builtin!($B, $A)
    @assert B ≈ baseB "Fort-builtin wrong?"; fill!(B, NaN);
    if INTEL_BENCH
        br[8+2INTEL_BENCH,i] = n_gflop / @belapsed ifAplusAt_builtin!($B, $A)
        @assert B ≈ baseB "ifort-builtin wrong?"; fill!(B, NaN);
    end
end

function randomaccess_bench!(br, s, i)
    A, C = totwotuple(s)
    P = rand(A, C) .+= 0.5;
    basis = rand(1:C, A, C);
    coefs = randn(C);
    n_gflop = 1e-9*(A*C + C)
    p = randomaccess(P, basis, coefs);
    br[1,i] = n_gflop / @belapsed randomaccessavx($P, $basis, $coefs)
    @assert p ≈ randomaccessavx(P, basis, coefs) "LoopVec wrong?"
    br[2,i] = n_gflop / @belapsed  randomaccess($P, $basis, $coefs)
    br[3,i] = n_gflop / @belapsed crandomaccess($P, $basis, $coefs)
    @assert p ≈ crandomaccess(P, basis, coefs) "Clang wrong?"
    br[4,i] = n_gflop / @belapsed frandomaccess($P, $basis, $coefs)
    @assert p ≈ frandomaccess(P, basis, coefs) "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icrandomaccess($P, $basis, $coefs)
        @assert p ≈ icrandomaccess(P, basis, coefs) "icc wrong?"
        br[6,i] = n_gflop / @belapsed ifrandomaccess($P, $basis, $coefs)
        @assert p ≈ ifrandomaccess(P, basis, coefs) "ifort wrong?"
    end
end

function logdettriangle_bench!(br, s, i)
    S = randn(s, 2s)
    U = cholesky(Symmetric(S * S')).U
    n_gflop = 1e-9*s
    ld = logdet(U)
    br[1,i] = n_gflop / @belapsed jlogdettriangleavx($U)
    @assert ld ≈ jlogdettriangleavx(U) "LoopVec wrong?"
    br[2,i] = n_gflop / @belapsed jlogdettriangle($U)
    @assert ld ≈ jlogdettriangle(U) "Julia wrong?"
    br[3,i] = n_gflop / @belapsed clogdettriangle($U)
    @assert ld ≈ clogdettriangle(U) "Clang wrong?"
    br[4,i] = n_gflop / @belapsed flogdettriangle($U)
    @assert ld ≈ flogdettriangle(U) "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed iclogdettriangle($U)
        @assert ld ≈ iclogdettriangle(U) "icc wrong?"
        br[6,i] = n_gflop / @belapsed iflogdettriangle($U)
        @assert ld ≈ iflogdettriangle(U) "ifort wrong?"
    end
    # br[7,i] = n_gflop / @belapsed elogdettriangle($U)
    # @assert ld ≈ elogdettriangle(U) "eigen wrong?"; fill!(B, NaN);
    # br[8,i] = n_gflop / @belapsed ielogdettriangle($U)
    # @assert ld ≈ ielogdettriangle(U) "i-eigen wrong?"; fill!(B, NaN);
    br[5+2INTEL_BENCH,i] = n_gflop / @belapsed logdet($U)
end

function filter2d_bench_run!(br, s, i, K)
    A = rand(s + 2, s + 2)
    B = OffsetArray(similar(A, (s,s)), 1, 1)
    Mk, Nk = size(K)
    n_gflop = 1e-9 * (2Mk * Nk - 1) * s^2
    br[1,i] = n_gflop / @belapsed filter2davx!($B, $A, $K)
    Bcopy = copy(B); fill!(B, NaN);
    br[2,i] = n_gflop / @belapsed filter2d!($B, $A, $K)
    @assert B ≈ Bcopy "LoopVec wrong?"
    br[3,i] = n_gflop / @belapsed cfilter2d!($B, $A, $K)
    @assert B ≈ Bcopy "Clang wrong?"
    br[4,i] = n_gflop / @belapsed ffilter2d!($B, $A, $K)
    @assert B ≈ Bcopy "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icfilter2d!($B, $A, $K)
        @assert B ≈ Bcopy "icc wrong?"
        br[6,i] = n_gflop / @belapsed iffilter2d!($B, $A, $K)
        @assert B ≈ Bcopy "ifort wrong?"
    end
end

function filter2dunrolled_bench_run!(br, s, i, K)
    A = rand(s + 2, s + 2)
    B = OffsetArray(similar(A, (s,s)), 1, 1)
    Mk, Nk = size(K)
    n_gflop = 1e-9 * (2Mk * Nk - 1) * s^2
    br[1,i] = n_gflop / @belapsed filter2dunrolledavx!($B, $A, $K)
    Bcopy = copy(B); fill!(B, NaN);
    br[2,i] = n_gflop / @belapsed filter2dunrolled!($B, $A, $K)
    @assert B ≈ Bcopy "LoopVec wrong?"
    br[3,i] = n_gflop / @belapsed cfilter2dunrolled!($B, $A, $K)
    @assert B ≈ Bcopy "Clang wrong?"
    br[4,i] = n_gflop / @belapsed ffilter2dunrolled!($B, $A, $K)
    @assert B ≈ Bcopy "Fort wrong?"
    if INTEL_BENCH
        br[5,i] = n_gflop / @belapsed icfilter2dunrolled!($B, $A, $K)
        @assert B ≈ Bcopy "icc wrong?"
        br[6,i] = n_gflop / @belapsed iffilter2dunrolled!($B, $A, $K)
        @assert B ≈ Bcopy "ifort wrong?"
    end
end


text:
#!/usr/bin/env julia

VERSION >= v"0.4" && __precompile__()

module ColorSchemesInvert

using ColorTypes, ColorSchemes, Colors

export invert, convertToScheme

"""
    invert(cscheme, c)

Compute the percentage value of the colors in cscheme.

# Examples
```julia-repl
    julia> invert(ColorSchemes.leonardo, RGB(1,0,0))
    0.625…
    julia> invert([RGB(0,0,0), RGB(1,1,1)], RGB(.5,.5,.5))
    0.543…
    julia> cs = linspace(RGB(0,0,0), RGB(1,1,1),5)
    julia> invert(cs, cs[3])
    0.500
```
"""
function invert(cscheme::Vector{C}, c, rangescale :: Tuple{Number, Number}=(0.0, 1.0)) where {C<:Colorant}
    if length(cscheme) <= 1 ; throw(InexactError()) ; end
    cdiffs = map(c_i->colordiff(promote(c,c_i)...), cscheme)
    closest = indmin(cdiffs)
    left = right = 0;
    if closest == 1 ; left = closest; right = closest + 1;
    elseif closest == length(cscheme) ; left = closest - 1; right = closest;
    else
        next_closest = cdiffs[closest-1] < cdiffs[closest+1] ? closest-1 : closest+1
        left = min(closest, next_closest)
        right = max(closest, next_closest)
    end

    v = left
    if cdiffs[left] != cdiffs[right] ;  # Prevents divide by 0.
        v += ( cdiffs[left] / (cdiffs[left] + cdiffs[right]))
     end
    return ColorSchemes.remap(v, 1, length(cscheme), rangescale...)
end

"""
    convertToScheme(cscheme, img)

Converts img from its current color values to use only the colors defined in cscheme.

```julia
image = nonTransparentImg
convertToScheme(ColorSchemes.leonardo, image)
convertToScheme(ColorSchemes.Paired_12, image)
```
"""
convertToScheme(cscheme::Vector{C},img) where {C<:Colorant} =
    map(c->get(cscheme, invert(cscheme, c)), img)

end


text:
abstract type AbstractWavefront end

"""
Represents a wavefront.
"""
struct Wavefront <: AbstractWavefront
    phase::Function
    amplitude::Function
end

"""
Represents the addition of a two wavefronts. It is created when `+` operator is
used on two [`Wavefront`](@ref).
"""
struct AddWavefront <: AbstractWavefront
    left::AbstractWavefront
    right::AbstractWavefront
end
Wavefront(phase::Function, amplitude::Real) = Wavefront(phase, (a,b)->amplitude)
Wavefront(ab::Aberration,coordinates::Symbol=:cartesian) = Wavefront(phase(ab, coordinates), 1)
Wavefront(ab::Aberration, amplitude::Function, coordinates::Symbol=:cartesian) = Wavefront(phase(ab, coordinates), amplitude)
Wavefront() = Wavefront(phase(Piston(0), :cartesian), 1)
Base.:*(num::Number, wavefront::Wavefront) = Wavefront(wavefront.phase, (a,b)->num*wavefront.amplitude(a,b))
Base.:*(wavefront::Wavefront, num::Number) = Wavefront(wavefront.phase, (a,b)->num*wavefront.amplitude(a,b))

Base.:+(wavefront_a::Wavefront, wavefront_b::Wavefront) = AddWavefront(wavefront_a, wavefront_b)

"""
    field(wavefront)

Create a function mapping the field in the coordinates of the wavefront.

See also: [`intensity`](@ref), [`amplitude`](@ref), [`phase`](@ref)
"""
field(wavefront::Wavefront) = (a,b) -> wavefront.amplitude(a,b) * exp(1im*wavefront.phase(a,b))
field(wavefront::AddWavefront) = begin
    field_a = field(wavefront.left)
    field_b = field(wavefront.right)
    (a,b) -> field_a(a,b) + field_b(a,b)
end
field(wavefront::AbstractWavefront, reference::AbstractWavefront) = field(wavefront+reference)

"""
    intensity(wavefront)

Create a function mapping the intensity in the coordinates of the wavefront.

See also: [`field`](@ref), [`amplitude`](@ref), [`phase`](@ref)
"""
intensity(wavefront::AbstractWavefront) = begin
    f = field(wavefront)
    (a,b) -> abs2(f(a,b))
end

"""
    intensity(wavefront, reference)

Create a function mapping the intensity of the interference pattern of the
wavefront with the reference. Both should be in the same coordinate system.

See also: [`field`](@ref), [`amplitude`](@ref), [`phase`](@ref)
"""
intensity(wavefront::AbstractWavefront, reference::AbstractWavefront) = intensity(wavefront+reference)

"""
    phase(wavefront)

Create a function mapping the phase of the wavefront.

See also: [`intensity`](@ref), [`amplitude`](@ref), [`field`](@ref)
"""
phase(wavefront::Wavefront) = wavefront.phase
phase(wavefront::AddWavefront) = begin
    sumf = field(wavefront)
    (a,b) -> angle(sumf(a,b))
end

"""
    amplitude(wavefront)

Create a function mapping the amplitude of the wavefront.

See also: [`intensity`](@ref), [`field`](@ref), [`phase`](@ref)
"""
amplitude(wavefront::Wavefront) = wavefront.amplitude
amplitude(wavefront::AddWavefront) = begin
    sumf = field(wavefront)
    (a,b) -> abs(sumf(a,b))
end


text:
# dependencies
using LFAToolkit
using LinearAlgebra
using CSV
using DataFrames

# setup
maxeigenvalues = DataFrame()

struct RunOption
    p::Int
    ω::Float64
    v::AbstractArray{Int,1}
    dimension::Int
end

runoptions = [
    RunOption(2, 1.000, [0, 1], 1),
    RunOption(2, 1.000, [1, 0], 1),
    RunOption(2, 1.000, [1, 1], 1),
    RunOption(2, 1.000, [1, 2], 1),
    RunOption(2, 1.000, [2, 1], 1),
    RunOption(2, 1.000, [2, 2], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [0, 1], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [1, 0], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [1, 1], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [1, 2], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [2, 1], 1),
    RunOption(2, 14 / (22 - 2sqrt(7)), [2, 2], 1),
    RunOption(2, 56 / 79, [0, 1], 1),
    RunOption(2, 56 / 79, [1, 0], 1),
    RunOption(2, 56 / 79, [1, 1], 1),
    RunOption(2, 56 / 79, [1, 2], 1),
    RunOption(2, 56 / 79, [2, 1], 1),
    RunOption(2, 56 / 79, [2, 2], 1),
    RunOption(3, 0.650, [0, 1], 1),
    RunOption(4, 0.640, [0, 1], 1),
    RunOption(2, 1.000, [0, 1], 2),
    RunOption(2, 1.000, [1, 1], 2),
    RunOption(2, 1.000, [1, 2], 2),
    RunOption(2, 1.000, [2, 1], 2),
    RunOption(2, 1.000, [2, 2], 2),
]

# test range
for runoption in runoptions
    println(
        "p: ",
        runoption.p,
        ", ω: ",
        round(runoption.ω; digits = 3),
        ", v: [",
        runoption.v[1],
        ", ",
        runoption.v[2],
        "], dimension: ",
        runoption.dimension,
    )
    # setup
    mesh = []
    if runoption.dimension == 1
        mesh = Mesh1D(1.0)
    elseif runoption.dimension == 2
        mesh = Mesh2D(1.0, 1.0)
    end
    numbercomponents = 1
    numberfineelements1d = 2
    ctofbasis = TensorH1UniformHProlongationBasis(
        runoption.p + 1,
        numbercomponents,
        runoption.dimension,
        numberfineelements1d,
    )

    # operators
    function diffusionweakform(du::Array{Float64}, w::Array{Float64})
        dv = du * w[1]
        return [dv]
    end
    # -- fine
    basis = TensorH1UniformMacroBasis(
        runoption.p + 1,
        runoption.p + 2,
        numbercomponents,
        runoption.dimension,
        numberfineelements1d,
    )
    inputs = [
        OperatorField(basis, [EvaluationMode.gradient]),
        OperatorField(basis, [EvaluationMode.quadratureweights]),
    ]
    outputs = [OperatorField(basis, [EvaluationMode.gradient])]
    finediffusion = Operator(diffusionweakform, mesh, inputs, outputs)
    # -- coarse
    basis = TensorH1UniformBasis(
        runoption.p + 1,
        runoption.p + 2,
        numbercomponents,
        runoption.dimension,
    )
    inputs = [
        OperatorField(basis, [EvaluationMode.gradient]),
        OperatorField(basis, [EvaluationMode.quadratureweights]),
    ]
    outputs = [OperatorField(basis, [EvaluationMode.gradient])]
    coarsediffusion = Operator(diffusionweakform, mesh, inputs, outputs)

    # -- Jacobi smoother
    jacobi = Jacobi(finediffusion)

    # -- h-multigrid preconditioner
    multigrid = HMultigrid(finediffusion, coarsediffusion, jacobi, [ctofbasis])

    # compute smoothing factor
    # -- setup
    numbersteps = 100
    maxeigenvalue = 0
    θ_min = -π / 2
    θ_max = 3π / 2
    θ_step = 2π / numbersteps
    θ_range = θ_min:θ_step:(θ_max-θ_step)

    θ_maxeigenvalue = -1
    # -- 1D --
    if runoption.dimension == 1
        for i = 1:numbersteps
            θ = [θ_range[i]]
            if abs(θ[1]) > π / 128
                M = computesymbols(multigrid, [runoption.ω], runoption.v, θ)
                eigenvalues = [abs(val) for val in eigvals(M)]
                currentmaxeigenvalue = max(eigenvalues...)
                if (currentmaxeigenvalue > maxeigenvalue)
                    maxeigenvalue = currentmaxeigenvalue
                    θ_maxeigenvalue = θ / π
                end
            end
        end
        # -- 2D --
    elseif runoption.dimension == 2
        for i = 1:numbersteps, j = 1:numbersteps
            θ = [θ_range[i], θ_range[j]]
            if sqrt(abs(θ[1])^2 + abs(θ[2])^2) > π / 128
                M = computesymbols(multigrid, [runoption.ω], runoption.v, θ)
                eigenvalues = [abs(val) for val in eigvals(M)]
                currentmaxeigenvalue = max(eigenvalues...)
                if (currentmaxeigenvalue > maxeigenvalue)
                    maxeigenvalue = currentmaxeigenvalue
                    θ_maxeigenvalue = θ / π
                end
            end
        end
    end
    append!(
        maxeigenvalues,
        DataFrame(
            p = runoption.p,
            ω = runoption.ω,
            θ = θ_maxeigenvalue,
            v = (10 * runoption.v[1] + runoption.v[2]),
            ρ = maxeigenvalue,
        ),
    )
end

CSV.write("two-grid-jacobi-hmultigrid.csv", maxeigenvalues)


text:
function calc_totaladaptationcosts(p, v, d, t, annual_year, r)
    # setting the year for entry in lists.
    yr = annual_year - 2015 + 1 # + 1 because of 1-based indexing in Julia

    v.act_adaptationcosts_total_ann[yr,r] = p.ac_adaptationcosts_economic_ann[yr,r] + p.ac_adaptationcosts_sealevelrise_ann[yr,r] + p.ac_adaptationcosts_noneconomic_ann[yr,r]
    v.act_percap_adaptationcosts_ann[yr,r] = v.act_adaptationcosts_total_ann[yr,r] / p.pop_population_ann[yr,r]
end

@defcomp TotalAdaptationCosts begin
    region = Index()
    year = Index()
    # Total Adaptation Costs
    pop_population = Parameter(index=[time, region], unit="million person")
    ac_adaptationcosts_economic = Parameter(index=[time, region], unit="\$million")
    ac_adaptationcosts_noneconomic = Parameter(index=[time, region], unit="\$million")
    ac_adaptationcosts_sealevelrise = Parameter(index=[time, region], unit="\$million")
    pop_population_ann = Parameter(index=[year, region], unit="million person")
    ac_adaptationcosts_economic_ann = Parameter(index=[year, region], unit="\$million")
    ac_adaptationcosts_noneconomic_ann = Parameter(index=[year, region], unit="\$million")
    ac_adaptationcosts_sealevelrise_ann = Parameter(index=[year, region], unit="\$million")

    act_adaptationcosts_total = Variable(index=[time, region], unit="\$million")
    act_adaptationcosts_total_ann = Variable(index=[year, region], unit="\$million")
    act_percap_adaptationcosts = Variable(index=[time, region], unit="\$/person")
    act_percap_adaptationcosts_ann = Variable(index=[year, region], unit="\$/person")

    function run_timestep(p, v, d, t)

        for r in d.region
            v.act_adaptationcosts_total[t,r] = p.ac_adaptationcosts_economic[t,r] + p.ac_adaptationcosts_sealevelrise[t,r] + p.ac_adaptationcosts_noneconomic[t,r]
            v.act_percap_adaptationcosts[t,r] = v.act_adaptationcosts_total[t,r] / p.pop_population[t,r]

            # calculate  for this specific year
            if is_first(t)
                for annual_year = 2015:(gettime(t))
                    calc_totaladaptationcosts(p, v, d, t, annual_year, r)
                end
            else
                for annual_year = (gettime(t - 1) + 1):(gettime(t))
                    calc_totaladaptationcosts(p, v, d, t, annual_year, r)
                end
            end

        end
    end
end


text:
using Documenter, StateEstimation

makedocs(
    format=:html,
    sitename="StateEstimation",
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "Working with States" => "state_examples.md",
            "Using Systems and Observers" => "system_observer_examples.md",
            "Linear Filtering" => "kalman_filter_examples.md",
            "Nonlinear Filtering" => "nonlinear_kalman_filter_examples.md",
            "Batch Estimation" => "batch_estimation_examples.md",
            #"Multi-Target Tracking" => "multi_target_examples.md"
            ],
        "All Documentation" => Any[
            "Core Components" => Any[
                "States" => "states.md",
                "Systems" => "systems.md",
                "Observers" => "observers.md"
            ],
            "Sequential Estimation" => Any[
                "Kalman Filtering" => "kalman_filter.md",
                "Multi-Target Tracking" => "multi_target.md",
            ],
            "Batch Estimation" => Any[
                "Least Squares" => "least_squares.md",
            ],
            "Simulation" => "simulation.md"
        ]
    ])


text:
# NOTE:  these are several hacks to make the forward declarations work
# this side steps dispatching when making ccalls with Ptr{Abstract} using a Ref{Concrete<:Abstract}
struct _Ref{To, From, R<:Base.RefValue{From}}
	r::R
end
Base.cconvert(::Type{Ptr{CO}}, r::R) where {CO<:Copaques, R<:Base.RefValue{CO}} = r
Base.cconvert(::Type{Ptr{CO}}, r::R) where {CO<:Copaques, T<:CO, R<:Base.RefValue{T}} = _Ref{CO, T, R}(r)
Base.unsafe_convert(::Type{Ptr{CO}}, r::_Ref{CO, T, R}) where {CO<:Copaques, T<:CO, R<:Base.RefValue{T}} = reinterpret(Ptr{CO}, Base.unsafe_convert(Ptr{T}, r.r))
Base.unsafe_load(p::Ptr{CO}, i::Integer = 1) where {CO<:Copaques} = Base.pointerref(reinterpret(Ptr{concrete(CO)}, p), Int(i), 1)
Base.unsafe_store!(p::Ptr{CO}, x, i::Integer = 1) where {CO<:Copaques} = Base.pointerset(reinterpret(Ptr{concrete(CO)}, p), convert(concrete(CO), x), Int(i), 1)

const _alignExprs = (Symbol("@calign"), :(CBinding.$(Symbol("@calign"))))
const _enumExprs = (Symbol("@cenum"), :(CBinding.$(Symbol("@cenum"))))
const _arrayExprs = (Symbol("@carray"), :(CBinding.$(Symbol("@carray"))))
const _structExprs = (Symbol("@cstruct"), :(CBinding.$(Symbol("@cstruct"))))
const _unionExprs = (Symbol("@cunion"), :(CBinding.$(Symbol("@cunion"))))
const _externExprs = (Symbol("@cextern"), :(CBinding.$(Symbol("@cextern"))))

# macros need to accumulate definition of sub-structs/unions and define them above the expansion of the macro itself
_expand(mod::Module, deps::Vector{Pair{Symbol, Expr}}, x, escape::Bool = true) = escape ? esc(x) : x
function _expand(mod::Module, deps::Vector{Pair{Symbol, Expr}}, e::Expr, escape::Bool = true)
	if Base.is_expr(e, :macrocall)
		if length(e.args) > 1 && e.args[1] in (_alignExprs..., _enumExprs..., _arrayExprs..., _structExprs..., _unionExprs...)
			if e.args[1] in _alignExprs
				return _calign(mod, deps, filter(x -> !(x isa LineNumberNode), e.args[2:end])...)
			elseif e.args[1] in _enumExprs
				return _cenum(mod, deps, filter(x -> !(x isa LineNumberNode), e.args[2:end])...)
			elseif e.args[1] in _arrayExprs
				return _carray(mod, deps, filter(x -> !(x isa LineNumberNode), e.args[2:end])...)
			elseif e.args[1] in _structExprs
				return _caggregate(mod, deps, :cstruct, filter(x -> !(x isa LineNumberNode), e.args[2:end])...)
			elseif e.args[1] in _unionExprs
				return _caggregate(mod, deps, :cunion, filter(x -> !(x isa LineNumberNode), e.args[2:end])...)
			end
		else
			return _expand(mod, deps, macroexpand(mod, e, recursive = false))
		end
	elseif Base.is_expr(e, :ref, 2)
		return _carray(mod, deps, e)
	else
		for i in eachindex(e.args)
			e.args[i] = _expand(mod, deps, e.args[i], escape)
		end
		return e
	end
end

macro calign(exprs...) return _calign(__module__, nothing, exprs...) end

function _calign(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, expr::Union{Integer, Expr})
	isOuter = isnothing(deps)
	deps = isOuter ? Pair{Symbol, Expr}[] : deps
	def = Expr(:align, _expand(mod, deps, expr))

	return isOuter ? quote $(map(last, deps)...) ; $(def) end : def
end

macro ctypedef(exprs...) return _ctypedef(__module__, nothing, exprs...) end

function _ctypedef(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, name::Symbol, expr::Union{Symbol, Expr})
	escName = esc(name)

	isOuter = isnothing(deps)
	deps = isOuter ? Pair{Symbol, Expr}[] : deps
	expr = _expand(mod, deps, expr)
	push!(deps, name => quote
		const $(escName) = $(expr)
	end)

	return isOuter ? quote $(map(last, deps)...) ; $(escName) end : escName
end


text:
## Testing the matlab interface
out=BAMF.matlab_DD_FlatBG(data.data,"gauss",1.3f0,θ_start,θ_step,len,mypdf,Int32(burnin),Int32(iterations),xystd,istd,split_std,bndpixels)

# mex interface

args=[MATLAB.mxarray(data.data),
MATLAB.mxarray("gauss"),MATLAB.mxarray(1.3f0),MATLAB.mxarray(θ_start),MATLAB.mxarray(θ_step),
MATLAB.mxarray(len),MATLAB.mxarray(mypdf),MATLAB.mxarray(Int32(burnin)),MATLAB.mxarray(Int32(iterations)),
MATLAB.mxarray(xystd),MATLAB.mxarray(istd),MATLAB.mxarray(split_std),MATLAB.mxarray(bndpixels)
]

BAMF.mextypes(args)
BAMF.mextest(args)

mapn=BAMF.matlab_DD_FlatBG_mex(args)

@time mapn=BAMF.matlab_DD_FlatBG_mex_lite(args);


text:
include("./simulate.jl")

# Compile
params = Params(nsteps=10)
simulate(params)

# 100k simulation
params = Params(
    nsteps=200,
    x0 = getcoor("./ne100k_initial.pdb"),
    cutoff = 12.,
    box = Box([ 99.88, 99.88, 99.88 ], 12., lcell=2)
)
@time simulate(params)


text:
import Base.setindex!, Base.getindex, Base.isapprox, Base.==, Base.!=,
        Base.+, Base.-, Base.*, Base./, Base.write

mutable struct SphericalHarmonicCoefficients
  c::Vector{T} where T<:Real
  L::Int
  R::Float64
  solid::Bool

  function SphericalHarmonicCoefficients(c::Vector{T} where T<:Real)
    if (!isinteger(sqrt(length(c))) || length(c)==0)
      throw(DomainError(sqrt(length(c)),"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
    end

    L = convert(Int,sqrt(length(c)))-1
    return new(c,L,1.0,false)
  end

  function SphericalHarmonicCoefficients(L::Int)
    if L<0
      throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
    end
    return new(zeros((L+1)^2),UInt(L),1.0,false)
  end

  function SphericalHarmonicCoefficients(L::Int,R::Float64,solid::Bool)
    if L<0
      throw(DomainError(L,"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
    end
    return new(zeros((L+1)^2),UInt(L),R,solid)
  end

  function SphericalHarmonicCoefficients(c::Vector{T} where T<:Real,R::Float64,solid::Bool)
    if (!isinteger(sqrt(length(c))) || length(c)==0)
       throw(DomainError(sqrt(length(c)),"Input vector needs to be of size (L+1)², where L ∈ ℕ₀."))
    end

    L = convert(Int,sqrt(length(c)))-1
    return new(c,L,R,solid)
  end

  function SphericalHarmonicCoefficients(c::Array{Vector{T}} where T<:Real,R::Array{Float64},solid::BitArray)
    if size(c) != size(R) || size(c) != size(solid)
       throw(DomainError(size(c),"Arrays do not have the same size."))
    end
    for co in c
	if (!isinteger(sqrt(length(co))) || length(co)==0)
           throw(DomainError(sqrt(length(co)),"Input vectors need to be of size (L+1)², where L ∈ ℕ₀."))
        end
    end

    L = [convert(Int,sqrt(length(co)))-1 for co in c]
    return reshape([new(c[n],L[n],R[n],solid[n]) for n in eachindex(c)],size(c)...)
  end

  function SphericalHarmonicCoefficients(c::Array{Vector{T}} where T<:Real,R::Float64,solid::Bool)
      for co in c
	  if (!isinteger(sqrt(length(co))) || length(co)==0)
             throw(DomainError(sqrt(length(co)),"Input vectors need to be of size (L+1)², where L ∈ ℕ₀."))
          end
      end

    L = [convert(Int,sqrt(length(co)))-1 for co in c]
    return reshape([new(c[n],L[n],R,solid) for n in eachindex(c)],size(c)...)
  end
end

# write and read coefficients to/from an HDF5-file
function SphericalHarmonicCoefficients(path::String)
  file = h5open(path,"r")
      coeffsArray = read(file, "/coeffs")
      R = read(file, "/normalization")
      solid = (read(file, "/solid") .== 1)

  if size(coeffsArray)[1:end-1] == ()
      coeffs = Array{Vector{Float64}}(undef,1)
  else
      coeffs = Array{Vector{Float64}}(undef,size(coeffsArray)[1:end-1])
  end
  coeffsArray = reshape(coeffsArray,(Int(length(coeffsArray)/size(coeffsArray)[end]),size(coeffsArray)[end]))
  for n=1:size(coeffsArray,1)
      coeffs[n] = coeffsArray[n,:]
  end

  return SphericalHarmonicCoefficients(coeffs, R, solid)
end

function write(path::String, coeffs::Array{SphericalHarmonicCoefficients})

  if size(coeffs) != (1,)
      coeffsArray = coeffs[1].c'
      for (n,co) in enumerate(coeffs[2:end])
          coeffsArray = vcat(coeffsArray,co.c')
      end
      coeffsArray = reshape(coeffsArray,(size(coeffs)...,(coeffs[1].L+1)^2))
  else
      coeffsArray = coeffs[1].c
  end

  R = Array{Float64}(undef,size(coeffs))
  solid = Array{Int}(undef,size(coeffs))
  for (n,co) in enumerate(coeffs)
      R[n] = co.R
      solid[n] = co.solid ? 1 : 0
  end

  h5open(path,"w") do file
      write(file, "/coeffs", coeffsArray)
      write(file, "/normalization", R)
      write(file, "/solid", solid)
  end
end

isapprox(shc1::SphericalHarmonicCoefficients, shc2::SphericalHarmonicCoefficients; kargs...) = isapprox(shc1.c,shc2.c;kargs...) && isapprox(shc1.R,shc2.R;kargs...) && isapprox(shc1.solid,shc2.solid;kargs...)
getindex(shc::SphericalHarmonicCoefficients,I) = getindex(shc.c,I)
getindex(shc::SphericalHarmonicCoefficients,l,m) = getindex(shc.c, l*(l+1)+m+1)
setindex!(shc::SphericalHarmonicCoefficients,value,I) = setindex!(shc.c,value,I)
setindex!(shc::SphericalHarmonicCoefficients,value,l,m) = setindex!(shc.c,value,l*(l+1)+m+1)
==(shca::SphericalHarmonicCoefficients, shcb::SphericalHarmonicCoefficients) = shca.c == shcb.c
!=(shca::SphericalHarmonicCoefficients, shcb::SphericalHarmonicCoefficients) = shca.c != shcb.c
+(shc::SphericalHarmonicCoefficients, value) = SphericalHarmonicCoefficients(shc.c .+ value,shc.R,shc.solid)
-(shc::SphericalHarmonicCoefficients, value) = SphericalHarmonicCoefficients(shc.c .- value,shc.R,shc.solid)
*(shc::SphericalHarmonicCoefficients, value) = SphericalHarmonicCoefficients(shc.c .* value,shc.R,shc.solid)
/(shc::SphericalHarmonicCoefficients, value) = SphericalHarmonicCoefficients(shc.c ./ value,shc.R,shc.solid)
+(value, shc::SphericalHarmonicCoefficients) = +(shc::SphericalHarmonicCoefficients, value)
-(value, shc::SphericalHarmonicCoefficients) = SphericalHarmonicCoefficients(value .- shc.c,shc.R,shc.solid)
*(value, shc::SphericalHarmonicCoefficients) = *(shc::SphericalHarmonicCoefficients, value)
LinearAlgebra.normalize(shc::SphericalHarmonicCoefficients,R::Float64) = SphericalHarmonicCoefficients([shc[l,m] * 1/(R^l) for l = 0:shc.L for m = -l:l],shc.R/R,shc.solid)
LinearAlgebra.normalize!(shc::SphericalHarmonicCoefficients,R::Float64) = SphericalHarmonicCoefficients([shc[l,m] *= 1/(R^l) for l = 0:shc.L for m = -l:l],shc.R /= R,shc.solid)

function +(shca::SphericalHarmonicCoefficients, shcb::SphericalHarmonicCoefficients)
    if shca.R != shcb.R
	throw(DomainError(shca.R,"Coefficients do not have the same normalization factor."))
    end
    if shca.solid != shcb.solid
        throw(DomainError(shca.solid,"Coefficients do not have the same type."))
    end
    return SphericalHarmonicCoefficients(shca.c + shcb.c,shca.R,shca.solid)
end
function -(shca::SphericalHarmonicCoefficients, shcb::SphericalHarmonicCoefficients)
    if shca.R != shcb.R
        throw(DomainError(shca.R,"Coefficients do not have the same normalization factor."))
    end
    if shca.solid != shcb.solid
        throw(DomainError(shca.solid,"Coefficients do not have the same type."))
    end
    return SphericalHarmonicCoefficients(shca.c - shcb.c,shca.R,shca.solid)
end

## Convertions ##
# convert spherical coefficients to solid coefficients
function solid!(shc::SphericalHarmonicCoefficients)
    if !(shc.solid)
        for l = 0:shc.L
            for m = -l:l
                shc[l,m] *= sqrt((2*l+1)/(4*pi))
            end
        end
        shc.solid = true
    end
    return shc
end
# convert solid coefficients to spherical coefficients
function spherical!(shc::SphericalHarmonicCoefficients)
    if shc.solid
        for l = 0:shc.L
            for m = -l:l
                shc[l,m] *= sqrt(4*pi/(2*l+1))
            end
        end
        shc.solid = false
    end
    return shc
end

"""
    sphericalHarmonicsExpansion(Clm::SphericalHarmonicCoefficients, x::Variable, y::Variable, z::Variable)
*Description:*  Calculation of the spherical or solid harmonics expansion in Cartesian coordinates
                for given coefficients which define the maximum degree of the spherical harmonics\\
*Input:*  `Clm`       - Array with coefficients (length = (l+1)², l = max. deg. of the spherical harmonics)\\
          `x, y, z`   - Cartesian coordinates\\
*Output:*  Spherical/Solid harmonics expansion
"""
function sphericalHarmonicsExpansion(Clm::SphericalHarmonicCoefficients, x::Variable, y::Variable, z::Variable)

  sum = 0

  for l in 0:Clm.L
    for m in -l:l
      if Clm[l,m] != 0
          if Clm.solid
              # solid expansion
              sum += Clm[l,m] * rlm(l,m,x,y,z)
          else
              # spherical expansion
              sum += Clm[l,m] * rlylm(l,m,x,y,z)
          end
      end
    end
  end
  return sum
end


text:
using Test

@testset "MultiKDE" begin
    include("kde.jl")
end


text:
# add some syntax to manually promote to full dummy coding
fulldummy(t::AbstractTerm) =
    throw(ArgumentError("can't promote $t (of type $(typeof(t))) to full dummy " *
                        "coding (only CategoricalTerms)"))

function StatsModels.apply_schema(
    t::FunctionTerm{typeof(fulldummy)},
    sch::StatsModels.FullRank,
    Mod::Type{<:RegressionModel},
)
    fulldummy(apply_schema.(t.args_parsed, Ref(sch), Mod)...)
end

function fulldummy(t::CategoricalTerm)
    new_contrasts = StatsModels.ContrastsMatrix(
        StatsModels.FullDummyCoding(),
        t.contrasts.levels,
    )
    t = CategoricalTerm(t.sym, new_contrasts)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ fa0ed272-bccc-42a1-b8d5-23015f7a71e2
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.Registry.update()
	Pkg.add("Yao")
	Pkg.add("YaoPlots")
end

# ╔═╡ b978fdae-05ac-11eb-21b9-a584bead4705
using Yao, YaoPlots

# ╔═╡ efa8e4a8-05ab-11eb-2faa-bb86b93de1fb
md"# Assignment Sheet 1"

# ╔═╡ 1002aa7c-05ac-11eb-3b42-87809ac716c4
md" **_Assignment 1_:**
1. _Make the following circuit in Yao_.
2. _Make the qubits with state `` |0000〉 ``, and feed it to the circuit._ Hint: You can use either of ` zero_state(number of qubits) ` function or ` ArrayReg(bit\" \") ` function
3. _Measure the state of all the 4 qubits 1024 times._"

# ╔═╡ 98ceb0e4-05ac-11eb-3058-e59aaaf58e82
begin
	a1_1 = chain(4, put(1=>X), repeat(H, 2:4), control(2, 1=>X), control(4, 3=>X), control(3, 1=>X), control(4, 3=>X), repeat(H, 1:4))
plot(a1_1)
end

# ╔═╡ 25293028-05ad-11eb-1e99-7f91895810fb
begin
	a1circuit = chain(4) #Complete the circuit
	plot(a1circuit)
end

# ╔═╡ 6db30ade-05ae-11eb-2343-e1a2986b6b54
Matrix(a1circuit) == Matrix(a1_1) ? md"✅" : md"❌"

# ╔═╡ 7995ad56-05af-11eb-118a-c3ee4c34d824
fourzerostate = 0              #Remove that 0 and make the state "0000"

# ╔═╡ b816e74c-05b0-11eb-3f8c-09c2e70d3ac9
fourzerostate == zero_state(4) ? md"✅" : md"❌"

# ╔═╡ 16d5e24e-05b1-11eb-2352-8d897cc95f75
a1measuredstate = 0             #Remove the 0 and feed the qubits to the circuit you created, then measure them.

# ╔═╡ ce40470a-05c6-11eb-18b2-2996f097502f
let
	flag1 = sum(a1measuredstate .== bit"0000")
	flag2 = sum(a1measuredstate .== bit"1111")
	(flag1/10.24>38 && flag1/10.24<61) && (flag2/10.24>38 && flag2/10.24<61) && (flag1+flag2==1024) ? md"✅. The number of times 0000 comes as value after qubit's measurement is $(flag1/10.24) and number of times 1111 comes is $(flag2/10.24)" : md"❌"
end

# ╔═╡ 4c971cb6-0665-11eb-17b7-df3aa73ce750
md"**_Assignment 2:_** "

# ╔═╡ 7ec00838-0665-11eb-1d58-193cf8fcc24b
md"Did you know that you can use a circuit in another circuit, in Yao? What does that mean? Well..."

# ╔═╡ cc07db8c-0665-11eb-0d54-5b74156a06fc
#For example
let
	circuit1 = chain(2, put(1=>X), put(2=>Y))
	circuit2 = chain(4, put(1:2 => circuit1), put(3=>Z), put(4=>H))
	plot(circuit2)
end

# ╔═╡ 5cdf4a0c-0666-11eb-2ca9-c36074271ac2
md"1. _Make the Bell Circuit_
2. _Make the Reverse Bell Circuit_
3. _Make a circuit which takes two qubits and passes it through, first the Bell Circuit, and then the Reverse Bell Circuit. Use the two circuits you created above._
4. _Create the qubits_ `` |00〉, |01〉, |10〉 `` and `` |11〉 ``, _and pass them through the circuit you created in step 3. Measure each time you pass them through the circuit, 1024 times_."

# ╔═╡ 7ed55e02-0667-11eb-1570-212652d5fe18
begin
	bellcircuit = chain(2) #complete the bell circuit
	plot(bellcircuit)
end

# ╔═╡ 766552de-067d-11eb-34fb-bbc0fc471071
Matrix(bellcircuit) == Matrix(chain(2, put(1=>H), control(1,2=>X))) ? md"✅" : md"❌"

# ╔═╡ a59c3f42-0667-11eb-34e2-c9ab7e9080f2
begin
	reversebellcircuit = chain(2) #complete the reverse bell circuit
	plot(bellcircuit)
end

# ╔═╡ c2a9bb30-067d-11eb-0594-5d416dc1e763
Matrix(reversebellcircuit) == Matrix(chain(2, control(1,2=>X), put(1=>H))) ? md"✅" : md"❌"

# ╔═╡ b0bfd050-0667-11eb-2ed7-1b6d4ea15816
begin
	bell_and_reverse_bell_circuit = chain(2) #complete the circuit as stated in step 3
	plot(bell_and_reverse_bell_circuit)
end

# ╔═╡ 5d52aca4-067f-11eb-3174-a966da065d76
Matrix(bell_and_reverse_bell_circuit) == Matrix(chain(2, put(1=>H), control(1, 2=>H), control(1, 2=>H), put(1=>H))) ? md"✅" : md"❌"

# ╔═╡ d87106c6-067d-11eb-3beb-8bfea8d168a3
begin
	qubit00 = 0 #Remove the 0. Create qubits with state 00
	qubit01 = 0 #Remove the 0. Create qubits with state 01
	qubit10 = 0 #Remove the 0. Create qubits with state 10
	qubit11 = 0 #Remove the 0. Create qubits with state 11
end

# ╔═╡ d1019e06-067f-11eb-2d40-4d15fd34a149
(qubit00 == ArrayReg(bit"00")) && (qubit01 == ArrayReg(bit"10")) && (qubit10 == ArrayReg(bit"01")) && (qubit11 == ArrayReg(bit"11")) ? md"✅" : md"❌"

# ╔═╡ 1db39250-0681-11eb-30d5-2143291302de
md"Think of what the output will be, after measurement... Is the answer you're thinking of, same as the output you're getting after measurement?"

# ╔═╡ 4ccfcc26-067e-11eb-14c0-0fe93728efc2
begin
	measurement00 = 0 #Remove the 0. Feed the qubits 00, to the circuit created in step 3, and measure them 1000 times
	measurement01 = 0 #Remove the 0. Feed the qubits 01, to the circuit created in step 3, and measure them 1000 times
	measurement10 = 0 #Remove the 0. Feed the qubits 10, to the circuit created in step 3, and measure them 1000 times
	measurement11 = 0 #Remove the 0. Feed the qubits 11, to the circuit created in step 3, and measure them 1000 times
end

# ╔═╡ 479143ca-0680-11eb-13f2-6b4d089f2f37
(sum(measurement00 .== bit"00") + sum(measurement01 .== bit"10") + sum(measurement10 .== bit"01")  +  sum(measurement11 .== bit"11")) == 4000 ? md"✅" : md"❌"

# ╔═╡ c6efaa82-067e-11eb-0ad1-2518b517ea4e
md"_**Note:**_ If you're using ` ArrayReg(bit\" \") ` to create the qubits, the qubits are read by the circuit in reverse order. For example, if you enter ` ArrayReg(\"1101\") `, the circuit will take it as ` 1011 `. Also, the value of measurement is read in reverse order."

# ╔═╡ ff7b74ba-0680-11eb-2688-3fb96fc978c5
md"_**Assignment 3:**_"

# ╔═╡ 131c831a-0681-11eb-2bab-3d4cf5b23006
md"Suppose that Alice and Bob have 2 pairs of entangled qubits. Both the pairs are in the state $ \frac{|00> \;+\; |11>}{\sqrt2} $. Suppose Alice has an extra qubit, with completely random state."

# ╔═╡ c89d84b4-0681-11eb-3c22-4def72fc3477
md"1. _Make the circuit for quantum teleportation_
2. _Pass the first pair of Alice and Bob's entangled qubit, with Alice's extra qubit, into the quantum teleportation circuit. Then collapse the first two qubit's states by measuring them._ Hint: Use the ` measure!(RemoveMeasured(), input_qubits, [location_of_qubits]) ` function for this.
3. _Make the circuit for Superdense coding_
4. _Use the other pair of entangled qubits and superdense coding to convey the information Alice got from measuring her qubits from quantum teleportation circuit, to Bob._
5. _Use the above information to make Bob's qubit from first pair's state, to jump to Alice's extra qubit state._"

# ╔═╡ e59122aa-0687-11eb-02b0-a7f72525c5e7
begin
	Alice_and_Bobs_first_entangled_qubit = ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize!
	Alice_and_Bobs_second_entangled_qubit = ArrayReg(bit"00") + ArrayReg(bit"11") |> normalize!
	Alices_extra_qubit = rand_state(1) |> normalize!
	state(Alices_extra_qubit)
end

# ╔═╡ b6670002-0688-11eb-3ecc-03476ce15431
begin
	quantumteleportationcircuit = chain(3) #Complete the circuit
	plot(quantumteleportationcircuit)
end

# ╔═╡ 458cf840-06e3-11eb-19fb-a5a4b0c70f0b
Matrix(quantumteleportationcircuit) == Matrix(chain(3, put(1=>H), control(1, 2=>X))) ? md"✅" : md"❌"

# ╔═╡ fc6e0da0-068c-11eb-1c31-e92bc7455530
begin
	input_to_teleportation_circuit = join(Alice_and_Bobs_first_entangled_qubit, Alices_extra_qubit)
	state(input_to_teleportation_circuit)
end

# ╔═╡ 88f1d456-06df-11eb-1827-912c481bb7a4
input = 0
#Remove the "0". Pass the above created qubit through the teleportation circuit, measuring the first two qubits using measure_remove!() function. Hint: qubit |> circuit |> r->measure_remove!(r, m:n)

# ╔═╡ 79280796-06e3-11eb-1e5c-b5cea995a88a
(input == bit"00" || input == bit"11" || input==bit"10" || input==bit"01") && (typeof(input) != Int64) ? md"✅" : md"❌"

# ╔═╡ e933ab46-06df-11eb-312a-6f2e934b059d
md"The above measurement will act as the information for superdense coding circuit."

# ╔═╡ 26cc4008-06e0-11eb-0cdf-51954af7254d
begin
	if(input==bit"00")
		superdense_coding_circuit = chain(2)
	elseif(input==bit"01")
		superdense_coding_circuit = chain(2)
	elseif(input==bit"10")
		superdense_coding_circuit = chain(2)
	elseif(input==bit"11")
		superdense_coding_circuit = chain(2)
	end
	#Remember the bits are read from right to left
	#Complete the circuits for Alice's qubit
	plot(superdense_coding_circuit)
end

# ╔═╡ dc0b86be-06e4-11eb-24db-4f56cb821e0f
begin
	if(input==bit"00")
		Matrix(superdense_coding_circuit) == Matrix(chain(2))
	elseif(input==bit"01")
		Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Z)))
	elseif(input==bit"10")
		Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>X)))
	elseif(input==bit"11")
		Matrix(superdense_coding_circuit) == Matrix(chain(2, put(1=>Y)))
	end && typeof(input) != Int64 ? md"✅" : md"❌"
end

# ╔═╡ 28ccce76-06e1-11eb-1c70-b3d7779a062e
Bobs_part = ((Alice_and_Bobs_second_entangled_qubit |> superdense_coding_circuit) |> reversebellcircuit) |> r->measure(r, nshots=1000)

# ╔═╡ 2aaa0910-06e2-11eb-1d6c-0fae35763d48
#Bobs_part now contains the information Alice wanted to convey to Bob regarding his qubits. Use this information to make Bob's qubit's state jump to Alice's random qubit's state
begin
	if(Bobs_part[1]==bit"00")
		Bobs_qubit = input_to_teleportation_circuit |> chain(1)
	elseif(Bobs_part[1]==bit"01")
		Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit
	elseif(Bobs_part[1]==bit"10")
		Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit
	elseif(Bobs_part[1]==bit"11")
		Bobs_qubit = input_to_teleportation_circuit |> chain(1) #Complete the circuit
	end
	state(Bobs_qubit)
end

# ╔═╡ b1392da2-06e5-11eb-3632-7d96a1b81a24
begin
	log(2,length(state(input_to_teleportation_circuit))) == 2 && ((Bobs_part[1]==bit"00" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1))) || (Bobs_part[1]==bit"01" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>Z)))) || ( Bobs_part[1]==bit"10" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>X))) ) || ( Bobs_part[1]==bit"11" && Bobs_qubit == (input_to_teleportation_circuit |> chain(1, put(1=>Y)))))  &&  (sum(Bobs_part .== Bobs_part[1]) == 1000) ? md"✅" : md"❌"
	#=elseif()
		Bobs_qubit ==
	elseif()
		Bobs_qubit ==
	elseif()
		Bobs_qubit == =#

end

# ╔═╡ 83df1840-06e7-11eb-2502-173dac5d4963
md"Seems to work... Although, wouldn't using the information from superdense coding make teleportation pointless."

# ╔═╡ Cell order:
# ╠═fa0ed272-bccc-42a1-b8d5-23015f7a71e2
# ╟─efa8e4a8-05ab-11eb-2faa-bb86b93de1fb
# ╠═b978fdae-05ac-11eb-21b9-a584bead4705
# ╟─1002aa7c-05ac-11eb-3b42-87809ac716c4
# ╟─98ceb0e4-05ac-11eb-3058-e59aaaf58e82
# ╠═25293028-05ad-11eb-1e99-7f91895810fb
# ╟─6db30ade-05ae-11eb-2343-e1a2986b6b54
# ╠═7995ad56-05af-11eb-118a-c3ee4c34d824
# ╟─b816e74c-05b0-11eb-3f8c-09c2e70d3ac9
# ╠═16d5e24e-05b1-11eb-2352-8d897cc95f75
# ╟─ce40470a-05c6-11eb-18b2-2996f097502f
# ╟─4c971cb6-0665-11eb-17b7-df3aa73ce750
# ╟─7ec00838-0665-11eb-1d58-193cf8fcc24b
# ╠═cc07db8c-0665-11eb-0d54-5b74156a06fc
# ╟─5cdf4a0c-0666-11eb-2ca9-c36074271ac2
# ╠═7ed55e02-0667-11eb-1570-212652d5fe18
# ╟─766552de-067d-11eb-34fb-bbc0fc471071
# ╠═a59c3f42-0667-11eb-34e2-c9ab7e9080f2
# ╟─c2a9bb30-067d-11eb-0594-5d416dc1e763
# ╠═b0bfd050-0667-11eb-2ed7-1b6d4ea15816
# ╟─5d52aca4-067f-11eb-3174-a966da065d76
# ╠═d87106c6-067d-11eb-3beb-8bfea8d168a3
# ╟─d1019e06-067f-11eb-2d40-4d15fd34a149
# ╟─1db39250-0681-11eb-30d5-2143291302de
# ╠═4ccfcc26-067e-11eb-14c0-0fe93728efc2
# ╟─479143ca-0680-11eb-13f2-6b4d089f2f37
# ╟─c6efaa82-067e-11eb-0ad1-2518b517ea4e
# ╟─ff7b74ba-0680-11eb-2688-3fb96fc978c5
# ╟─131c831a-0681-11eb-2bab-3d4cf5b23006
# ╟─c89d84b4-0681-11eb-3c22-4def72fc3477
# ╠═e59122aa-0687-11eb-02b0-a7f72525c5e7
# ╠═b6670002-0688-11eb-3ecc-03476ce15431
# ╟─458cf840-06e3-11eb-19fb-a5a4b0c70f0b
# ╠═fc6e0da0-068c-11eb-1c31-e92bc7455530
# ╠═88f1d456-06df-11eb-1827-912c481bb7a4
# ╟─79280796-06e3-11eb-1e5c-b5cea995a88a
# ╟─e933ab46-06df-11eb-312a-6f2e934b059d
# ╠═26cc4008-06e0-11eb-0cdf-51954af7254d
# ╟─dc0b86be-06e4-11eb-24db-4f56cb821e0f
# ╠═28ccce76-06e1-11eb-1c70-b3d7779a062e
# ╠═2aaa0910-06e2-11eb-1d6c-0fae35763d48
# ╟─b1392da2-06e5-11eb-3632-7d96a1b81a24
# ╟─83df1840-06e7-11eb-2502-173dac5d4963


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    StructuredGrid(X, Y, Z, ...)

A structured grid with coordinates `X`, `Y`, `Z`, ...

## Examples

A 2D structured grid can be constructed from coordinates stored in
2D matrices. For example, we can georeference locations in the Earth
surface using `LAT` and `LON` coordinates:

```julia
julia> StructuredGrid(LAT, LON)
```
"""
struct StructuredGrid{T,N} <: SpatialDomain{T,N}
  coords::Matrix{T}
  dims::Dims{N}
end

function StructuredGrid(coordarrays::Vararg{<:AbstractArray{T},N}) where {N,T}
  sizes = [size(array) for array in coordarrays]
  @assert length(unique(sizes)) == 1 "coordinates arrays must have the same dimensions"

  coords = Matrix{T}(undef, N, prod(sizes[1]))
  for (i, array) in enumerate(coordarrays)
    coords[i,:] = array
  end

  StructuredGrid{T,N}(coords, sizes[1])
end

Base.size(grid::StructuredGrid) = grid.dims

nelms(grid::StructuredGrid) = prod(grid.dims)

function coordinates!(buff::AbstractVector{T}, grid::StructuredGrid{T,N},
                      location::Int) where {N,T}
  buff .= grid.coords[:,location]
end

==(grid₁::StructuredGrid, grid₂::StructuredGrid) =
  grid₁.coords == grid₂.coords && grid₁.dims == grid₂.dims

# ------------
# IO methods
# ------------
function Base.show(io::IO, grid::StructuredGrid{T,N}) where {N,T}
  dims = join(grid.dims, "×")
  print(io, "$dims StructuredGrid{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", grid::StructuredGrid{T,N}) where {N,T}
  println(io, grid)
  Base.print_array(io, grid.coords)
end


text:
using Test, DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus
using LinearAlgebra: norm

@testset "Mesh obtuse triangle" begin
    m = Metric(2)
    s = Simplex(Point(0,0), Point(2,0), Point(1,.5))
    for center in [centroid, circumcenter(m)]
        primal = TriangulatedComplex([s])
        mesh = Mesh(primal, center)
        @test map(length, mesh.primal.complex.cells) ==
            reverse(map(length, mesh.dual.complex.cells))
        cell_1 = mesh.dual.complex.cells[1][1]
        @test cell_1.K == 1
        @test cell_1.points[1] == Point(center(s))
        @test length(cell_1.parents) == 3
        @test length(cell_1.children) == 0
        cells_2 = mesh.dual.complex.cells[2]
        @test all([length(c.parents) == 2 for c in cells_2])
        @test all([length(c.children) == 1 for c in cells_2])
        cells_3 = mesh.dual.complex.cells[3]
        @test all([length(c.parents) == 0 for c in cells_3])
        @test all([length(c.children) == 2 for c in cells_3])
    end
end

@testset "signed_volume" begin
    points = [Point(0, 0), Point(2, 0), Point(1, .5), Point(1, -5)]
    obtuse = Simplex(points[1:3])
    acute = Simplex(points[1], points[2], points[4])
    m = Metric(2)
    b, _ = circumsphere_barycentric(m, obtuse)
    @test any(b.coords .< 0) # obtuse
    b, _ = circumsphere_barycentric(m, acute)
    @test !any(b.coords .< 0) # acute
    tcomp = TriangulatedComplex([acute, obtuse])
    center = circumcenter(m)
    mesh = Mesh(tcomp, center)
    # check that the length of the edge between the cells is correct
    edge = subcomplex(mesh.primal.complex, [Simplex(points[1:2])]).cells[2][1]
    dual_edge = dual(mesh, edge)
    ccs = [Point(center(x)) for x in [Simplex(edge), obtuse, acute]]
    d1 = norm(m, ccs[1].coords - ccs[2].coords)
    d2 = norm(m, ccs[1].coords - ccs[3].coords)
    @test abs(volume(m, mesh.dual, dual(mesh, edge))) == abs(d1 - d2)
    # check that the choice of sign was correct
    simplices, bools = zip(mesh.dual.simplices[dual_edge]...)
    @test bools[1] != bools[2]
    @test Set(simplices[bools[1] ? 1 : 2].points) == Set([ccs[1], ccs[3]])
    @test Set(simplices[bools[1] ? 2 : 1].points) == Set([ccs[1], ccs[2]])
end

@testset "honeycomb mesh" begin
    n = 5
    _, tcomp = DEC.triangulated_lattice(n * [1,0], n * [.5, .5 * sqrt(3)], n, n)
    m = Metric(2)
    mesh = Mesh(tcomp, centroid)
    for c in mesh.primal.complex.cells[1]
        @test volume(m, mesh.primal, c) == 1
    end
    for c in mesh.primal.complex.cells[2]
        @test volume(m, mesh.primal, c) ≈ 1
    end
    for c in mesh.primal.complex.cells[3]
        @test volume(m, mesh.primal, c) ≈ sqrt(3)/4
    end
    for c in mesh.dual.complex.cells[1]
        @test volume(m, mesh.dual, c) == 1
    end
    for c in mesh.dual.complex.cells[2]
        @test volume(m, mesh.dual, c) ≈ sqrt(3)/6 * length(c.children)
    end
    for c in mesh.dual.complex.cells[3]
        k = length(c.children)
        @test volume(m, mesh.dual, c) ≈ (sqrt(3)/12) * (k == 6 ? k : k-1)
    end
end