text:
"""
    nowts()

the current timestamp in seconds from epoch
"""
nowts = () -> Libc.TimeVal().sec

#
# Return the universe index of `ts`: an integer between 1 and `ages`.
#
# `interval` is the correlation interval and `ages` is the number
# of universes.
#
mvindex(ts, interval, ages) = div(ts,interval) % ages + 1

"""
    lastspan(interval::Int64=900)::Int64

Returns the epoch start time of the last span.
The last span is the nearest in the present time interval that satisfies the
condition: `span.endtime < now`.
"""
function lastspan(interval::Int64=900)::Int64
    current = Int(round(tim()))
    current - current%interval - interval
end

"""
    span(ts::Int64, interval::Int64=900)

Returns the start_time of the time interval that contains `ts`.
"""
function span(ts::Int64, interval::Int64=900)
    ts - ts%interval
end


text:
function entropy(model::EoSModel, p, T,   z=SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_entropy(model,v,T,z)
end

function chemical_potential(model::EoSModel, p, T, z= SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_chemical_potential(modle,v,T,z)
end

function internal_energy(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_internal_energy(model,v,T,z)
end

function enthalpy(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_enthalpy(model,v,T,z)
end

function gibbs_free_energy(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_gibbs_free_energy(model,v,T,z)
end

function helmholtz_free_energy(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v      = volume(model, p, T, z; phase=phase)
    return vt_helmholtz_free_energy(model,v,T,z)
end

function isochoric_heat_capacity(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v       = volume(model, p, T, z; phase=phase)
    return vt_isochoric_heat_capacity(model,v,T,z)
end

function isobaric_heat_capacity(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v       = volume(model, p, T, z; phase=phase)
    return vt_isobaric_heat_capacity(model,v,T,z)
end

function isothermal_compressibility(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v       = volume(model, p, T, z; phase=phase)
    return vt_isothermal_compressibility(model,v,T,z)
end

function isentropic_compressibility(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v       = volume(model, p, T, z; phase=phase)
    return vt_isobaric_expansivity(model,v,T,z)
end

function speed_of_sound(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v       = volume(model, p, T, z; phase=phase)
    return vt_speed_of_sound(model,v,T,z)
end

function isobaric_expansivity(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v = volume(model, p, T, z; phase=phase)
    return vt_isobaric_expansivity(model,v,T,z)
end

function joule_thomson_coefficient(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v  = volume(model, p, T, z; phase=phase)
    return vt_joule_thomson_coefficient(model,v,T,z)
end

function compressibility_factor(model::EoSModel, p, T,  z=SA[1.]; phase = :unknown)
    v  = volume(model, p, T, z; phase=phase)
    return p*v/(R̄*T)
end

function inversion_temperature(model::EoSModel,p,z=SA[1.0])
    T0 = 6.75*T_scale(model,z)
    μⱼₜ(T) = joule_thomson_coefficient(model,p,T,z)
    return Roots.find_zero(μⱼₜ,T0)
end

export entropy, chemical_potential, internal_energy, enthalpy, gibbs_free_energy
export helmholtz_free_energy, isochoric_heat_capacity, isobaric_heat_capacity
export isothermal_compressibility, isentropic_compressibility, speed_of_sound
export isobaric_expansivity, joule_thomson_coefficient, compressibility_factor, inversion_temperature

ext: jl


text:
# User settings.
# State grids.
    grids = [   range(0.0, 5.0, length = 500), # none of these are front-loaded to the degree we use in production, which means numerical quantities will differ a bit from algebraic
                range(0.0, 7.0, length = 700),
                unique([range(0.0, 5.0, length = 500)' range(5.0, 7.0, length = 200)']), # irregular
            ]
# Overall parameters.
    params = @with_kw (ρ = 0.02, σ = 4.2508, N = 10, θ = 5.1269, γ = 1.01, d = 2.3701, κ = 0.013, ζ = 1, η = 0, Theta = 1, χ = 1/(2.1868), υ = 0.0593, μ = 0, δ = 0.053)
    baseline = params()
# Solver settings.
    initial_values = [  [0.25, 3.0, 1.0],
                        [0.0190, 1.434969, 1.06517] # ~ equilibrium value for these parameters.
                     ]

# Execute tests.
@testset "Tests with Vanilla Parameters" begin
    # Compute algebraic solutions.
    algebraic_sols = [stationary_algebraic(baseline, settings_defaults_tests(z_ex = grid, stationary_x0 = (x, y) -> iv)) for grid in grids, iv in initial_values]
    # Test algebraic equilibrium quantities.
    # Growth rate tests.
    algebraic_gs = (x -> x.g).(algebraic_sols)
    @test var(algebraic_gs) < 1e-6 # Tests that the solutions are similar to one another.
    @test algebraic_gs[1] ≈ 0.01900455065415125 # Tests proximity to true value.
    # Ω tests.
    algebraic_Ωs = (x -> x.Ω).(algebraic_sols)
    @test var(algebraic_Ωs) < 1e-6
    @test algebraic_Ωs[1] ≈ 1.0651775565541244
    # z_hat tests.
    algebraic_zs = (x -> x.z_hat).(algebraic_sols)
    @test var(algebraic_zs) < 1e-6
    @test all(algebraic_zs .≈ 1.434969541725385)

    # Compute numerical solutions.
    numerical_sols = [stationary_numerical(baseline, settings_defaults_tests(z_ex = grid)) for grid in grids]
    # Test numerical equilibrium quantities.
    numerical_gs = (x -> x.g).(numerical_sols)
    @test var(numerical_gs) < 1e-6 # Tests that the solutions are similar to one another.
    # Ω tests.
    numerical_Ωs = (x -> x.Ω).(numerical_sols)
    @test var(numerical_Ωs) < 1e-6
    # z_hat tests.
    numerical_zs = (x -> x.z_hat).(numerical_sols)
    @test var(numerical_zs) < 1e-6
    # Numerical residuals tests.
    for i in 1:length(numerical_sols)
        # Get values to test.
        grid = grids[i]
        sol = numerical_sols[i]
        # Compute interim quantities.
        ω = ω_weights(grid, baseline.θ, baseline.σ-1)
        value_matching = sol.v_tilde[1] - dot(sol.v_tilde, ω) + sol.x
        free_entry = sol.v_tilde[1] - sol.x*(1-baseline.χ)/baseline.χ
        # Test.
        @test value_matching < 1e-8
        @test free_entry < 1e-8
        @test sol.L_tilde ≈ sol.L_tilde_a + sol.L_tilde_E + sol.L_tilde_x
    end
end

@testset "Tests with Perturbed Parameters" begin
    # Generate new parameters.
    newparams = params(σ = 4.0) # Arbitrary change.
    # Compute algebraic solutions.
    algebraic_sols = [stationary_algebraic(newparams, settings_defaults_tests(z_ex = grid, stationary_x0 = (x, y) -> iv)) for grid in grids, iv in initial_values]
    # Test consistency.
    algebraic_gs = (x -> x.g).(algebraic_sols)
    algebraic_Ωs = (x -> x.Ω).(algebraic_sols)
    algebraic_zs = (x -> x.z_hat).(algebraic_sols)
    @test var(algebraic_gs) < 1e-6 # Tests that the solutions are similar to one another.
    @test var(algebraic_Ωs) < 1e-6
    @test var(algebraic_zs) < 1e-6

    # Compute numerical solutions.
    numerical_sols = [stationary_numerical(newparams, settings_defaults_tests(z_ex = grid)) for grid in grids]
    # Test consistency.
    numerical_gs = (x -> x.g).(numerical_sols)
    numerical_Ωs = (x -> x.Ω).(numerical_sols)
    numerical_zs = (x -> x.z_hat).(numerical_sols)
    @test var(numerical_gs) < 1e-6 # Tests that the solutions are similar to one another.
    @test var(numerical_Ωs) < 1e-6
    @test var(numerical_zs) < 1e-6

    # Numerical residuals tests.
    for i in 1:length(numerical_sols)
        # Get values to test.
        grid = grids[i]
        sol = numerical_sols[i]
        # Compute interim quantities.
        ω = ω_weights(grid, baseline.θ, baseline.σ-1)
        value_matching = sol.v_tilde[1] - dot(sol.v_tilde, ω) + sol.x
        free_entry = sol.v_tilde[1] - sol.x*(1-baseline.χ)/baseline.χ
        # Test.
        @test value_matching < 1e-8
        @test free_entry < 1e-8
        @test sol.L_tilde ≈ sol.L_tilde_a + sol.L_tilde_E + sol.L_tilde_x
    end
end

# === test utilities === #
p = load_parameters("full/calibration_params.csv") # because this is included in runtests.jl, full path required
p = merge(p, (d = p.d_0,))
@test stationary_algebraic(p, settings_defaults()) isa NamedTuple # just needs to run


text:
# load functions and packages

using JuMP, Ipopt, Gurobi, Combinatorics, JLD, HDF5, DelimitedFiles, Statistics, LinearAlgebra;

include("def.jl");
include("readin.jl");
include("main.jl");
include("detForm.jl");
include("auxiliary.jl");
include("forwardPass.jl");
include("backwardPass.jl");


text:
add_verbose_scope(:qAdic)
add_assert_scope(:qAdic)

export setprecision!, defining_polynomial

function Base.setprecision(q::qadic, N::Int)
  r = parent(q)()
  r.N = N
  ccall((:padic_poly_set, :libflint), Nothing, (Ref{qadic}, Ref{qadic}, Ref{FlintQadicField}), r, q, parent(q))
  return r
end

function Base.setprecision(q::padic, N::Int)
  r = parent(q)()
  r.N = N
  ccall((:padic_set, :libflint), Nothing, (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), r, q, parent(q))
  return r
end

function setprecision!(q::qadic, N::Int)
  @assert N >= q.N
  q.N = N
  return q
end

function setprecision!(Q::FlintQadicField, n::Int)
  Q.prec_max = n
end

function setprecision!(Q::FlintPadicField, n::Int)
  Q.prec_max = n
end

function setprecision!(f::Generic.Poly{qadic}, N::Int)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end

function Base.setprecision(f::Generic.Poly{qadic}, N::Int)
  f = deepcopy(f)
  for i=1:length(f)
    f.coeffs[i].N = N
  end
  return f
end

function setprecision!(a::AbstractArray{qadic}, N::Int)
  for x = a
    setprecision!(x, N)
  end
end

function Base.setprecision(a::AbstractArray{qadic}, N::Int)
  return map(x->setprecision(x, N), a)
end

function setprecision!(a::Generic.MatSpaceElem{qadic}, N::Int)
  setprecision!(a.entries, N)
end

function Base.setprecision(a::Generic.MatSpaceElem{qadic}, N::Int)
  b = deepcopy(a)
  setprecision!(b, N)
  return B
end

function trace(r::qadic)
  t = base_ring(parent(r))()
  ccall((:qadic_trace, :libflint), Nothing, (Ref{padic}, Ref{qadic}, Ref{FlintQadicField}), t, r, parent(r))
  return t
end

function norm(r::qadic)
  t = base_ring(parent(r))()
  ccall((:qadic_norm, :libflint), Nothing, (Ref{padic}, Ref{qadic}, Ref{FlintQadicField}), t, r, parent(r))
  return t
end

function setcoeff!(x::fq_nmod, n::Int, u::UInt)
  ccall((:nmod_poly_set_coeff_ui, :libflint), Nothing,
                (Ref{fq_nmod}, Int, UInt), x, n, u)
end

function coeff(x::qadic, i::Int)
  R = FlintPadicField(prime(parent(x)), parent(x).prec_max)
  c = R()
  ccall((:padic_poly_get_coeff_padic, :libflint), Nothing,
           (Ref{padic}, Ref{qadic}, Int, Ref{FlintQadicField}), c, x, i, parent(x))
  return c
end

function setcoeff!(x::qadic, i::Int, y::padic)
  ccall((:padic_poly_set_coeff_padic, :libflint), Nothing,
           (Ref{qadic}, Int, Ref{padic}, Ref{FlintQadicField}), x, i, y, parent(x))
end

function setcoeff!(x::qadic, i::Int, y::UInt)
  R = FlintPadicField(prime(parent(x)), parent(x).prec_max)
  Y = R(fmpz(y))
  ccall((:padic_poly_set_coeff_padic, :libflint), Nothing,
           (Ref{qadic}, Int, Ref{padic}, Ref{FlintQadicField}), x, i, Y, parent(x))
end

function ResidueField(Q::FlintQadicField)
  k = GF(Int(prime(Q)), degree(Q))[1]
  pro = function(x::qadic)
    v = valuation(x)
    v < 0 && error("elt non integral")
    v > 0 && return k(0)
    z = k()
    for i=0:degree(Q)
      setcoeff!(z, i, UInt(lift(coeff(x, i))%prime(Q)))
    end
    return z
  end
  lif = function(x::fq_nmod)
    z = Q()
    for i=0:degree(Q)-1
      setcoeff!(z, i, coeff(x, i))
    end
    return z
  end
  return k, MapFromFunc(pro, lif, Q, k)
end

function ResidueField(Q::FlintPadicField)
  k = GF(Int(prime(Q)))
  pro = function(x::padic)
    v = valuation(x)
    v < 0 && error("elt non integral")
    v > 0 && return k(0)
    z = k(lift(x))
    return z
  end
  lif = function(x::gfp_elem)
    z = Q(lift(x))
    return z
  end
  return k, MapFromFunc(pro, lif, Q, k)
end

function coefficient_ring(Q::FlintQadicField)
  return FlintPadicField(prime(Q), precision(Q))
end
coefficient_field(Q::FlintQadicField) = coefficient_ring(Q)

function prime(R::PadicField, i::Int)
  p = fmpz()
  ccall((:padic_ctx_pow_ui, :libflint), Cvoid, (Ref{fmpz}, Int, Ref{PadicField}), p, i, R)
  return p
end

function getUnit(a::padic)
  u = fmpz()
  ccall((:fmpz_set, :libflint), Cvoid, (Ref{fmpz}, Ref{Int}), u, a.u)
  return u, a.v, a.N
end

function lift_reco(::FlintRationalField, a::padic; reco::Bool = false)
  if reco
    u, v, N = getUnit(a)
    R = parent(a)
    fl, c, d = rational_reconstruction(u, prime(R, N-v))
    !fl && return nothing

    x = FlintQQ(c, d)
    if v < 0
      return x//prime(R, -v)
    else
      return x*prime(R, v)
    end
  else
    return lift(FlintQQ, a)
  end
end

function *(A::fmpz_mat, B::MatElem{padic})
  return matrix(base_ring(B), A) * B
end

uniformizer(Q::FlintQadicField) = Q(prime(Q))
Base.precision(Q::FlintQadicField) = Q.prec_max

uniformizer(Q::FlintPadicField) = Q(prime(Q))
Base.precision(Q::FlintPadicField) = Q.prec_max

nrows(A::Array{T, 2}) where {T} = size(A)[1]
ncols(A::Array{T, 2}) where {T} = size(A)[2]

import Base.^
^(a::qadic, b::qadic) = exp(b*log(a))
^(a::padic, b::padic) = exp(b*log(a))

import Base.//
//(a::qadic, b::qadic) = divexact(a, b)
//(a::padic, b::qadic) = divexact(a, b)
//(a::qadic, b::padic) = divexact(a, b)

function defining_polynomial(Q::FlintQadicField, P::Ring = coefficient_ring(Q))
  Pt, t = PolynomialRing(P, cached = false)
  f = Pt()
  for i=0:Q.len-1
    j = unsafe_load(reinterpret(Ptr{Int}, Q.j), i+1)
    a = fmpz()
    ccall((:fmpz_set, :libflint), Nothing, (Ref{fmpz}, Int64), a, Q.a+i*sizeof(Ptr))
    setcoeff!(f, j, P(a))
  end
  return f
end

function defining_polynomial(Q::FqNmodFiniteField, P::Ring = GF(characteristic(Q)))
  Pt, t = PolynomialRing(P, cached = false)
  f = Pt()
  for i=0:Q.len-1
    j = unsafe_load(reinterpret(Ptr{Int}, Q.j), i+1)
    a = fmpz()
    ccall((:fmpz_set, :libflint), Nothing, (Ref{fmpz}, Int64), a, Q.a+i*sizeof(Ptr))
    setcoeff!(f, j, P(a))
  end
  return f
end


text:
module CausticsEngineering

using DocStringExtensions
using Random, Images

using Plots
gr()

using Meshes, FileIO, MeshIO

include("../examples/personal/original.jl")

include("parameters.jl")
include("utilities.jl")
include("math_utilities.jl")
include("mesh.jl")
include("io.jl")
include("plots.jl")

include("create_mesh.jl")

"""
$(SIGNATURES)
"""
function main()
    @assert size(ARGS) == (1,) "Intented usage is: julia create_mesh.jl image.png"

    img = Images.load(ARGS[1])
    return engineer_caustics(img)
end

export
    # Constants
    Meters_Per_Pixel,
    Bottom_Offset,
    Top_Offset,

    # Types
    Vertex3D,
    FieldVertex3D,
    FaceMesh,

    # math
    average,
    average_absolute,
    laplacian,
    ∇,

    # Mesh
    create_solid,
    create_solid_as_dict,
    get_lens_pixels_area,

    # Plotting
    plot_as_quiver,

    # I/O
    save_obj!,

    # Main procedures
    propagate_poisson,
    march_mesh,
    engineer_caustics,
    main

end # End module


text:
module TensorCalculus

export Tensor, δ, ϵ, ⊗, ⋅, contract, trace, pushover

include("tensor.jl")
include("broadcast.jl")
include("arithmetic.jl")
include("algebra.jl")

end


text:
using BinaryBuilder

# Collection of sources required to build p7zip
sources = [
    "https://downloads.sourceforge.net/project/p7zip/p7zip/16.02/p7zip_16.02_src_all.tar.bz2?r=https%3A%2F%2Fsourceforge.net%2Fprojects%2Fp7zip%2Ffiles%2Fp7zip%2F16.02%2Fp7zip_16.02_src_all.tar.bz2%2Fdownload&ts=1524045807" =>
    "5eb20ac0e2944f6cb9c2d51dd6c4518941c185347d4089ea89087ffdd6e2341f",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
tar xjvf p7zip*
cd p7zip_16.02
sed -i 's/\usr\/local/$WORKSPACE\/destdir/g' install.sh
mv makefile.linux_amd64 makefile.linux
make 7z
./install.sh
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Linux(:powerpc64le, :glibc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "7z", :lib7z),
    LibraryProduct(prefix, "Rar", :librar),
    ExecutableProduct(prefix, "", :p7zip)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "p7zip", sources, script, platforms, products, dependencies)


text:
"""
   shi(x)

   hyperbolic sine integral function

   shi(x) = ₀∫ˣ (sinh(t) / t) δt
""" shi

"""
    chi(x)

    hyperbolic cosine integral

chi(x) = γ + ln(x) + ₀∫ˣ ((cosh(t) - 1) / t) δt
where γ == 0.5772_1566_4901_5328_6060_6512_0901-
           (the Euler-Mascheroni constant)
""" chi

"""
    li(x)

    logarithmic integral

li(x) =  ᵤ∫ˣ (δt / ln(t))

         where u = 1.4513_6923_4883_3810_5028+
                   (Soldner's constant)
""" li

"""
    si(x)

    sine integral

si(x) = ₀∫ˣ (sin(t) / t) δt
""" si


text:
#these are all the control data, they hold the flags that guide the calculation
mutable struct InputData
    name   ::Symbol
    option ::Symbol
    data   ::Any
end
name(data::InputData) = data.name

@with_kw_noshow mutable struct DFInput{P <: Package}
    name     ::String
    dir      ::String
    flags    ::SymAnyDict
    data     ::Vector{InputData} = InputData[]
    execs    ::Vector{Exec}
    run      ::Bool
    outdata  ::SymAnyDict=SymAnyDict()
end
DFInput{P}(name, dir, flags, data, execs, run) where P<:Package = DFInput{P}(name, abspath(dir), flags, data, execs, run, SymAnyDict())

"""
    DFInput(template::DFInput, name, newflags...; excs=deepcopy(execs(template)), run=true, data=nothing, dir=copy(template.dir))

Creates a new `DFInput` from a template `DFInput`, setting the newflags in the new one.
"""
function DFInput(template::DFInput, name, newflags...; excs=deepcopy(execs(template)), run=true, data=nothing, dir=deepcopy(template.dir))
    newflags = Dict(newflags...)

    input          = deepcopy(template)
    input.name     = name
    input.execs    = excs
    input.run      = run
    input.dir      = dir
    setflags!(input, newflags..., print=false)

    if data != nothing
        for (name, (option, data)) in data
            setdata!(input, name, data, option=option, print=false)
        end
    end
    return input
end

name(input::DFInput)  = input.name
dir(input::DFInput)   = input.dir
flags(input::DFInput) = input.flags
setdir!(input::DFInput, dir) = (input.dir = dir)
name_ext(input::DFInput, ext)          = name(input) * ext
infilename(input::DFInput{QE})         = name_ext(input, ".in")
infilename(input::DFInput{Wannier90})  = name_ext(input, ".win")
infilename(input::DFInput{Elk})        = "elk.in"
outfilename(input::DFInput{QE})        = name_ext(input, ".out")
outfilename(input::DFInput{Wannier90}) = name_ext(input, ".wout")
outfilename(input::DFInput{Elk})       = "elk.out"
inpath(input::DFInput)                 = joinpath(dir(input),  infilename(input))
outpath(input::DFInput)                = joinpath(dir(input),  outfilename(input))

hasflag(i::DFInput, s::Symbol) = haskey(flags(i), s)

function flag(input::DFInput, flag::Symbol)
    if hasflag(input, flag)
        return input.flags[flag]
    end
end

Base.eltype(::DFInput{P}) where P = P
package(::DFInput{P}) where P = P

data(input::DFInput)  = input.data

execs(input::DFInput) = input.execs
hasexec(input::DFInput, ex::AbstractString) = exec(input, ex) != nothing
setflow!(input::DFInput, run) = input.run = run

"Runs through all the set flags and checks if they are allowed and set to the correct value"
function cleanflags!(input::DFInput)
    for (flag, value) in flags(input)
        flagtype_ = flagtype(input, flag)
        if flagtype_ == Nothing
            @warn "Flag $flag was not found in allowed flags for exec $(execs(input)[2]). Removing flag."
            rmflags!(input, flag)
            continue
        end
        if !(isa(value, flagtype_) || eltype(value) <: flagtype_)
            try
                if isbitstype(eltype(value))
                    if length(value) > 1
                        flags(input)[flag] = convert(flagtype_, value)
                    else
                        flags(input)[flag] = convert(eltype(flagtype_), value)
                    end
                else
                    flags(input)[flag] = convert.(flagtype_, value)
                end
            catch
                error("Input $(name(input)): Could not convert :$flag of value $value to the correct type ($flagtype_), please set it to the correct type.")
            end
        end
    end
end

#TODO implement abinit and wannier90
"""
    sanitizeflags!(input::DFInput)

Cleans up flags, i.e. remove flags that are not allowed and convert all
flags to the correct types.
Tries to correct common errors for different input types.
"""
function sanitizeflags!(input::DFInput)
    cleanflags!(input)
end

function sanitizeflags!(input::DFInput{QE})
    setflags!(input, :outdir => "$(joinpath(dir(input), "outputs"))", print=false)
    cleanflags!(input)
    if isvcrelaxcalc(input)
	    #this is to make sure &ions and &cell are there in the input
	    !hasflag(input, :ion_dynamics)  && setflags!(input, :ion_dynamics  => "bfgs", print=false)
	    !hasflag(input, :cell_dynamics) && setflags!(input, :cell_dynamics => "bfgs", print=false)
    end
    #TODO add all the required flags
end

function setoradd!(datas::Vector{InputData}, data::InputData)
    found = false
    for (i, d) in enumerate(datas)
        if d.name == data.name
            datas[i] = data
            found = true
            break
        end
    end
    if !found
        push!(datas, data)
    end
end

"""
    setdata!(input::DFInput, data::InputData)

Adds the given data to the input. Should put it in the correct arrays.
"""
function setdata!(input::DFInput, data::InputData)
    setoradd!(input.data, data)
    return input
end

isbandscalc(input::DFInput{QE})    = flag(input, :calculation) == "bands"
isbandscalc(input::DFInput{Elk})   = input.name == "20"
isbandscalc(input::DFInput)        = false

isnscfcalc(input::DFInput{QE})     = flag(input, :calculation) == "nscf"
isnscfcalc(input::DFInput{Elk})    = input.name == "elk2wannier" #nscf == elk2wan??
isnscfcalc(input::DFInput)         = false

isscfcalc(input::DFInput{QE})      = flag(input, :calculation) == "scf"
isscfcalc(input::DFInput{Elk})     = input.name ∈ ["0", "1"]
isscfcalc(input::DFInput)          = false

iscolincalc(input::DFInput{QE})    = flag(input, :nspin) == 2
iscolincalc(input::DFInput)        = false

isnoncolincalc(input::DFInput{QE}) = flag(input, :noncolin) == true
isnoncolincalc(input::DFInput)     = false

isvcrelaxcalc(input::DFInput{QE})  = flag(input, :calculation) == "vc-relax"
isvcrelaxcalc(input::DFInput)      = false

isprojwfccalc(input::DFInput{QE})  = hasexec(input, "projwfc.x")
isprojwfccalc(input::DFInput)      = false

isdftucalc(input::DFInput{QE})     = flag(input, :lda_plus_u) == true
isdftucalc(input::DFInput)         = false

ismagneticcalc(input::DFInput{QE}) = flag(input, :nspin) ∈ [2, 4] || (flag(input, :lda_plus_u) == true && flag(input, :noncolin) == true)
ismagneticcalc(input::DFInput)     = false

#TODO review this!
outdata(input::DFInput) = input.outdata
hasoutput(input::DFInput) = !isempty(outdata(input))

hasoutfile(input::DFInput) = ispath(outpath(input))
hasnewout(input::DFInput, time) = mtime(outpath(input)) > time

readoutput(input::DFInput{QE}) = qe_read_output(input)
readoutput(input::DFInput{Wannier90}) = wan_read_output(outpath(input))

pseudodir(input::DFInput{QE}) = flag(input, :pseudo_dir)

setcutoffs!(input::DFInput, args...) = @warn "Setting cutoffs is not implemented for package $(package(input))"
setcutoffs!(input::DFInput{QE}, ecutwfc, ecutrho) = setflags!(input, :ecutwfc => ecutwfc, :ecutrho=>ecutrho)

function Emin_from_projwfc(structure::AbstractStructure, projwfc::DFInput{QE}, threshold::Number)
    hasoutput_assert(projwfc)
    hasexec_assert(projwfc, "projwfc.x")
    if !haskey(outdata(projwfc), :states)
        states, bands = qe_read_projwfc(outpath(projwfc))
        outdata(projwfc)[:states] = states
        outdata(projwfc)[:bands]  = bands
    else
        states, bands = outdata(projwfc)[:states], outdata(projwfc)[:bands]
    end

    mask = zeros(length(states))
    for (atid, at) in enumerate(atoms(structure))
        projs = projections(at)

        if isempty(projs)
            continue
        end

        stateids = Int[]
        for proj in projs
            orb = orbital(proj)
            push!.((stateids,), findall(x -> x.atom_id == atid && x.l == orb.l, states))
        end
        mask[stateids] .= 1.0

    end
    Emin = -10000.0
    for b in bands
        ψ = mean(b.extra[:ψ])
        tot_relevant_occupation = dot(mask, ψ)

        if tot_relevant_occupation > threshold
            Emin = minimum(b.eigvals)
            break
        end

    end
    if Emin == -10000.0
        error("Couldn't find any band with occupation of relevant projections above $threshold, were any set in the structure?")
    end
    return Emin
end

#asserts
function iscalc_assert(i::DFInput{QE}, calc)
    @assert flag(i, :calculation) == calc "Please provide a valid '$calc' calculation."
end

function hasoutput_assert(i::DFInput)
    @assert hasoutfile(i) "Please specify an input that has an outputfile."
end

function hasexec_assert(i::DFInput, exec::String)
    @assert hasexec(i, exec) "Please specify an input with $exec as it's executable."
end

#TODO Temporary handling of HubbardU situation
function set_hubbard_flags!(input::DFInput{QE}, str::AbstractStructure{T}) where {T}
	u_ats = unique(atoms(str))
    isdftucalc = any(x -> dftu(x).U != 0 || dftu(x).J0 != 0.0 || sum(dftu(x).J) != 0, u_ats)

	if isdftucalc
		Jmap = map(x -> dftu(x).J, u_ats)
		Jdim = maximum(length.(Jmap))
		Jarr = zeros(Jdim, length(u_ats))
		for (i, J) in enumerate(Jmap)
			diff = Jdim - length(J)
			if diff > 0
				for d in 1:diff
					push!(J, zero(eltype(J)))
				end
			end
			Jarr[:, i] .= J
		end
    	setflags!(input,
    	          :lda_plus_u    => true,
    	          :Hubbard_U     => map(x -> dftu(x).U, u_ats),
    	          :Hubbard_alpha => map(x -> dftu(x).α , u_ats),
    	          :Hubbard_beta  => map(x -> dftu(x).β , u_ats),
    	          :Hubbard_J     => Jarr,
    	          :Hubbard_J0    => map(x -> dftu(x).J0, u_ats);
	              print=false)
	end
end

function set_starting_magnetization_flags!(input::DFInput{QE}, str::AbstractStructure{T}) where {T}
	u_ats = unique(atoms(str))
	mags  = magnetization.(u_ats)
	starts= T[]
	θs    = T[]
	ϕs    = T[]
    ismagcalc      = any(x -> sum(magnetization(x)) != 0, u_ats)
    isnc =  isnoncolincalc(input) || any(x -> sum(magnetization(x)) !=0 && !isapprox(abs(normalize(magnetization(x)) ⋅ [0,0,1]), 1.0), u_ats)
	if ismagcalc && isnc
		for m in mags
			tm = normalize(m)
			if norm(m) == 0
				push!.((starts, θs, ϕs), 0.0)
			else
				θ = acos(tm[3]) * 180/π
				ϕ = atan(tm[2], tm[1])   * 180/π
				start = norm(m)
				push!(θs, θ)
				push!(ϕs, ϕ)
				push!(starts, start)
			end
		end
		setflags!(input, :noncolin=>true; print=false)
	elseif ismagcalc
		for m in mags
			push!.((θs, ϕs), 0.0)
			if norm(m) == 0
				push!(starts, 0)
			else
				push!(starts, sign(sum(m))*norm(m))
			end
		end
		setflags!(input, :nspin => 2; print=false)
	end
	setflags!(input, :starting_magnetization => starts, :angle1 => θs, :angle2 => ϕs; print=false)
end


text:
@testset "loss" begin

    arrtype = KnetLayers.arrtype
    x = arrtype(randn(10,10))
    indices = rand(1:10,10)
    l1 = CrossEntropyLoss(dims=2)
    l1(x,indices);
    @test true

    #https://github.com/denizyuret/Knet.jl/pull/374
    #x = arrtype(rand(10))
    #indices = rand([0,1],10)
    #l2 = BCELoss()
    #l2(x,indices)
    #@test true

    #x = arrtype(randn(10))
    #indices = rand([-1,1],10)
    #l3 = LogisticLoss()
    #l3(x,indices)
    #@test true

end


text:
"""
"""
struct DiscreteModelPortion{Dc,Dp} <: DiscreteModel{Dc,Dc}
  model::DiscreteModel{Dc,Dp}
  oldmodel::DiscreteModel{Dc,Dp}
  d_to_dface_to_old_dface::Vector{Vector{Int}}
end

get_grid(model::DiscreteModelPortion) = get_grid(model.model)

get_grid_topology(model::DiscreteModelPortion) = get_grid_topology(model.model)

get_face_labeling(model::DiscreteModelPortion) = get_face_labeling(model.model)

get_face_to_oldface(model::DiscreteModelPortion,d::Integer) = model.d_to_dface_to_old_dface[d+1]

get_cell_to_oldcell(model::DiscreteModelPortion) = get_face_to_oldface(model,num_cell_dims(model))

get_oldmodel(model::DiscreteModelPortion) = model.oldmodel

"""
"""
function DiscreteModelPortion(model::DiscreteModel, cell_to_oldcell::Vector{Int})
  grid_p =  GridPortion(get_grid(model),cell_to_oldcell)
  topo_p = GridTopology(grid_p)
  labels_p, d_to_dface_to_old_dface = _setup_labels_p(model,topo_p,cell_to_oldcell)
  model_p = DiscreteModel(grid_p,topo_p,labels_p)
  DiscreteModelPortion(model_p,model,d_to_dface_to_old_dface)
end

function DiscreteModelPortion(model::DiscreteModel, cell_to_is_in::Vector{Bool})
  cell_to_oldcell = findall(cell_to_is_in)
  DiscreteModelPortion(model,cell_to_oldcell)
end

function _setup_labels_p(model,topo_p,cell_to_oldcell)

  topo = get_grid_topology(model)
  labels = get_face_labeling(model)
  D = num_cell_dims(model)

  T = eltype(cell_to_oldcell)
  d_to_dface_to_old_dface = Vector{Vector{T}}(undef,D+1)
  for d in 0:(D-1)
    cell_lface_to_face = Table(get_faces(topo_p,D,d))
    ocell_lface_to_oface = Table(get_faces(topo,D,d))
    nfaces = num_faces(topo_p,d)
    d_to_dface_to_old_dface[d+1] = _find_face_to_oldface(
      cell_lface_to_face,
      ocell_lface_to_oface,
      cell_to_oldcell,
      nfaces)
  end
  d_to_dface_to_old_dface[D+1] = cell_to_oldcell

  d_to_dface_to_entity = [
    zeros(eltype(labels.d_to_dface_to_entity[d+1]),num_faces(topo_p,d)) for  d in 0:D]
  for d in 0:D
    face_to_label = d_to_dface_to_entity[d+1]
    oface_to_label = labels.d_to_dface_to_entity[d+1]
    face_to_oface = d_to_dface_to_old_dface[d+1]
    _update_labels_dim!(face_to_label,oface_to_label,face_to_oface)
  end

  labels_p = FaceLabeling(d_to_dface_to_entity,labels.tag_to_entities,labels.tag_to_name)

  labels_p, d_to_dface_to_old_dface
end

function _find_face_to_oldface(
  cell_lface_to_face::Table, ocell_lface_to_oface::Table, cell_to_ocell, nfaces)

  # Assumes that the local faces have the same order in new and old

  T = eltype(cell_to_ocell)
  face_to_oface = zeros(T,nfaces)

  for cell in 1:length(cell_lface_to_face)
    ocell = cell_to_ocell[cell]
    a = cell_lface_to_face.ptrs[cell]-1
    b = cell_lface_to_face.ptrs[cell+1]
    nlfaces = b - (a+1)
    c = ocell_lface_to_oface.ptrs[ocell]-1
    for lface in 1:nlfaces
      face = cell_lface_to_face.data[a+lface]
      oface = ocell_lface_to_oface.data[c+lface]
      face_to_oface[face] = oface
    end
  end

  face_to_oface
end

function _update_labels_dim!(face_to_label,oface_to_label,face_to_oface)
  for face in 1:length(face_to_label)
    oface = face_to_oface[face]
    label = oface_to_label[oface]
    face_to_label[face] = label
  end
end


text:
export auto_play

"""
    auto_play(answer::String; scorer::Function = min_max_score, verbose::Bool = true)
    auto_play(; scorer::Function = min_max_score, verbose::Bool = true)
    auto_play(a::Int; scorer::Function = min_max_score)

Computer solves a Wordle puzzle without user inputs. The number of steps to solve is returned.

* In the first case, `answer` is a five-letter word given by the user.
* In the second case, the five-letter word is randomly chosen by the computer.
* In the third case, an integer index into `ANS_LIST` is the input, and `verbose` is set to `false`.

The `scorcer` function is one of:
* `min_max_score`
* `entropy_score`
* `equi_score`

The first guess is chosen from the dictionary `FIRST_GUESS` depending on the `scorer` function.
"""
function auto_play(answer::String; scorer::Function = min_max_score, verbose::Bool = true)
    a = ANS_DICT[uppercase(answer)]  # this is the word we want to guess; throws an error if bad word

    stop = code_2_byte((2, 2, 2, 2, 2))
    history = Dict{Int,Byte}()

    g = FIRST_GUESS[scorer]
    if scorer == equi_score
        g = mod1(rand(Int), NG)
    end

    x = fast_wordle_score(g, a)
    history[g] = x
    if verbose
        wordle_print(g, x)
    end

    while x != stop
        possibles = filter_answers(history)
        g = best_guess(possibles, scorer)
        x = fast_wordle_score(g, a)
        if verbose
            wordle_print(g, x)
        end
        history[g] = x
    end
    return length(history)
end

function auto_play(; scorer::Function = min_max_score, verbose::Bool = true)
    a = mod1(rand(Int), NA)
    answer = ANS_LIST[a]
    auto_play(answer, scorer = scorer, verbose = verbose)
end

function auto_play(a::Int; scorer::Function = min_max_score)
    answer = ANS_LIST[a]
    auto_play(answer, scorer = scorer, verbose = false)

end

export score_all

"""
    score_all(scorer::Function = min_max_score)
Use the selected `scorer` function to solve all possible Wordle words
once each and return a vector containing the number of steps for each
solution.
"""
function score_all(scorer::Function = min_max_score)
    result = zeros(Int, NA)
    PM = Progress(NA)
    for a = 1:NA
        result[a] = auto_play(a, scorer = scorer)
        next!(PM)
    end
    return result
end

export all_play

scorer_list = [min_max_score, entropy_score, L2_score, equi_score]

export all_play

"""
    all_play(word::String)
Compare all solving methods to solve for the word.
"""
function all_play(answer::String)
    for s in scorer_list
        println("\nSolving with $s")
        steps = auto_play(answer, scorer=s)
        println("Solution in $steps steps")
    end
end


text:
function updateBuoyancy!(
    mb :: ModelBlock
)

    _b    = mb.fi._b
    _TEMP = mb.fi.sv[:_TEMP]
    _SALT = mb.fi.sv[:_SALT]

    @. _b = TS2b(_TEMP, _SALT)

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
struct SemCor{S}
    filepaths::Vector{S}
end

function init_datadeps(::Type{SemCor})
    for (ver, checksum) in [("1.6", "16814254fe194d55a2fcc24858aa76d71de3c49e495bd98478cc7345e766d8b7"),
                ("1.7", "0495577ac3a87c2a64fe6189798ea046de0f44943dfb7b60fe38cf648d34c421"),
                ("1.7.1", "70b9eb7ca0dc9d67655f9d671d40be10aeff490f0bea4f10cb1946127b74c102"),
                ("2.0", "93fbae725f0125dedb7369403fda1dace85b2dcd8a523ed80af23e863b18ef2c"),
                ("2.1", "0714f07dbcb84a215d668f3ee85892fa8fa4a8154439662eb7529413367b8f56"),
                ("3.0", "a8000014d6fc864f8bd9d83c62be601151cadd617c6554a39a1ad38b4b3f017b")]

        register(DataDep("SemCor $ver",
            """
            Website: http://web.eecs.umich.edu/%7Emihalcea/downloads.html#semcor
            Orignal Author: George A. Miller et al.
            Maintainer: Rada Mihalcea
            For WordNet version $ver

            This is SemCor orginally developed along side WordNet.
            It was automatically updated with the new numbering beyond WordNet 1.6,  by Rada Mihalcea

            Please cite the following publication if you use the corpora:
            George A. Miller, Claudia Leacock, Randee Tengi, and Ross T. Bunker. (1993). "A Semantic Concordance." In: Proceedings of the 3 DARPA Workshop on Human Language Technology.
            """,
            "http://web.eecs.umich.edu/~mihalcea/downloads/semcor/semcor$(ver).tar.gz",
            checksum;
            post_fetch_method = fn -> begin
                unpack(fn)
                innerdir = "semcor$(ver)"
                innerfiles = readdir(innerdir)
                # Move everything to current directory, under same name
                mv.(joinpath.(innerdir, innerfiles), innerfiles)
                rm(innerdir)
            end
        ))
    end
end

"""
    SemCor()

Creates a SemCor instance for lazy loading the corpus.

# Example Usage

```jldoctest
julia> corp = load(SemCor())
Channel{CorpusLoaders.Document{Array{Array{Array{CorpusLoaders.TaggedWord,1},1},1},String}}(sz_max:16,sz_curr:16)

julia> index = Dict{String, Vector{Vector{String}}}()
Dict{String,Array{Array{String,1},1}} with 0 entries

julia> for sentence in flatten_levels(corp, (!lvls)(SemCor, :sent, :word))
           multisense_words = filter(x->x isa SenseAnnotatedWord,  sentence)
           context_words = word.(sentence)
           for target in multisense_words
               word_sense = sensekey(target)
               uses = get!(Vector{Vector{String}}, index, word_sense)
               push!(uses, context_words)
           end
       end

julia> index["abandon%1:07:00::"]
1-element Array{Array{String,1},1}:
 String["in", "his", "motherland", ";", "in", "the", "spacious", "hunting_grounds", "of", "``"  …  ",", "color", ",", "dance", ",", "design", ",", "and", "thought", "."]

julia> index["abandon%2:38:00::"]
2-element Array{Array{String,1},1}:
 String["Jonathan", "wrote", "grimly", "of", "the", "destruction", "of", "Harpers_Ferry", "before", "they"  …  "First_Brigade", "had", "destroyed", "all", "the", "rolling_stock", "of", "the", "B_+_O_Railroad", "."]
 String["The", "expense", "of_this", "type", "of", "organization", "in", "religious", "life", ","  …  "in", "congregations", "illumines", "the", "nature", "of", "the", "Protestant", "development", "."]

```

SemCor is a classical Sense Annotated Corpus.

It has a structure of documents, paragraphs, sentences, words.

The words are either tagged with part of speech, or tagged with full lemma, part of speech and sensekey.

We expose the word function to get back the words, and the sensekey function to extract the sense key.

Texts semantically annotated with WordNet 1.6 senses (created at Princeton University), and automatically mapped to WordNet 1.7, WordNet 1.7.1, WordNet 2.0, WordNet 2.1, WordNet 3.0 to get:
SemCor 1.7
SemCor 1.7.1
SemCor 2.0
SemCor 2.1
SemCor 3.0

Please cite the following publication if you use the corpora:
        George A. Miller, Claudia Leacock, Randee Tengi, and Ross T. Bunker. (1993). "A Semantic Concordance." In: Proceedings of the 3 DARPA Workshop on Human Language Technology.

"""
function SemCor(dirpath)
    @assert(isdir(dirpath), dirpath)

    if "brown1" ∈ readdir(dirpath) && "brown2" ∈ readdir(dirpath)
        tagfile_paths = joinpath.(dirpath, ["brown1", "brown2"], "tagfiles")

        paths = joinpath.(tagfile_paths[1], readdir(tagfile_paths[1]))
        append!(paths, joinpath.(tagfile_paths[2], readdir(tagfile_paths[2])))

    elseif "tagfiles" ∈ readdir(dirpath)
        paths = joinpath.(dirpath,"tagfiles", readdir(joinpath(dirpath, "tagfiles")))
    else
        paths = joinpath.(dirpath, readdir(dirpath))
    end
    SemCor(paths)
end

SemCor() = SemCor(datadep"SemCor 3.0")

MultiResolutionIterators.levelname_map(::Type{SemCor}) = [
    :doc=>1, :contextfile=>1, :context=>1, :document=>1,
    :para=>2, :paragraph=>2,
    :sent=>3, :sentence=>3,
    :word=>4, :token=>4,
    :char=>5, :character=>5
	]

function parse_sense_annotated_word(line::AbstractString)
    captures = match(r"<wf cmd=done.* pos=(.*) lemma=(.*) wnsn=(.*) lexsn=(\d.*:\d*).*>(.*).*</wf>", line).captures
	pos, lemma, wnsn, lexsn, word = captures
	if ';' in wnsn
		# Discard Extra Senses
		wnsn = split(wnsn, ';') |> first
		lexsn = split(lexsn, ';') |> first
	end

	SenseAnnotatedWord(pos, lemma, parse(Int, wnsn), lexsn, word)
end

function parse_tagged_word(line::AbstractString)
    captures = match(r"<wf cmd=.* pos=(.*).*>(.*)</wf>", line).captures
    PosTaggedWord(captures...)
end

function parse_punc(line::AbstractString)
    captures = match(r"<punc>(.*)</punc>", line).captures
    PosTaggedWord("PUNC", first(captures))
end

function parse_semcorfile(filename)
    local sent
    local para
    paras = @NestedVector(TaggedWord,3)()
    context = Document(intern(basename(filename)), paras)

    # structure
    function new_paragraph(line)
        para = @NestedVector(TaggedWord,2)()
        push!(paras, para)
    end

    function new_sentence(line)
        sent = @NestedVector(TaggedWord,1)()
        push!(para, sent)
    end

    # words
    get_tagged(line) = push!(sent, parse_tagged_word(line))
    get_sense_annotated(line) = push!(sent, parse_sense_annotated_word(line))
	get_punc(line) = push!(sent, parse_punc(line))

    # Parse
    subparsers = [
        "<wf cmd=tag"=>          get_tagged,
        "<wf cmd=ignore"=>       get_tagged,
        "<wf cmd=done"=>         line -> occursin("lemma=", line) ? get_sense_annotated(line) : get_tagged(line),
        "<punc>"=>               get_punc,
        "<context"=>             ignore,
        "</context" =>           ignore,
        "<p" =>                  new_paragraph,
        "</p" =>                 ignore,
        "<s" =>                  new_sentence,
        "</s" =>                 ignore
    ]

    apply_subparsers(filename,subparsers)

    return context
end

function load(corpus::SemCor, doc_buffersize=16)
    Channel(;ctype=Document{@NestedVector(TaggedWord, 3), String}, csize=doc_buffersize) do ch
        for fn in corpus.filepaths
            doc = parse_semcorfile(fn)
            put!(ch, doc)
        end
    end
end


text:
using LinearAlgebra

function massmatrix(S::PBSpline{T}, q = S.p+1) where {T}
    M = zeros(T, S.nₕ, S.nₕ)
    for i in 1:S.nₕ
        for k in 0:S.p
            M[1,i] += integrate_gausslegendre(x -> (eval_PBSpline.(S, 1, x) .* eval_PBSpline.(S, i, x)),
                                                k*S.h, (k+1)*S.h, q)
        end
    end
    for j in 2:S.nₕ
        for i in 1:S.nₕ
            M[j,i] = M[j-1,(i-1+S.nₕ-1)%S.nₕ + 1]
        end
    end
    return M
end

function stiffnessmatrix(S::PBSpline{T}, q = S.p) where {T}
    K = zeros(T, S.nₕ, S.nₕ)
    for i in 1:S.nₕ
        for k in 0:S.p
            K[1,i] += integrate_gausslegendre(x -> (eval_deriv_PBSpline.(S, 1, x) .* eval_deriv_PBSpline.(S, i, x) ),
                                                k*S.h, (k+1)*S.h, q)
        end
    end
    for j in 2:S.nₕ
        for i in 1:S.nₕ
            K[j,i] = K[j-1,(i-1+S.nₕ-1)%S.nₕ + 1]
        end
    end
    return K
end

struct PoissonSolverPBSplines{DT <: Real} <: PoissonSolver{DT}
    p::Int
    nx::Int
    Δx::DT
    xgrid::Vector{DT}
    L::DT

    bspl::PBSpline{DT}

    M::Matrix{DT}
    S::Matrix{DT}
    Ŝ::Matrix{DT}

    P::Matrix{DT}
    R::Matrix{DT}

    ρ::Vector{DT}
    ϕ::Vector{DT}
    rhs::Vector{DT}

    Mfac::LU{DT, Matrix{DT}}
    Sfac::LU{DT, Matrix{DT}}

    function PoissonSolverPBSplines{DT}(p::Int, nx::Int, L::DT) where {DT}
        Δx = L/nx
        xgrid = collect(0:Δx:1)

        bspl = PBSpline(p, nx, L)
        M = massmatrix(bspl)
        S = stiffnessmatrix(bspl)

        A = ones(nx)
        R = A * A' / (A' * A)
        P = Matrix(I, nx, nx) .- R
        Ŝ = S .+ R

        new(p, nx, Δx, xgrid, L, bspl, M, S, Ŝ, P, R,
            zeros(DT,nx), zeros(DT,nx), zeros(DT,nx),
            lu(M), lu(Ŝ))
    end
end

PoissonSolverPBSplines(p::Int, nx::Int, L::DT) where {DT} = PoissonSolverPBSplines{DT}(p, nx, L)

Base.length(p::PoissonSolverPBSplines) = p.nx

function solve!(p::PoissonSolverPBSplines{DT}, x::AbstractVector{DT}, w::AbstractVector{DT} = one.(x) ./ length(x)) where {DT}
    rhs_particles_PBSBasis(x, w, p.bspl, p.rhs)
    ldiv!(p.ϕ, p.Sfac, - p.P * p.rhs)
    ldiv!(p.ρ, p.Mfac, p.rhs)
    return p
end

function eval_density(p::PoissonSolverPBSplines{DT}, x::DT) where {DT}
    eval_PBSBasis(p.ρ, p.bspl, x)
end

function eval_potential(p::PoissonSolverPBSplines{DT}, x::DT) where {DT}
    eval_PBSBasis(p.ϕ, p.bspl, x)
end

function eval_field(p::PoissonSolverPBSplines{DT}, x::DT) where {DT}
    eval_deriv_PBSBasis(p.ϕ, p.bspl, x)
end


text:
module HyperECNNs

using LinearAlgebra
using Quaternions
using Random
rng = MersenneTwister(1234);

### Bilinear form used in the paper!

function LambdaInner(U,x,Params)
    w = zeros(size(U,3),1)
    for i = 1:size(U,2)
        w = w + Params[i]*(U[:,i,:]'*x[:,i])
    end
    return w
end

### Some activation functions

function csign(x,K)
    z = x[:,1]+x[:,2]*im
    phase_quanta = (round.(K*(2*pi.+angle.(z))./(2*pi))).%K
    z = exp.(2.0*pi*phase_quanta*im/K)
    return hcat(real.(z),imag.(z))
end

function twincsign(x, Params = 16)
    c1 = csign(x[:,1:2],Params)
    c2 = csign(x[:,3:4],Params)
    return hcat(c1,c2)
end

function SplitSign(a, Params = nothing)
    x = zeros(size(a))
    for i = 1:size(a,3)
        x[:,:,i] = sign.(a[:,:,i])
    end
    return x
end

##########################################################################
### Generic Hypercomplex-Valued Networks
##########################################################################

function Sync(BilinearForm, BilinearFormParams, ActFunction, ActFunctionParams, U, xinput, alpha = 1, beta = 0, it_max = 1.e3)
    ### Exponential Hypercomplex RCNN
    Name = "Hypercomplex ECNN (Synchronous)"

    N = size(U,1)
    hyperN = size(U,2)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau
    # Compute the weights;
    w = exp.(alpha*BilinearForm(U,x,BilinearFormParams).+beta);

    # Compute the energy
    Energy = zeros(1)
    Energy[1]=-sum(w)/alpha
    a = zeros(N,hyperN)
    while (Error>tau)&&(it<it_max)
        it = it+1

        # Compute the quaternion-valued activation potentials;
        for j = 1:hyperN
            a[:,j] = U[:,j,:]*w
        end

        # Compute the next state;
        x = ActFunction(a, ActFunctionParams)

        # Compute the weights;
        w = exp.(alpha*BilinearForm(U,x,BilinearFormParams).+beta);

        append!(Energy,-sum(w)/alpha)
        Error = norm(x-xold)
        xold = copy(x)
    end
    if it_max<=it
        println(Name," failed to converge in ",it_max," iterations.")
    end
    return x, Energy
end

function Seq(BilinearForm, BilinearFormParams, ActFunction, ActFunctionParams, U, xinput, alpha = 1, beta = 0, it_max = 1.e3)
    ### Exponential Hypercomplex RCNN
    Name = "Hypercomplex ECNN (Asynchronous)"

    N = size(U,1)
    hyperN = size(U,2)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau
    Energy = zeros(1)
    EnergyAux = zeros(N)
    # Compute the weights;
    w = exp.(alpha*BilinearForm(U,x,BilinearFormParams).+beta);
    Energy[1]=-sum(w)/alpha
    a = zeros(1,hyperN)
    while (Error>tau)&&(it<it_max)
        it = it+1
        ind = randperm(rng, N)
        for i = 1:N
            # Compute the quaternion-valued activation potentials;
            for j = 1:hyperN
                a[1,j] = dot(U[ind[i],j,:],w)
            end

            # Compute the next state;
            x[ind[i],:] = ActFunction(a, ActFunctionParams)

            # Compute the weights;
            w = w.*exp.(alpha*BilinearForm(U[[ind[i],ind[i]],:,:],vcat(-xold[ind[i],:]',x[ind[i],:]'),BilinearFormParams));
            EnergyAux[i]=-sum(w)/alpha
        end
        append!(Energy,EnergyAux)
        Error = norm(x-xold)
        xold = copy(x)
    end
    if it_max<=it
        println(Name," failed to converge in ",it_max," iterations.")
    end
    return x, Energy
end

function Seq_slow(BilinearForm, BilinearFormParams, ActFunction, ActFunctionParams, U, xinput, alpha = 1, beta = 0, it_max = 1.e3)
    ### Exponential Hypercomplex RCNN
    Name = "Hypercomplex ECNN (Asynchronous)"

    N = size(U,1)
    hyperN = size(U,2)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau
    Energy = zeros(1)
    EnergyAux = zeros(N)
    # Compute the weights;
    w = exp.(alpha*BilinearForm(U,x,BilinearFormParams).+beta);
    Energy[1]=-sum(w)/alpha
    a = zeros(1,hyperN)
    while (Error>tau)&&(it<it_max)
        it = it+1
        ind = randperm(rng, N)
        for i = 1:N
            # Compute the quaternion-valued activation potentials;
            for j = 1:hyperN
                a[1,j] = dot(U[ind[i],j,:],w)
            end

            # Compute the next state;
            x[ind[i],:] = ActFunction(a, ActFunctionParams)

            # Compute the weights;
            w = exp.(alpha*BilinearForm(U,x,BilinearFormParams).+beta);
            EnergyAux[i]=-sum(w)/alpha
        end
        append!(Energy,EnergyAux)
        Error = norm(x-xold)
        xold = copy(x)
    end
    if it_max<=it
        println(Name," failed to converge in ",it_max," iterations.")
    end
    return x, Energy
end

##########################################################################
### Unit Quaternion Networks
##########################################################################
function UnitQ_Sync(U, xinput, alpha = 1, beta=0, it_max = 1.e3)
    ### Exponential Unit Quaternion RCNN
    Name = "Unit Quaternion ECNN (Synchronous)"

    N = size(U,1)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau

    # Compute the weights
    w = exp.(alpha*Array{Float64}(real((U'*x))).+beta);

    # Compute the energy
    Energy = zeros(Float64,(1,))
    Energy[1] = -sum(w)/alpha
    while (Error>tau)&&(it<it_max)
        it = it+1

        # Compute the next state
        a = U*w
        x = a ./ abs.(a)

        # Compute the weights and the energy
        w = exp.(alpha*Array{Float64}(real((U'*x))).+beta);
        append!(Energy,-sum(w)/alpha)

        Error = norm(x-xold)
        xold = copy(x)
    end
    if it_max<=it
        println(Name," failed to converge in ",it_max," iterations.")
    #else
    #    println(Name," converged in ",it," iterations using synchronous update.")
    end
    return x, Energy
end

function UnitQ_Seq(U, xinput, alpha = 1, beta=0, it_max = 1.e3)
    ### Exponential Unit Quaternion RCNN
    Name = "Unit Quaternion ECNN (Asynchronous)"

    N = size(U,1)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau

    # Compute the weights;
    w = exp.(alpha*Array{Float64}(real((U'*x))).+beta);

    # Compute the energy
    Energy = zeros(Float64,(1,))
    EnergyAux = zeros(N)
    Energy[1]=-sum(w)/alpha
    while (Error>tau)&&(it<it_max)
        it = it+1
        ind = randperm(rng, N)
        for i = 1:N
            # Compute the quaternion-valued activation potentials;
            a = dot(conj(U[ind[i],:]),w)

            # Compute the next state;
            x[ind[i]] = a/abs(a)

            # Compute the weights;
            w = exp.(alpha*Array{Float64}(real((U'*x))).+beta);
            EnergyAux[i]=-sum(w)/alpha
        end
        append!(Energy,EnergyAux)
        Error = norm(x-xold)
        xold = copy(x)
    end
    if it_max<=it
        println(Name," failed to converge in ",it_max," iterations.")
    end
    return x, Energy
end

end


text:
abstract type IntegratorCache{DT,D} end

abstract type ODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type DAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type IDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PODEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end
abstract type PDAEIntegratorCache{DT,D} <: IntegratorCache{DT,D} end

IntegratorCache(params::Parameters{DT}) where {DT} = IntegratorCache{DT}(params)
IntegratorCache{ST}(params::Parameters) where {ST} = error("IntegratorCache{ST}(params) not implemented for ", typeof(params))

IntegratorCache(integrator::Integrator{DT}) where {DT} = IntegratorCache{DT}(integrator)
IntegratorCache{ST}(integrator::Integrator) where {ST} = error("IntegratorCache{ST}(int)! not implemented for ", typeof(integrator))

CacheType(T, params::Parameters) = error("CacheType(T, params) not implemented for ", typeof(params))

struct CacheDict{PT}
    params::PT
    caches::Dict{UInt64, IntegratorCache}

    function CacheDict(params::Parameters)
        new{typeof(params)}(params, Dict{UInt64, IntegratorCache}())
    end
end

@inline function Base.getindex(c::CacheDict, ST::DataType)
    key = hash(Threads.threadid(), hash(ST))
    if haskey(c.caches, key)
        c.caches[key]
    else
        c.caches[key] = IntegratorCache{ST}(c.params)
    end::CacheType(ST, c.params)
end


text:
using Artifacts
using BinaryProvider
using CxxWrap

vizdoom_artifact_dir = artifact"ViZDoom"
vizdoom_dir = joinpath(vizdoom_artifact_dir, "ViZDoom-1.1.8")

# Copy our lib_julia to the src dir and remove any old ones
const deps_dir = @__DIR__
run(`rm -rf $(joinpath(vizdoom_dir, "src", "lib_julia"))`)
run(`cp -r $(joinpath(deps_dir, "lib_julia")) $(joinpath(vizdoom_dir, "src"))`)

const JLCxx_DIR = get(ENV, "JLCXX_DIR", joinpath(CxxWrap.prefix_path(), "lib", "cmake", "JlCxx"))
const cmake_cmd = `cmake -DJlCxx_DIR=$(JLCxx_DIR) -DBUILD_JULIA=ON -DJulia_EXECUTABLE=$(joinpath(Sys.BINDIR, "julia")) .`
@info cmake_cmd
cd(() -> run(cmake_cmd), vizdoom_dir)
# build with num_cores - 1. Saving 1 for the user
num_cores = parse(Int, read(`nproc`, String))
cd(() -> run(`make -j $(num_cores - 1)`), vizdoom_dir)

# Copy freedoom2.wad into bin where everything else is.
# The default config looks for ./freedoom2.wad so this allows for simpler configs
run(`cp $(joinpath(vizdoom_dir, "src", "freedoom2.wad")) $(joinpath(vizdoom_dir, "bin"))`)


text:
function linmin2( pw, Vpot, α_t, psi, d, g )

    psic = ortho_gram_schmidt(psi + α_t*d)
    gt = calc_grad( pw, Vpot, psic )
    denum = real(sum(conj(g-gt).*d))
    if denum != 0.0
      α = abs(α_t*real(sum(conj(g).*d))/denum)
    else
      α = 0.0
    end

    return α

end


text:
@cell "thick arc" arc(Point2f(0), 10f0, 0f0, pi, linewidth=20)

@cell "stroked rect poly" poly(Recti(0, 0, 200, 200), strokewidth=20, strokecolor=:red, color=(:black, 0.4))

@cell "array of rects poly" begin
    f, ax, pl = poly([Rect(0, 0, 20, 20)])
    scatter!(Rect(0, 0, 20, 20), color=:red, markersize=20)
    f
end

@cell "char marker scenespace" begin
    f, ax, pl = lines(Rect(0, 0, 1, 1), linewidth=4)
    scatter!([Point2f(0.5, 0.5)], markersize=1, markerspace=:data, marker='I')
    f
end

@cell "lines number color" lines(RNG.rand(10), RNG.rand(10), color=RNG.rand(10), linewidth=10)
@cell "lines array of colors" lines(RNG.rand(10), RNG.rand(10), color=RNG.rand(RGBAf, 10), linewidth=10)
@cell "scatter interval" scatter(0..1, RNG.rand(10), markersize=RNG.rand(10) .* 20)
@cell "scatter linrange" scatter(LinRange(0, 1, 10), RNG.rand(10))

@cell "scatter rotation" begin
    angles = range(0, stop=2pi, length=20)
    pos = Point2f.(sin.(angles), cos.(angles))
    f, ax, pl = scatter(pos, markersize=0.2, markerspace=:data, rotations=-angles, marker='▲', axis=(;aspect = DataAspect()))
    scatter!(pos, markersize=10, color=:red)
    f
end

@cell "heatmap transparent colormap" heatmap(RNG.rand(50, 50), colormap=(:RdBu, 0.2))

@cell "contour small x" contour(RNG.rand(10, 100))
@cell "contour small y" contour(RNG.rand(100, 10))
@cell "contour with levels" contour(RNG.randn(100, 90), levels=3)

@cell "contour with levels array" contour(RNG.randn(100, 90), levels=[0.1, 0.5, 0.8])
@cell "contour with color per level" contour(RNG.randn(33, 30), levels=[0.1, 0.5, 0.9], color=[:black, :green, (:blue, 0.4)], linewidth=2)

@cell "contour with colorrange" contour(
    RNG.rand(33, 30) .* 6 .- 3, levels=[-2.5, 0.4, 0.5, 0.6, 2.5],
    colormap=[(:black, 0.2), :red, :blue, :green, (:black, 0.2)],
    colorrange=(0.2, 0.8)
)

@cell "circle line" lines(Circle(Point2f(0), Float32(1)))

@cell "streamplot with func" begin
    v(x::Point2{T}) where T = Point2{T}(x[2], 4 * x[1])
    streamplot(v, -2..2, -2..2, arrow_size=10)
end

@cell "lines with func" lines(-1..1, x -> x^2)
@cell "scatter with func" scatter(-1..1, x -> x^2)

@cell "volume translated" begin
    r = range(-3pi, stop=3pi, length=100)
    fig, ax, vplot = Makie.volume(r, r, r, (x, y, z) -> cos(x) + sin(y) + cos(z), algorithm=:iso, isorange=0.1f0, show_axis=false)
    v2 = volume!(ax, r, r, r, (x, y, z) -> cos(x) + sin(y) + cos(z), algorithm=:mip,
                 transformation=(translation=Vec3f(6pi, 0, 0),))
    fig
end

@cell "meshscatter color numbers" meshscatter(RNG.rand(10), RNG.rand(10), RNG.rand(10), color=RNG.rand(10))
@cell "meshscatter color array" meshscatter(RNG.rand(10), RNG.rand(10), RNG.rand(10), color=RNG.rand(RGBAf, 10), transparency=true)

@cell "transparent mesh texture" begin
    s1 = uv_mesh(Sphere(Point3f(0), 1f0))
    f, ax, pl = mesh(uv_mesh(Sphere(Point3f(0), 1f0)), color=RNG.rand(50, 50))
    # ugh, bug In GeometryTypes for UVs of non unit spheres.
    s2 = uv_mesh(Sphere(Point3f(0), 1f0))
    s2.position .= s2.position .+ (Point3f(0, 2, 0),)
    mesh!(s2, color=RNG.rand(RGBAf, 50, 50))
    f
end

@cell "Unequal x and y sizes in surface" begin
    NL = 15
    NR = 31
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    lspace = range(-10, stop=10, length=NL)
    rspace = range(-10, stop=10, length=NR)

    z = Float32[xy_data(x, y) for x in lspace, y in rspace]
    l = range(0, stop=3, length=NL)
    r = range(0, stop=3, length=NR)
    surface(
        l, r, z,
        colormap=:Spectral
    )
end

@cell "Matrices of data in surfaces" begin
    NL = 30
    NR = 31
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    lspace = range(-10, stop=10, length=NL)
    rspace = range(-10, stop=10, length=NR)

    z = Float32[xy_data(x, y) for x in lspace, y in rspace]
    l = range(0, stop=3, length=NL)
    r = range(0, stop=3, length=NR)
    surface(
        [l for l in l, r in r], [r for l in l, r in r], z,
        colormap=:Spectral
    )
end

@cell "heatmaps & surface" begin
    data =
        hcat(LinRange(2, 3, 4), LinRange(2, 2.5, 4), LinRange(2.5, 3, 4), [1, NaN, NaN, 5])

    fig = Figure()
    heatmap(
        fig[1, 1],
        data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
    )
    surface(
        fig[1, 2],
        zeros(size(data)),
        color = data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
        shading = false,
    )
    surface!(
        Axis(fig[2, 2]),
        data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
        shading = false,
    )
    fig
end

@cell "log10 heatmap" begin
    heatmap(RNG.rand(10, 5), axis = (yscale = log10, xscale=log10))
end

@cell "reverse range heatmap" begin
    x = [1 0
         2 3]
    heatmap(1:2, 1:-1:0, x)
end

@cell "lines linesegments width test" begin
    res = 200
    s = Scene(camera=campixel!, resolution=(res, res))
    half = res / 2
    linewidth = 10
    xstart = half - (half/2)
    xend = xstart + 100
    half_w = linewidth/2

    lines!(s, Point2f[(xstart, half), (xend, half)], linewidth=linewidth)
    scatter!(s, Point2f[(xstart, half + half_w), (xstart, half - half_w), (xend, half + half_w), (xend, half - half_w)], color=:red, markersize=2)

    l2 = linesegments!(s, Point2f[(xstart, half), (xend, half)], linewidth=linewidth, color=:gray)
    s2 = scatter!(s, Point2f[(xstart, half + half_w), (xstart, half - half_w), (xend, half + half_w), (xend, half - half_w)], color=:red, markersize=2)

    for p in (l2, s2)
        translate!(p, 0, 20, 0)
    end

    s
end

@cell "multipoly with multi strokes" begin
    P = Polygon.([Point2f[[0.45, 0.05], [0.64, 0.15], [0.37, 0.62]],
         Point2f[[0.32, 0.66], [0.46, 0.59], [0.09, 0.08]]])
    poly(P, color = [:red, :green], strokecolor = [:blue, :red], strokewidth = 2)
end

@cell "fast pixel marker" begin
    scatter(RNG.rand(Point2f, 10000), marker=Makie.FastPixel())
end

@cell "axsi3" begin
    meshscatter(RNG.rand(Point3f, 10), axis=(type=Axis3,))
end


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_input(id = "input-3", value = "initial value", type = "text"),
    html_div(id="output-1")
end

callback!(app, Output("output-1", "children"), Input("input-3", "value")) do input_value
    "You've entered $(input_value)"
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using Pkg
pkg"activate ."
using Revise
using Layered
using Colors

function example()
    c, l = canvas(4.5, 3)
    n = 5

    defaultfont("Helvetica Neue Light")

    sls = layer!.((r, i) -> begin
        margin = 20
        avail_w = (r.width - 2margin) / n
        avail_h = (r.height - 2margin) / n
        Transform(translation = topleft(r) + P(margin + (i-0.5) * avail_w, margin + (i-0.5) * avail_h))
    end, l, c.rect, 1:n)

    rs = rect!.((r, i) -> begin
        Rect(P(0, 0), 80, 50, deg(0))
    end, sls, c.rect, 0:n-1) .+ Fill(Gray(0.5)) .+ Visible(true)

    crosses = polygon!.(r -> begin
        ncross(r.center, 4, 3, 0.3)
    end, sls, rs) .+ Fill("black") .+ Stroke("transparent")

    circlesleft = circle!.(r -> begin
        Circle(r.center - X(r.width/4), 10)
    end, sls, rs) .+ Linestyle(:dashed) .+
        Stroke.([i == 2 ? "lime" : "black" for i in 1:n])

    circlesright = circle!.(r -> begin
        Circle(r.center + X(r.width/4), 10)
    end, sls, rs) .+ Linestyle(:dashed) .+
        Stroke.([i == 4 ? "tomato" : "black" for i in 1:n])

    focuscircle = circle!(r -> begin
        Circle(topright(r) + P(-35, 35), 30)
    end, l, c.rect) + Linestyle(:dashed)

    textl = layer!(l, Transform())

    phases = txt!.((r, phase) -> begin
        Txt(bottomleft(r) + Y(0), "phase $phase", 10, :l, :t)
    end, textl, rs, 1:5)

    path!((c1, c2) -> begin
        Path(false, outertangents(c1, c2)...)
    end, l, focuscircle, circlesright[4]) + Opacity(0.3)#+ Linestyle(:dashed)

    arr = path!(l, rs[1], rs[end]) do r1, r2
        a = arrow(topright(r1), topright(r2), 9, 9, 1, 1, 0)
        a + normal(a.segments[1], -5)
    end + Fill("tomato") + Stroke("transparent")

    te = txt!(textl, arr) do arr
        pos = normfrom(arr.segments[1], 0.5, -4)
        Txt(pos, "time", 10, :c, :b, angle(arr.segments[1]))
    end

    r = rect_pre!(textl, te, te) do t
        Rect(t, 4)
    end + Fill("tomato") + Stroke("transparent")

    txt!(textl, focuscircle) do fc
        Txt(fc.center, "focus circle", 10, :c, :c)
    end

    txt!(textl, c.rect) do r
        Txt(bottomleft(r) + P(5, -5), "experimental paradigm", 12, :l, :b)
    end

    c

    # draw_svg(c, "example.svg")

end; example()


text:
using DashBootstrapComponents, DashHtmlComponents

card_content = [
    dbc_cardheader("Card header"),
    dbc_cardbody([
        html_h5("Card title", className = "card-title"),
        html_p("This is some card content that we'll reuse", className = "card-text"),
    ]),
];

cards = html_div([
    dbc_row(
        [
            dbc_col(dbc_card(card_content, color = "primary", inverse = true)),
            dbc_col(dbc_card(card_content, color = "secondary", inverse = true)),
            dbc_col(dbc_card(card_content, color = "info", inverse = true)),
        ],
        className = "mb-4",
    ),
    dbc_row(
        [
            dbc_col(dbc_card(card_content, color = "success", inverse = true)),
            dbc_col(dbc_card(card_content, color = "warning", inverse = true)),
            dbc_col(dbc_card(card_content, color = "danger", inverse = true)),
        ],
        className = "mb-4",
    ),
    dbc_row([
        dbc_col(dbc_card(card_content, color = "light")),
        dbc_col(dbc_card(card_content, color = "dark", inverse = true)),
    ]),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module EmittersConsts
    ARRROW_LENGTH = 0.3
    ARRROW_SIZE = 0.01
    MARKER_SIZE = 1
end

#-------------------------------------
# draw debug information - local axes and positions
#-------------------------------------
function maybe_draw_debug_info(scene::Scene, so::AbstractSceneObject, o::Origins.AbstractOriginDistribution; transform::Geometry.Transform = Transform(), debug::Bool=false, kwargs...) where {T<:Real}

    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    if (debug)

        debug_so = EmptySceneObject(name="Debug Info")
        OpticSimVis.parent!(debug_so , so)

        # # this is a stupid hack to force makie to render in 3d - for some scenes, makie decide with no apperent reason to show in 2d instead of 3d
        # Makie.scatter!(scene, [pos[1], pos[1]+0.1], [pos[2], pos[2]+0.1], [pos[3], pos[3]+0.1], color=:red, markersize=0)

        # draw the origin and normal of the surface
        # Makie.scatter!(scene, pos, color=:blue, markersize = MARKER_SIZE * visual_size(o))
        # draw!(scene, [pos], size=0.02, color=:blue)

        # draw a point in the origin as a point cloud
        pc_mat = OpticSimVis.Material(color=RGBA(0.1, 0.1, 9.1, 1.0), size=0.05)
        pc = OpticSimVis.PointCloud(points=[pos], material=pc_mat, name="Origin Point")
        OpticSimVis.parent!(pc , debug_so)

        x_axis = uv * 0.5 * EmittersConsts.ARRROW_LENGTH * visual_size(o)
        y_axis = vv * 0.5 * EmittersConsts.ARRROW_LENGTH * visual_size(o)
        z_axis = dir * EmittersConsts.ARRROW_LENGTH * visual_size(o)
        axes = OpticSimVis.Axes(
            tr=Glimmer.transform(Transform(pos)),
            x_axis=x_axis,
            y_axis=y_axis,
            z_axis=z_axis,
            # material=axes_mat,
            name="Local Axes")
        OpticSimVis.parent!(axes , debug_so)

    end

end

#-------------------------------------
# draw point origin
#-------------------------------------
function draw!(scene::Scene, o::Origins.Point; parent_so::AbstractSceneObject = root(scene), kwargs...) where {T<:Real}
    transform = local_tr(parent_so)

    pos = origin(transform)
    obj = draw!(scene, [pos]; size=0.01, kwargs...)

    maybe_draw_debug_info(scene, obj[:scene_object], o; transform=transform, kwargs...)
    return obj
end

#-------------------------------------
# draw RectGrid and RectUniform origins
#-------------------------------------
function draw!(scene::Scene, o::Union{Origins.RectGrid, Origins.RectUniform, Origins.RectJitterGrid}; parent_so::AbstractSceneObject = root(scene), kwargs...) where {T<:Real}
    transform = local_tr(parent_so)
    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    plane = OpticSim.Plane(dir, pos)
    rect = OpticSim.Rectangle(plane, o.width / 2, o.height / 2, uv, vv)

    obj = draw!(scene, rect;  parent_so=parent_so, kwargs...)

    maybe_draw_debug_info(scene, obj[:scene_object], o; transform=transform, kwargs...)
    return obj
end

#-------------------------------------
# draw hexapolar origin
#-------------------------------------
function draw!(scene::Scene, o::Origins.Hexapolar{T}; parent_so::AbstractSceneObject = root(scene), kwargs...) where {T<:Real}
    transform = local_tr(parent_so)
    dir = forward(transform)
    uv = SVector{3}(right(transform))
    vv = SVector{3}(up(transform))
    pos = origin(transform)

    plane = OpticSim.Plane(dir, pos)
    ellipse = OpticSim.Ellipse(plane, o.halfsizeu, o.halfsizev, uv, vv)
    obj = draw!(scene, ellipse;  kwargs...)

    maybe_draw_debug_info(scene, obj[:scene_object], o; transform=transform, kwargs...)
    return obj
end

#-------------------------------------
# draw source
#-------------------------------------
function draw!(scene::Scene, s::Sources.Source{T}; parent_so::AbstractSceneObject = root(scene), debug::Bool=false, kwargs...) where {T<:Real}

    name = "Source-$(UUIDs.uuid1())"
    t = s.transform;
    root_so = EmptySceneObject(tr=t, name=name)
    OpticSimVis.parent!(root_so , parent_so)

    obj = draw!(scene, s.origins;  color=RGBA(1.0, 1.0, 0.0, 0.5), parent_so=root_so, debug=debug, kwargs...)

    if (debug)
        parent_transform = toTransform(local_tr(parent_so))

        m = zeros(T, length(s), 7)
        for (index, optical_ray) in enumerate(s)
            ray = OpticSim.ray(optical_ray)
            ray = parent_transform * ray
            m[index, 1:7] = [ray.origin... ray.direction... OpticSim.power(optical_ray)]
        end

        m[:, 4:6] .*= m[:, 7] * EmittersConsts.ARRROW_LENGTH * visual_size(s.origins)

        base_so = obj[:scene_object]

        points = Vector{SVector{3, Float64}}(undef, size(m)[1] * 2)
        for i in 1:size(m)[1]
            point = SVector(m[i,1], m[i, 2], m[i, 3])
            dir = SVector(m[i,4], m[i, 5], m[i, 6])
            point2 = point + dir
            index = (i-1) * 2 + 1
            points[index] = point
            points[index+1] = point2
        end

        # @info typeof(points)
        segments_mat = OpticSimVis.Material(color=RGBA(0.9, 0.9, 0.1, 0.5))
        segments = OpticSimVis.LineSegments(
            points=points,
            tr=transform(Transform(Vec3(0.0, 0.0, 0.0))),
            material=segments_mat,
            name="Debug Rays")
        OpticSimVis.parent!(segments, base_so)

        # debug_so = EmptySceneObject(name="Debug")
        # OpticSimVis.set_parent!(debug_so , base_so)

        # arrow_mat = OpticSimVis.Material(color=RGBA(0.7, 0.7, 0.1, 0.5))
        # for i in 1:size(m)[1]
        #     point = SVector(m[i,1], m[i, 2], m[i, 3])
        #     dir = SVector(m[i,4], m[i, 5], m[i, 6])
        #     point2 = point + dir
        #     arrow = OpticSimVis.Arrow(point, point2; material=arrow_mat, name="Arrow_$i")
        #     OpticSimVis.set_parent!(arrow , debug_so)
        # end

        # Makie.arrows!(scene, [Makie.Point3f0(origin(ray))], [Makie.Point3f0(rayscale * direction(ray))]; kwargs..., arrowsize = min(0.05, rayscale * 0.05), arrowcolor = color, linecolor = color, linewidth = 2)
        # color = :yellow
        # arrow_size = ARRROW_SIZE * visual_size(s.origins)
        # Makie.arrows!(scene, m[:,1], m[:,2], m[:,3], m[:,4], m[:,5], m[:,6]; kwargs...,  arrowcolor=color, linecolor=color, arrowsize=arrow_size, linewidth=arrow_size*0.5)
    end

    return obj
end

#-------------------------------------
# draw composite source
#-------------------------------------
function draw!(scene::Scene, s::Sources.CompositeSource{T}; parent_so::AbstractSceneObject = root(scene), kwargs...) where {T<:Real}

    name = "Composite Source" # -$(UUIDs.uuid1())"
    t = s.transform;
    root_so = EmptySceneObject(name=name)
    OpticSimVis.parent!(root_so , parent_so)

    # axes = OpticSimVis.Axes(
    #     # tr=Transform(pos),
    #     axes_scale=20,
    #     shaft_scale=5,
    #     name="Debug Axes")
    # OpticSimVis.parent!(axes , root_so)

    for source in s.sources
        draw!(scene, source; parent_so=root_so, kwargs...)
    end
end


text:
# ======================================================= #
#                                                         #
#                  TRIQS Covnersion                       #
#                                                         #
# ======================================================= #
function triqs_read_gf(file, name)
    freqSize = h5read(file, name)["mesh"]["size"]
    pos_only = h5read(file, name)["mesh"]["positive_freq_only"]
    arr_tmp = h5read(file, name)["data"]
    n0 = (1-pos_only)*(trunc(Int, freqSize/2))+1
    res = arr_tmp[1,n0:end] + 1im .* arr_tmp[2,n0:end]
end

function triqs_mesh_to_range(mesh)
	s = mesh["size"]
	p = mesh["positive_freq_only"]
	n0 = (1-p)*(trunc(Int, s/2))
	return (0:s-1) .- n0
end

function triqs_build_freqGrid(mesh)
	nu_range = triqs_mesh_to_range(mesh["MeshComponent2"])
	nup_range = triqs_mesh_to_range(mesh["MeshComponent1"])
	w_range = triqs_mesh_to_range(mesh["MeshComponent0"])
	return collect(Base.product(nu_range, nup_range, w_range))
end

function triqs_linearize_mesh(mesh)
	new_arr = Array{eltype(mesh),1}(undef, length(mesh))
	ind = 1
	for wi in 1:size(mesh,3)
		for nui in 1:size(mesh,1)
			for nupi in 1:size(mesh,2)
				new_arr[ind] = (mesh[nui, nupi, wi][3], mesh[nui, nupi, wi][1], mesh[nui, nupi, wi][2])
				ind += 1
			end
		end
	end
	return new_arr
end

function triqs_linearize(raw_arr)
	new_arr = Array{Complex{Float64},1}(undef, size(raw_arr,4)*size(raw_arr,3)*size(raw_arr,2))
	ind = 1
	for wi in 1:size(raw_arr,4)
		for nui in 1:size(raw_arr,2)
			for nupi in 1:size(raw_arr,3)
				new_arr[ind] = raw_arr[1,nui, nupi, wi] + 1im .* raw_arr[2,nui, nupi, wi]
				ind += 1
			end
		end
	end
	return new_arr
end

# ======================================================= #
#                                                         #
#                      GF Stuff                           #
#                                                         #
# ======================================================= #
@inline get_symm_f(f::Array{Complex{Float64},1}, i::Int64) = (i < 0) ? conj(f[-i]) : f[i+1]

function computeχ0(ω_range::AbstractArray{Int,1}, ν_range::AbstractArray{Int,1}, gImp::Array{Complex{Float64}, 1}, β::Float64)
    χ0 = Dict{Tuple{Int,Int},Complex{Float64}}()
    for ω in ω_range, ν in ν_range
        χ0[(ω,ν)] = -β*get_symm_f(gImp, ν)*get_symm_f(gImp, ν+ω)
    end
    return χ0
end

function computeΓ(freqList::Array, χ::Array{T,3}, χ0::Dict{Tuple{Int,Int},Complex{Float64}}, bGrid, fGrid) where T
    res = Array{T}(undef, length(bGrid), length(fGrid), length(fGrid))
    for (ωn,ω) in enumerate(bGrid)
        res[ωn,:,:] = inv(χ[ωn,:,:])
        for (νn,ν) in enumerate(fGrid)
            res[ωn,νn,νn] -= 1.0/χ0[(ω,ν)]
        end
    end
    return res
end

function gen_config(U::Float64, μ::Float64, β::Float64, nden::Float64, t::Float64)
cfg="""[Model]
U    = $U
mu   = $μ
beta = $β
nden = $nden
kGrid = "2Dsc-$t"

[Simulation]
Nk        = [4,20]                      # number of k points in each dimension
tail_correction = "Nothing"             # improvement algorithm for sum extrapolations. possible values: Nothing, Richardson, Shanks
lambda_correction = "sp"                # nothing, sp, sp_ch
bosonic_sum = "common"                  # common (intersection of individual ranges), individual (max range in each channel; fortran default), full (all frequencies), fixed:N:M (always sum from N to (including) M, indexing starts at 0)
force_full_bosonic_chi = true           # compute all omega frequencies for chi and trilex
chi_unusable_fill_value = "chi_lambda"  # can be "0", "chi_lambda" or "chi". sets either 0, lambda corrected or non lambda corrected values outside usable omega range
rhs  = "native"                         # native (fixed for tc, error_comp for naive), fixed (n/2 (1 - n/2) - sum(chi_ch)), error_comp (chi_loc_ch + chi_loc_sp - chi_ch)
fermionic_tail_coeffs = [0,1,2,3,4]     # internal parameter for richardson sum extrapolation
bosonic_tail_coeffs = [0,1,2,3,4]       # internal parameter for richardson sum extrapolation
usable_prct_reduction = 0.1             # safety cutoff for usable ranges, 0.0 means all values where chi is positive and strictly decreasing

[Environment]
inputDataType = "jld2"                  # LEGACY, do not change
writeFortran = false                    # LEGACY, do not change
loadAsymptotics = false                 # LEGACY, do not change
inputDir = ""                           # absolute path to input dir
freqFile = ""                           # absolute path to freqList.jld2 file
inputVars = "triqs_out.jld2"
asymptVars = "vars_asympt_sums.jld"     # LEGACY, do not change
cast_to_real = false                    # TODO: not implemented. cast all arrays with vanishing imaginary part to real
loglevel = "debug"                      # error, warn, info, debug
logfile = "stderr"                      # STDOUT, STDERR, filename
progressbar = false                     # LEGACY, do not change

[legacy]
lambda_correction = true                # Should a lambda-correction be performed only in the spin-channel?

[Debug]
read_bubble = false                     # LEGACY, do not change
"""
    return cfg
end


text:
# DBSCAN Clustering
#
#   References:
#
#       Martin Ester, Hans-peter Kriegel, Jörg S, and Xiaowei Xu
#       A density-based algorithm for discovering clusters
#       in large spatial databases with noise. 1996.
#

type DbscanResult <: ClusteringResult
    seeds::Vector{Int}          # starting points of clusters, size (k,)
    assignments::Vector{Int}    # assignments, size (n,)
    counts::Vector{Int}         # number of points in each cluster, size (k,)
end

## main algorithm

function dbscan{T<:Real}(D::DenseMatrix{T}, eps::Real, minpts::Int)
    # check arguments
    n = size(D, 1)
    size(D, 2) == n || error("D must be a square matrix.")
    n >= 2 || error("There must be at least two points.")
    eps > 0 || error("eps must be a positive real value.")
    minpts >= 1 || error("minpts must be a positive integer.")

    # invoke core algorithm
    _dbscan(D, convert(T, eps), minpts, 1:n)
end

function _dbscan{T<:Real}(D::DenseMatrix{T}, eps::T, minpts::Int, visitseq::AbstractVector{Int})
    n = size(D, 1)

    # prepare
    seeds = Int[]
    counts = Int[]
    assignments = zeros(Int, n)
    visited = zeros(Bool, n)
    k = 0

    # main loop
    for p in visitseq
        if assignments[p] == 0 && !visited[p]
            visited[p] = true
            nbs = _dbs_region_query(D, p, eps)
            if length(nbs) >= minpts
                k += 1
                cnt = _dbs_expand_cluster!(D, k, p, nbs, eps, minpts, assignments, visited)
                push!(seeds, p)
                push!(counts, cnt)
            end
        end
    end

    # make output
    return DbscanResult(seeds, assignments, counts)
end

## key steps

function _dbs_region_query{T<:Real}(D::DenseMatrix{T}, p::Int, eps::T)
    n = size(D,1)
    nbs = Int[]
    dists = view(D,:,p)
    for i = 1:n
        @inbounds if dists[i] < eps
            push!(nbs, i)
        end
    end
    return nbs::Vector{Int}
end

function _dbs_expand_cluster!{T<:Real}(D::DenseMatrix{T},           # distance matrix
                                       k::Int,                      # the index of current cluster
                                       p::Int,                      # the index of seeding point
                                       nbs::Vector{Int},            # eps-neighborhood of p
                                       eps::T,                      # radius of neighborhood
                                       minpts::Int,                 # minimum number of neighbors of a density point
                                       assignments::Vector{Int},    # assignment vector
                                       visited::Vector{Bool})       # visited indicators
    assignments[p] = k
    cnt = 1
    while !isempty(nbs)
        q = shift!(nbs)
        if !visited[q]
            visited[q] = true
            qnbs = _dbs_region_query(D, q, eps)
            if length(qnbs) > minpts
                for x in qnbs
                    if assignments[x] == 0
                        push!(nbs, x)
                    end
                end
            end
        end
        if assignments[q] == 0
            assignments[q] = k
            cnt += 1
        end
    end
    return cnt
end


text:
using ProgressMeter

"""
    likelihood(yhat, experimental_ecdf, experimental, I, ν0)

Find the likelihood of a specific set of parameters (yhat), which is based on
the ordinary least squares of the experimental vs. theoretical CDFs

# Arguments
- `yhat`: an array containing parameters for qcc, η, and the standard deviation
- `experimental_ecdf`: the estimated ECDF of the experimental data
- `experimental`: the experimental data
- `I`: spin (3 in the case of Boron-10)
- `ν0`: the Larmor frequency

"""
function likelihood(yhat, experimental_ecdf, experimental, I, ν0)
    # We must choose a distribution for the error!
    yhat[5] < 0 && return 0
    likelihood_dist = Normal(0, yhat[5])

    # Here we're generating the theoretical sample and converting it to a CDF
    quad = Quadrupolar(yhat[1:4])
    ismissing(quad) && return 0
    powder_pattern = estimate_powder_pattern(quad, 100_000, ν0, I)
    theoretical_ecdf = ecdf(powder_pattern).(experimental[:, 1])
    theoretical_ecdf .-= theoretical_ecdf[1]
    theoretical_ecdf ./= theoretical_ecdf[end]
    # Then we return the likelihood, based on two CDFs' differences
    return sum(logpdf.(likelihood_dist, experimental_ecdf .- theoretical_ecdf))
end

function likelihood(yhat, experimental_ecdf, experimental; sites = 1)
    if sites == 1
        # We must choose a distribution for the error!
        yhat[7] < 0 && return 0
        likelihood_dist = Normal(0, yhat[7])

        # Here we're generating the theoretical sample and converting it to a CDF
        csa = ChemicalShift(vcat(yhat[1:6], [1.0]))
        ismissing(csa) && return 0
        powder_pattern = estimate_powder_pattern(csa, 100_000)
        theoretical_ecdf = ecdf(powder_pattern).(experimental[:, 1])

        # Then we return the likelihood, based on two CDFs' differences
        return sum(logpdf.(likelihood_dist, experimental_ecdf .- theoretical_ecdf))
    elseif sites == 2
        # We must choose a distribution for the error!
        yhat[14] < 0 && return 0
        likelihood_dist = Normal(0, yhat[14])

        # Here we're generating the theoretical sample and converting it to a CDF
        csa = ChemicalShift(vcat(yhat[1:6], [yhat[13]], yhat[7:12], [1 - yhat[13]]))
        ismissing(csa) && return 0
        powder_pattern = estimate_powder_pattern(csa, 100_000)
        theoretical_ecdf = ecdf(powder_pattern).(experimental[:, 1])

        # Then we return the likelihood, based on two CDFs' differences
        return sum(logpdf.(likelihood_dist, experimental_ecdf .- theoretical_ecdf))
    end
end

"""
    metropolis_hastings(experimental)

An implementation of the Metropolis-Hastings algorithm for Monte Carlo Markov
Chain (MCMC) Bayesian estimates, with optional keyword arguments for the number
of samples (default 1,000,000) and tolerances for each parameter (default 0.1
for qcc, 0.1 for η, 0.05 for σ)

"""
function metropolis_hastings(
    experimental;
    interaction = "quadrupolar",
    N = 1_000_000,
    sites = 1,
    tol = interaction == "quadrupolar" ? [0.1, 0.2, 0.1, 0.1, 0.05] : sites == 1 ? [100.0, 20.0, 30.0, 5.0, 0.1, 0.05, 0.05] : [10.0, 5.0, 10.0, 5.0, 0.1, 0.05, 10.0, 5.0, 10.0, 5.0, 0.1, 0.05, 0.01, 0.05],
    I = 3
)
    experimental_ecdf = get_experimental_ecdf(experimental)
    if interaction == "quadrupolar"
        ν0 = get_ν0(experimental, experimental_ecdf)

        samples = zeros(N, 5)  # initialize zero array for samples

        #We need to define prior distributions for each parameter
        prior_dist_qcc = Uniform(0, 9)
        prior_dist_σqcc = Uniform(0, 3)
        prior_dist_η = Uniform(0, 1)
        prior_dist_ση = Uniform(0, 1)
        prior_dist_σ = Uniform(0, 1)

        prior_quad(x) = logpdf(prior_dist_qcc, x[1]) + logpdf(prior_dist_σqcc, x[2]) +
            logpdf(prior_dist_η, x[3]) + logpdf(prior_dist_ση, x[4]) +
            logpdf(prior_dist_σ, x[5]) # Assume variables are independent

        a = [rand(prior_dist_qcc), rand(prior_dist_σqcc), rand(prior_dist_η),
            rand(prior_dist_ση),rand(prior_dist_σ)]
        samples[1, :] = a

        # Repeat N times
        @showprogress for i = 2:N
            b = a + tol .* (rand(5) .- 0.5)  # Compute new state randomly
            # Calculate density
            prob_old = likelihood(a, experimental_ecdf, experimental, I, ν0) +
                       prior_quad(a)
            prob_new = likelihood(b, experimental_ecdf, experimental, I, ν0) +
                       prior_quad(b)
            r = prob_new - prob_old # Compute acceptance ratio
            if log(rand()) < r
                a = b  # Accept new state and update
            end
            samples[i, :] = a  # Update state
        end
    else
        if sites == 1
            samples = zeros(N, 7)  # initialize zero array for samples

            #We need to define prior distributions for each parameter
            prior_dist_σᵢₛₒ = Uniform(minimum(experimental[:, 1]), maximum(experimental[:, 1]))
            prior_dist_σσᵢₛₒ = Uniform(0.00001, 800)
            prior_dist_Δσ = Uniform(-4000, 4000)
            prior_dist_σΔσ = Uniform(0.00001, 400)
            prior_dist_ησ = Uniform(0, 1)
            prior_dist_σησ = truncated(Normal(0.00001, 0.2), 0.0, Inf)
            prior_dist_σ = truncated(Normal(0.00001, 0.5), 0.0, Inf)

            prior_csa(x) = logpdf(prior_dist_σᵢₛₒ, x[1]) + logpdf(prior_dist_σσᵢₛₒ, x[2]) +
                logpdf(prior_dist_Δσ, x[3]) + logpdf(prior_dist_σΔσ, x[4]) +
                logpdf(prior_dist_ησ, x[5]) + logpdf(prior_dist_σησ, x[6]) +
                logpdf(prior_dist_σ, x[7]) # Assume variables are independent

            a = [rand(prior_dist_σᵢₛₒ), rand(prior_dist_σσᵢₛₒ), rand(prior_dist_Δσ),
                rand(prior_dist_σΔσ), rand(prior_dist_ησ), rand(prior_dist_σησ),
                rand(prior_dist_σ)]
            samples[1, :] = a

            # Repeat N times
            @showprogress for i = 2:N
                b = a + tol .* (rand(7) .- 0.5)  # Compute new state randomly
                # Calculate density
                prob_old = likelihood(a, experimental_ecdf, experimental) + prior_csa(a)
                prob_new = likelihood(b, experimental_ecdf, experimental) + prior_csa(b)
                r = prob_new - prob_old # Compute acceptance ratio
                if log(rand()) < r
                    a = b  # Accept new state and update
                end
                samples[i, :] = a  # Update state
            end
        elseif sites == 2
            samples = zeros(N, 14)  # initialize zero array for samples

            #We need to define prior distributions for each parameter
            prior_dist_σᵢₛₒ1 = Uniform(minimum(experimental[:, 1]), maximum(experimental[:, 1]))
            prior_dist_σσᵢₛₒ1 = Uniform(0.00001, 800)
            prior_dist_Δσ1 = Uniform(-4000, 4000)
            prior_dist_σΔσ1 = Uniform(0.00001, 400)
            prior_dist_ησ1 = Uniform(0, 1)
            prior_dist_σησ1 = truncated(Normal(0.00001, 0.2), 0.0, Inf)
            prior_dist_σᵢₛₒ2 = Uniform(minimum(experimental[:, 1]), maximum(experimental[:, 1]))
            prior_dist_σσᵢₛₒ2 = Uniform(0.00001, 800)
            prior_dist_Δσ2 = Uniform(-4000, 4000)
            prior_dist_σΔσ2 = Uniform(0.00001, 400)
            prior_dist_ησ2 = Uniform(0, 1)
            prior_dist_σησ2 = truncated(Normal(0.00001, 0.2), 0.0, Inf)
            prior_dist_w = Uniform(0, 1)
            prior_dist_σ = truncated(Normal(0.00001, 0.5), 0.0, Inf)

            a = [rand(prior_dist_σᵢₛₒ1), rand(prior_dist_σσᵢₛₒ1), rand(prior_dist_Δσ1),
                rand(prior_dist_σΔσ1), rand(prior_dist_ησ1), rand(prior_dist_σησ1),
                rand(prior_dist_σᵢₛₒ2), rand(prior_dist_σσᵢₛₒ2), rand(prior_dist_Δσ2),
                rand(prior_dist_σΔσ2), rand(prior_dist_ησ2), rand(prior_dist_σησ2),
                rand(prior_dist_w),
                rand(prior_dist_σ)]
            samples[1, :] = a

            # Repeat N times
            @showprogress for i = 2:N
                b = a + tol .* (rand(14) .- 0.5)  # Compute new state randomly

                # Calculate density
                prob_old = likelihood(a, experimental_ecdf, experimental, sites = 2) +
                logpdf(prior_dist_σᵢₛₒ1, a[1]) + logpdf(prior_dist_σσᵢₛₒ1, a[2]) +
                   logpdf(prior_dist_Δσ1, a[3]) + logpdf(prior_dist_σΔσ1, a[4]) +
                   logpdf(prior_dist_ησ1, a[5]) + logpdf(prior_dist_σησ1, a[6]) +
                   logpdf(prior_dist_σᵢₛₒ2, a[7]) + logpdf(prior_dist_σσᵢₛₒ2, a[8]) +
                   logpdf(prior_dist_Δσ2, a[9]) + logpdf(prior_dist_σΔσ2, a[10]) +
                   logpdf(prior_dist_ησ2, a[11]) + logpdf(prior_dist_σησ2, a[12]) +
                   logpdf(prior_dist_w, a[13]) +
                   logpdf(prior_dist_σ, a[14])
                prob_new = likelihood(b, experimental_ecdf, experimental, sites = 2) +
                logpdf(prior_dist_σᵢₛₒ1, b[1]) + logpdf(prior_dist_σσᵢₛₒ1, b[2]) +
                   logpdf(prior_dist_Δσ1, b[3]) + logpdf(prior_dist_σΔσ1, b[4]) +
                   logpdf(prior_dist_ησ1, b[5]) + logpdf(prior_dist_σησ1, b[6]) +
                   logpdf(prior_dist_σᵢₛₒ2, b[7]) + logpdf(prior_dist_σσᵢₛₒ2, b[8]) +
                   logpdf(prior_dist_Δσ2, b[9]) + logpdf(prior_dist_σΔσ2, b[10]) +
                   logpdf(prior_dist_ησ2, b[11]) + logpdf(prior_dist_σησ2, b[12]) +
                   logpdf(prior_dist_w, b[13]) +
                   logpdf(prior_dist_σ, b[14])
                r = prob_new - prob_old # Compute acceptance ratio
                if log(rand()) < r
                    a = b  # Accept new state and update
                end
                samples[i, :] = a  # Update state
            end
        end
    end

    return samples
end;


text:
# Boundary conditions for OceanTurb.jl

import Base: convert, setproperty!

abstract type BCType end
struct Flux <: BCType end
struct Gradient <: BCType end
struct Value <: BCType end

"Return c₁, where c₁ = c₀ + ∂c/∂z * (z₁ - z₀) = c₀ + ∂c/∂z * Δz."
@inline linear_extrap(c₀, ∂c∂z, Δz) = c₀ + ∂c∂z * Δz

mutable struct BoundaryCondition{C, T}
    condition :: T
end

function Base.convert(::Type{BoundaryCondition{C, T1}}, bc::BoundaryCondition{C}) where {C, T1}
    return BoundaryCondition{C, T1}(convert(T1, bc.condition))
end

BoundaryCondition(bctype, bc) = BoundaryCondition{bctype, typeof(bc)}(bc)

GradientBoundaryCondition(bc) = BoundaryCondition(Gradient, bc)
FluxBoundaryCondition(bc) = BoundaryCondition(Flux, bc)
ValueBoundaryCondition(bc) = BoundaryCondition(Value, bc)

const BC = BoundaryCondition

@inline getbc(model, bc::BC{C, <:Number}) where C = bc.condition
@inline getbc(model, bc::BC{C, <:Function}) where C = bc.condition(model)

@inline gradient(bc::BC{<:Gradient}, model, args...) = getbc(model, bc)
@inline gradient(bc::BC{<:Flux},     model, κ, args...) = -getbc(model, bc) / κ
@inline gradient(bc::BC{<:Value},    model, κ, cᴺ, Δf, args...) = 2 * (getbc(model, bc) - cᴺ) / Δf

"""
    fill_bottom_ghost_cell!(c, κ, model, bc)

Update the bottom ghost cell of c given the boundary condition `bc`, `model`, and
diffusivity `kappa`. `kappa` is used only if a flux boundary condition
is specified.
"""
@inline function fill_bottom_ghost_cell!(bc, c, κ, model)
    @inbounds c[0] = linear_extrap(c[1], gradient(bc, model, κ, c[1], -Δf(c, 1)), -Δc(c, 1))
    return nothing
end

"""
    fill_top_ghost_cell!(c, κ, model, bc)

Update the top ghost cell of `c` given boundary condition `bc`, `model`, and
diffusivity `kappa`
"""
@inline function fill_top_ghost_cell!(bc, c, κ, model, N)
    @inbounds c[N+1] = linear_extrap(c[N], gradient(bc, model, κ, c[N], Δf(c, N)), Δc(c, N+1))
    return nothing
end

"""
    FieldBoundaryConditions(; top=TopBC, bottom=BottomBC)

Create an instance of `FieldBoundaryConditions` with `top` and `bottom`
"""
mutable struct FieldBoundaryConditions{B, T}
    bottom :: B
    top    :: T
end

"""
    ZeroFluxBoundaryConditions(T=Float64)

Construct `FieldBoundaryConditions` with a zero `FluxBoundaryCondition`
at top and bottom.
"""
function ZeroFluxBoundaryConditions(T=Float64)
    FieldBoundaryConditions(
        FluxBoundaryCondition(-zero(T)),
        FluxBoundaryCondition(-zero(T))
    )
end

"""
    DefaultBoundaryConditions(T=Float64)

Returns default oceanic boundary conditions: a zero
`GradientBoundaryCondition` on bottom and a zero
`FluxBoundaryCondition` on top.
"""
function DefaultBoundaryConditions(T=Float64)
    return FieldBoundaryConditions(
                GradientBoundaryCondition(-zero(T)),
                FluxBoundaryCondition(-zero(T))
                )
end

"""
    BoundaryConditions([T=Float64;] bottom = GradientBoundaryCondition(-zero(T)),
                                       top = FluxBoundaryCondition(-zero(T)))

Returns `FieldBoundaryConditions` with a `bottom` and `top` boundary condition.
The type `T` is only relevant for the default values of `bottom` and `top`.
"""
function BoundaryConditions(T=Float64; bottom = GradientBoundaryCondition(-zero(T)),
                                          top = FluxBoundaryCondition(-zero(T))
                           )
    return FieldBoundaryConditions(bottom, top)
end

"""
    set_bcs!(model; bcspecs...)

Set boundary conditions of model solution fields.
The keyword argument name must be the name of a model solution
and its value is a (bottom_bc, top_bc) tuple.

Example
========

julia> set_bcs!(model, c=(FluxBoundaryCondition(-1),
                          FluxBoundaryCondition(0))
                )
"""
function set_bcs!(model; bcspecs...)
    bcs = model.bcs
    for (ϕsym, bcs) in bcspecs
        ϕbcs = getproperty(model.bcs, ϕsym)
        setproperty!(ϕbcs, :bottom, bcs[1])
        setproperty!(ϕbcs, :top, bcs[2])
    end
    return nothing
end