text:
using Test
using Distributions
using KernelDensity

import KernelDensity: kernel_dist, default_bandwidth, kde_boundary, kde_range, tabulate

for D in [Normal,Uniform,Logistic]
    d = kernel_dist(D,0.5)
    @test isa(d,D)
    @test mean(d) == 0.0
    @test std(d) ≈ 0.5
end

r = kde_range((-2.0,2.0), 128)
@test step(r) > 0
r2 = kde_range((0.12698109160784082, 0.9785547869337731), 256)
@test length(r2) == 256

for X in ([0.0], [0.0,0.0], [0.0,0.5], [-0.5:0.1:0.5;])
    w = default_bandwidth(X)
    @test w > 0
    lo, hi = kde_boundary(X,w)
    @test lo < hi
    kr = kde_range((lo,hi), 10)
    @test step(kr) > 0

    for D in (Normal, )
        k1 = tabulate(X,r)
        @test isa(k1,UnivariateKDE)
        @test length(k1.density) == length(k1.x)
        @test all(k1.density .>= 0.0)
        @test sum(k1.density)*step(k1.x) ≈ 1.0

        k2 = KernelDensity.conv(k1,kernel_dist(D,0.1))
        @test isa(k2,UnivariateKDE)
        @test length(k2.density) == length(k2.x)
        @test all(k2.density .>= 0.0)
        @test sum(k2.density)*step(k2.x) ≈ 1.0

        k3 = kde(X;kernel=D)
        @test isa(k3,UnivariateKDE)
        @test length(k3.density) == length(k3.x)
        @test all(k3.density .>= 0.0)
        @test sum(k3.density)*step(k3.x) ≈ 1.0

        k4 = kde(X,r;kernel=D)
        @test isa(k4,UnivariateKDE)
        @test length(k4.density) == length(k4.x)
        @test all(k4.density .>= 0.0)
        @test sum(k4.density)*step(k4.x) ≈ 1.0

        k5 = kde_lscv(X)
        @test isa(k5,UnivariateKDE)
        @test length(k5.density) == length(k5.x)
        @test all(k5.density .>= 0.0)
        @test sum(k5.density)*step(k5.x) ≈ 1.0

        k6 = kde(X,r;kernel=D, weights=fill(1.0/length(X),length(X)))
        @test k4.density ≈ k6.density
    end
end

k11 = kde([0.0, 1.], r, bandwidth=1, weights=[0,1])
k12 = kde([1.], r, bandwidth=1)
@test k11.density ≈ k12.density


text:
using Documenter
using LPSjl
using Graphs

makedocs(
    sitename = "Julia Bayesian Networks Library",
    format = Documenter.HTML(),
    modules = [LPSjl, LPSjl.BayesianNetworks, LPSjl.ParameterLearning],
    pages = [
        "index.md",
        "structureLearning.md",
        "Examples" => ["EX_BayesianNetworks.md"],
        "API" => ["API_BayesianNetworks.md", "API_ParameterLearning.md"]
    ]
)

deploydocs(
    repo = "github.com/vbob/LPSjl.jl.git"
)

ext: jl


text:
module HyperCubicRoots

# # Obtain the Gauss hypergeometric function ₂F₁(a, b; c; z):
import HypergeometricFunctions
HgF = HypergeometricFunctions._₂F₁

export HgF,
       solve_real_cubic_roots, solve_all_cubic_roots,
       solve_quadratic_roots,
       solve_all_quartic_roots

# Tolerance to maintain accuracy on the evaluation of the root (p(root) ≈ 0)
const CUBIC_ATOL = 10.0*sqrt(eps(Float64))
const QUARTIC_ATOL = 10.0*sqrt(eps(Float64))

poor_conditioning(v, compare_tol) = any(abs.(v) .> compare_tol) ? true : false

include("quadratic.jl")
include("cubic.jl")
include("quartic.jl")

end

max_stars_repo_name: jd-foster/HyperCubicRoots

max_issues_repo_name: jd-foster/HyperCubicRoots

max_forks_repo_name: jd-foster/HyperCubicRoots


text:
x = 6
if x < 10
    println("Too small")
end

y = 15
if y < 10
    println("y is too small")
elseif y > 10 && y < 25
    println("y is medium")
else
    println("y is big")
end

z = 7
println(z < 10 ? "z is small" : "z is big")

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
export GramMatrix

import MultivariateMoments: trimat, SymMatrix, getmat
export getmat

"""
    struct GramMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: MP.APL{T}
        Q::SymMatrix{T}
        x::MVT
    end

Gram matrix ``x^\\top Q x`` where `Q` is a symmetric matrix indexed by the
vector of monomials `x`.
"""
struct GramMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: MP.APL{T} # should be MP.APL{eltype(T)} but it doesn't work
    Q::SymMatrix{T}
    x::MVT
end
# When taking the promotion of a GramMatrix of JuMP.Variable with a Polynomial JuMP.Variable, it should be a Polynomial of AffExpr
MP.coefficienttype(::Type{<:GramMatrix{T}}) where {T} = Base.promote_op(+, T, T)
function MP.constantmonomial(p::GramMatrix)
    if isempty(p.x)
        return MP.constantmonomial(MP.monomialtype(p))
    else
        return MP.constantmonomial(p.x[1])
    end
end
function MP.monomialtype(::Union{GramMatrix{T, MT},
                                 Type{<:GramMatrix{T, MT}}}) where {T, MT}
    return MT
end
MP.polynomialtype(::Type{GramMatrix{T, MT, MVT}}) where {T, MT, MVT} = MP.polynomialtype(MT, MP.coefficienttype(GramMatrix{T, MT, MVT}))
MP.polynomialtype(::Type{GramMatrix{T, MT, MVT}}, ::Type{S}) where {S, T, MT, MVT} = MP.polynomialtype(MT, S)
MP.variables(p::GramMatrix) = MP.variables(p.x)
MP.nvariables(p::GramMatrix) = MP.nvariables(p.x)

Base.zero(::Type{GramMatrix{T, MT, MVT}}) where {T, MT, MVT} = GramMatrix{T, MT, MP.monovectype(MT)}(SymMatrix{T}(T[], 0), MP.emptymonovec(MT))
Base.iszero(p::GramMatrix) = iszero(MP.polynomial(p))

Base.getindex(p::GramMatrix, I...) = getindex(p.Q, I...)
Base.copy(p::GramMatrix) = GramMatrix(copy(p.Q), copy(p.x))

MultivariateMoments.getmat(p::GramMatrix{T}) where {T} = p.Q

function GramMatrix{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:MP.AbstractMonomial}
    GramMatrix{T, MT, MP.monovectype(x)}(trimat(T, f, length(x), σ), x)
end
GramMatrix{T}(f::Function, x::AbstractVector, σ) where T = GramMatrix{T}(f, MP.monomial.(x), σ)
function GramMatrix{T}(f::Function, x::AbstractVector) where T
    σ, X = MP.sortmonovec(x)
    GramMatrix{T}(f, X, σ)
end
GramMatrix(f::Function, x) = GramMatrix{Base.promote_op(f, Int, Int)}(f, x)

function GramMatrix(Q::AbstractMatrix{T}, x, σ) where T
    GramMatrix{T}((i,j) -> Q[σ[i], σ[j]], x)
end
function GramMatrix(Q::AbstractMatrix, x)
    σ, X = MP.sortmonovec(x)
    GramMatrix(Q, X, σ)
end

#function Base.convert{T, PT <: AbstractPolynomial{T}}(::Type{PT}, p::GramMatrix)
#    # coefficienttype(p) may be different than T and MP.polynomial(p) may be different than PT (different module)
#    convert(PT, MP.polynomial(p))
#end
function MP.polynomial(p::GramMatrix)
    MP.polynomial(getmat(p), p.x)
end
function MP.polynomial(p::GramMatrix, ::Type{S}) where {S}
    MP.polynomial(getmat(p), p.x, S)
end

Base.:(+)(x::MP.APL, y::GramMatrix) = x + MP.polynomial(y)
Base.:(+)(x::GramMatrix, y::MP.APL) = MP.polynomial(x) + y
Base.:(+)(x::GramMatrix, y::GramMatrix) = MP.polynomial(x) + MP.polynomial(y)
Base.:(-)(x::MP.APL, y::GramMatrix) = x - MP.polynomial(y)
Base.:(-)(x::GramMatrix, y::MP.APL) = MP.polynomial(x) - y
Base.:(-)(x::GramMatrix, y::GramMatrix) = MP.polynomial(x) - MP.polynomial(y)
Base.:(==)(p::MP.APL, q::GramMatrix) = p == MP.polynomial(q)
Base.:(==)(p::GramMatrix, q::GramMatrix) = iszero(p - q)
function Base.isapprox(p::GramMatrix, q::GramMatrix; kwargs...)
    p.x == q.x && isapprox(p.Q, q.Q; kwargs...)
end


text:
# User settings.
  # Time and space grid.
    z_min = 0.0
    z_max = 5.0
    P = 75
    T_default = 20.0
  # Experiment settings.
    d_0 = 5.0
    d_T = 2.3701
    Δ_E = 1e-06
  # Overall parameters.
    params = (ρ = 0.02, σ = 4.2508, N = 10, θ = 5.1269, γ = 1.00, κ = 0.013, ζ = 1, η = 0, Theta = 1, χ = 1/(2.1868), υ = 0.0593, μ = 0, δ = 0.053)
  # Solver settings.

# Construct intermediate objects.
  z_ex = range(z_min, stop = z_max, length = P+2)
  z = z_ex[2:end-1]
  settings = settings_defaults_tests(z_ex = z_ex, T = T_default)
  params_0 = merge(params, (d = d_0,))
  params_T = merge(params, (d = d_T,))

# Get (numerical) stationary solution.
  stationary_0 = stationary_numerical(params_0, settings_defaults_tests(z_ex = z_ex))
  stationary_T = stationary_numerical(params_T, settings_defaults_tests(z_ex = z_ex))

# Process those
  Ω_0 = stationary_0.Ω
  Ω_T = stationary_T.Ω
  v_T = stationary_T.v_tilde
  g_T = stationary_T.g
  z_hat_T = stationary_T.z_hat

# Define more interim quantities.
  Ω(t) = Ω_T # This is constant.

@testset "Regression Tests with constant Ω" begin
  # Run the solver.
    sol = solve_dynamics(params_T, stationary_T, settings, Ω)
    sol = merge(sol, (results = prepare_results(sol, stationary_T, stationary_0),))
    filter!(row -> row.t >= 0, sol.results) # get rid of pre-shock values

    # Spot-checks.
    @test sol.sol.t[5] ≈ 19.996000000000002 atol = 1e-5
    @test sol.results[!, :λ_ii][end] ≈ 0.8205696009690461 atol = 1e-5
    @test sol.sol.u[4][3] ≈ 0.713179042868647 atol = 1e-5
    @test sol.sol.prob.u0[1] ≈ 0.9329809580656181 atol = 1e-5
  # Consistency checks.
    @test mean(sol.results[!, :g]) ≈ 0.006332516508877774  atol = 1e-5 # Probably the most important of these checks.
    @test mean(sol.results[!, :z_hat]) ≈ 1.626173967797714 atol = 1e-5
    @test mean(sol.results[!, :Ω]) ≈ 1.6928323073432825 atol = 1e-5
    # @test all(sol.results[:E] .== 0.053)
    @test mean(sol.results[!, :S]) ≈ -0.013749473107398985 atol = 1e-5
    @test mean(sol.results[!, :L_tilde]) ≈ 0.1892988756097014 atol = 1e-5
  # Full coverage of each column.
    @test sol.results[!, :t][4] ≈ 0.75 atol = 1e-5
    @test sol.results[!, :g][5] ≈ 0.006332516539866723 atol = 1e-5
    @test sol.results[!, :z_hat][6] ≈ 1.6261739679041103 atol = 1e-5
    @test sol.results[!, :Ω][7] ≈ 1.6928323073432825 atol = 1e-5
    @test sol.results[!, :E][8] ≈ 0.05300001047668845 atol = 1e-5
    @test sol.results[!, :v_1][9] ≈ 0.9329809601676073 atol = 1e-5
    @test sol.results[!, :L_tilde][10] ≈ 0.18929887587415023 atol = 1e-5
    @test sol.results[!, :λ_ii][11] ≈ 0.8205696004618322 atol = 1e-5
    @test sol.results[!, :c][12] ≈ 1.3801007368941405 atol = 1e-5
    @test sol.results[!, :S][4] ≈ -0.013749472948335131 atol = 1e-5
    @test sol.results[!, :z_bar][3] ≈ 1.702354537311347 atol = 1e-5
    @test sol.results[!, :π_min][2] ≈ 0.04423572434773227 atol = 1e-5
    @test sol.results[!, :entry_residual][12] ≈ 2.7166018323754315e-9 atol = 1e-5
    @test sol.results[!, :r][19] ≈ 0.0793325164566891 atol = 1e-5
end

@testset "Regression Tests with time-varying Ω" begin
    # Run the solver.
    T = sqrt(2*(log(Ω_0) - log(Ω_T)) / params.δ)
    Ω_t(t) = t < T ? Ω_0 * exp(-params.δ*T*t + params.δ*t*t/2) : Ω_T # Exponential Ω with time smoothing
    settings = settings_defaults_tests(z_ex = z_ex, T = T)
    sol = solve_dynamics(params_T, stationary_T, settings, Ω_t)
    sol = merge(sol, (results = prepare_results(sol, stationary_T, stationary_0),))
    filter!(row -> row.t >= 0, sol.results)

    # Spot-checks.
    @test sol.sol.t[5] ≈ 3.77533130102673 atol = 1e-5
    @test sol.results[!, :λ_ii][end] ≈ 0.8205696009690461 atol = 1e-5
    @test sol.sol.u[4][3] ≈ 0.713186458707734 atol = 1e-5
    @test sol.sol.prob.u0[1] ≈ 0.9329809580656181 atol = 1e-5
    # Consistency checks.
    @test mean(sol.results[!, :g]) ≈ 0.005968185219011899 atol = 1e-5 # Probably the most important of these checks.
    @test mean(sol.results[!, :z_hat]) ≈ 1.6551587501780332 atol = 1e-5
    @test mean(sol.results[!, :Ω]) ≈ 1.8166904089164235 atol = 1e-5
    @test mean(sol.results[!, :S]) ≈ -0.01561736319741233 atol = 1e-5
    @test mean(sol.results[!, :L_tilde]) ≈ 0.18682781261445597 atol = 1e-5
    # Full coverage of each column.
    @test sol.results[!, :t][4] ≈ 0.75 atol = 1e-5
    @test sol.results[!, :g][5] ≈ 0.005164669443013719 atol = 1e-5
    @test sol.results[!, :z_hat][6] ≈ 1.7014106198547416 atol = 1e-5
    @test sol.results[!, :Ω][7] ≈ 1.9418121379049862 atol = 1e-5
    @test sol.results[!, :v_1][9] ≈ 0.9329809578947367 atol = 1e-5
    @test sol.results[!, :L_tilde][10] ≈ 0.18787808053459026 atol = 1e-5
    @test sol.results[!, :λ_ii][11] ≈ 0.8237900176167677 atol = 1e-5
    @test sol.results[!, :c][12] ≈ 1.3916274048069284 atol = 1e-5
    @test sol.results[!, :S][4] ≈ -0.020681135472245087 atol = 1e-5
    @test sol.results[!, :z_bar][3] ≈ 1.844294549955518 atol = 1e-5
    @test sol.results[!, :π_min][2] ≈ 0.03330744238041946 atol = 1e-5
    @test sol.results[!, :entry_residual][12] ≈ 0.0 atol = 1e-5
    @test sol.results[!, :r][19] ≈ 0.07705459878582416 atol = 1e-5
    # Sub-pieces of L_tilde
    @test sol.results[!, :L_tilde_a] + sol.results[!, :L_tilde_x] + sol.results[!, :L_tilde_E] ≈ sol.results[!, :L_tilde]
end

@testset "Correctness Tests" begin # Here, we compare the DAE output to known correct values, such as MATLAB output or analytical results.
  # First case.
  sol = solve_dynamics(params_T, stationary_T, settings, Ω)
  sol = merge(sol, (results = prepare_results(sol, stationary_T, stationary_0),))
  filter!(row -> row.t >= 0, sol.results)
  @test all([isapprox(x, 0.0, atol = 1e-6) for x in sol.results[!, :entry_residual]]) # Free-entry condition holds ∀ t.
end

@testset "Regression Tests for f! at T" begin
  # Instantiate the f! arguments
    du = zeros(P+3)
    u = [v_T..., g_T, z_hat_T, params.δ]
    resid = zeros(P+3)
    t = T_default
    p = solve_dynamics(params_T, stationary_T, settings, Ω).dae_parameters
  # Fill the residuals vector
    f!(resid, du, u, p, t)
  # Tests
    # Accuracy
      @test mean(resid[1:P]) ≈ 0 atol = 1e-6
      @test mean(resid[P+1]) ≈ 0 atol = 1e-6
      @test mean(resid[P+2]) ≈ 0 atol = 1e-6
      @test mean(resid[P+3]) ≈ 0 atol = 1e-6
      @test all(abs.(resid) .<= 1e-6) # Test that we have small residuals across the board.
    # Regression
      @test resid[4] ≈ 0 atol = 1e-6
      @test resid[end] ≈ 0 atol = 1e-6
    # Benchmarking results
      println("Kernel Benchmarking")
      @show @btime f!($resid, $du, $u, $p, $t)
end


text:
module test_pkgs_macrotools

using Base.Test
using MacroTools # @capture
using MacroTools: prewalk, postwalk

@capture :[1, 2, 3, 4, 5, 6, 7]  [1, a_, 3, b__, c_]
@test (a, b, c) == (2, Any[4, 5, 6], 7)

end  # module test_pkgs_macrotools


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Non Linear Burger Equation

using NeuralPDE, Flux, ModelingToolkit, Optim, DiffEqFlux
using GalacticOptim
using Base
using Plots, PyPlot
#GalacticOptim
# 3D PDE
@parameters x t θ
@variables u(..)
@derivatives Dxx''~x
#@derivatives Dyy''~y
@derivatives Dx'~x
#@derivatives Dy'~y
@derivatives Dt'~t

Re = 2500.0

# 3D PDE
eq  = Dt(u(x,t,θ)) + u(x,t,θ) * Dx(u(x,t,θ)) ~ 1/Re * Dxx(u(x,t,θ))
#eq  = Dt(u(x,y,t,θ)) + u(x,y,t,θ) * Dx(u(x,y,t,θ)) ~ 1/Re * (Dxx(u(x,y,t,θ)) + Dyy(u(x,y,t,θ)))

# Initial and boundary conditions

A0   = 1.0
muX0 = 0.
muY0 = 0.
sgX0 = 1.0
sgY0 = 1.0

bcs = [u(x,0,θ) ~ A0 * exp(-(((x-muX0)^2/(2*sgX0^2))))]#+((y-muY0)^2/(2*sgY0^2))))]

#bcs = [u(x,y,0,θ) ~ exp(x+y)*cos(x+y) ,
#       u(0,y,t,θ) ~ exp(y)*cos(y+4t),
#       u(2,y,t,θ) ~ exp(2+y)*cos(2+y+4t) ,
#       u(x,0,t,θ) ~ exp(x)*cos(x+4t),
#       u(x,2,t,θ) ~ exp(x+2)*cos(x+2+4t)]

       # Space and time domains
domains = [x ∈ IntervalDomain(-5.0,5.0),
           t ∈ IntervalDomain(-5.0,5.0)]

# Discretization
dx = 0.5; dt = 0.5

# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = NeuralPDE.PhysicsInformedNN([dx,dt],
                                             chain,
                                             strategy = NeuralPDE.StochasticTraining(include_frac=0.9))
pde_system = PDESystem(eq,bcs,domains,[x,t],[u])
prob = NeuralPDE.discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, ADAM(0.1), progress = true; cb = cb, maxiters=500)
phi = discretization.phi

xs,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

#xs,ys,ts = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

#u_predict = [reshape([first(phi([x,y,t],res.minimizer)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts]

analytic_sol_func(x,t) = [exp(-t)*sin(pi*x), exp(-t)*cos(pi*x), (1+pi^2)*exp(-t)]
#u_real  = [reshape([first(analytic_sol_func(x,t)[i] for x in xs for y in ys], (length(xs),length(ys))) for t in ts]
u_real  = [[analytic_sol_func(x,t)[i] for t in ts for x in xs] for i in 1:3]

u_predict  = [[phi([x,t],res.minimizer)[i] for t in ts for x in xs] for i in 1:1]
u_predict = u_predict[1]
#u_predict = [reshape([first(phi([x,t],res.minimizer)) for x in xs for t in ts])]

diff_u = [abs.(u_real[i] .- u_predict[i] ) for i in 1:3]

p2 = Plots.plot(xs, ts, u_predict, st=:surface,title = "predict")
Plots.plot(p2)

for i in 1:3
    p1 = Plots.plot(xs, ts, u_real[i], st=:surface,title = "u$i, analytic");
    p2 = Plots.plot(xs, ts, u_predict[i], st=:surface,title = "predict");
    #p3 = Plots.plot(xs, ts, diff_u[i],linetype=:contourf,title = "error");
    #plot(p1,p2,p3)
    Plots.plot(p2)
    Plots.savefig("sol_U$i")
end

maxlim = maximum(maximum(u_predict[t]) for t = 1:length(ts))
minlim = minimum(minimum(u_predict[t]) for t = 1:length(ts))

result = @animate for time = 1:length(ts)
    Plots.plot(xs, ys, u_predict[time],st=:surface,camera=(30,30), zlim=(minlim,maxlim), clim=(minlim,maxlim),
                title = string("ψ: max = ",round(maxlim, digits = 3)," min = ", round(minlim, digits = 3),"\\n"," t = ",time))
end

gif(result, "burger_surface.gif", fps = 6)

result = @animate for time = 1:length(ts)
    Plots.plot(xs, ys, u_predict[time],st=:contour,camera=(30,30), zlim=(minlim,maxlim), clim=(minlim,maxlim),
                title = string("ψ: max = ",round(maxlim, digits = 3)," min = ", round(minlim, digits = 3),"\\n"," t = ",time))
end

gif(result, "burger_contour.gif", fps = 6)


text:
include("setup_1d.jl")

println("idwt by lifting (N=",N,"), ", L, " levels")
f(x0, L, wl) = for i = 1:tn; idwt(x0, wl, L); end
f(x0, L, wl);
@time f(x0, L, wl);

#@profile f(x0,L,wl);
#Profile.print()


text:
using MolecularGraph: ATOMTABLE, ATOMSYMBOLMAP, ATOM_COVALENT_RADII, ATOM_VANDERWAALS_RADII

ATOMSYMBOLKEYS = ATOMSYMBOLMAP |> collectkeys
ATOMSYMBOLVALS = ATOMSYMBOLMAP |> collectvals
sp = sortperm(ATOMSYMBOLMAP |> collectvals)
atomicmasses = OrderedDict(ATOMSYMBOLKEYS[sp].=>[ATOMTABLE[ATOMSYMBOLMAP[x]]["Weight"] for x in ATOMSYMBOLKEYS[sp]])

_covrad = []
for x in ATOMSYMBOLKEYS[sp]
	try
		push!(_covrad,OrderedDict(x=>ATOM_COVALENT_RADII[ATOMTABLE[ATOMSYMBOLMAP[x]]["Number"]]))
	catch

	end
end
covrad = merge(_covrad...)

_vdwrad = []
for x in ATOMSYMBOLKEYS[sp]
	try
		push!(_vdwrad,OrderedDict(x=>ATOM_VANDERWAALS_RADII[ATOMTABLE[ATOMSYMBOLMAP[x]]["Number"]]))
	catch

	end
end
vdwrad = merge(_vdwrad...)

res3letters = ["ARG", "MET", "ASN", "GLU", "PHE",
	"ILE", "ASP", "LEU", "ALA", "GLN",
	"GLY", "CYS", "TRP", "TYR", "LYS",
	"PRO", "THR", "SER", "VAL", "HIS",
	"ASX", "GLX", "XAA", "XLE", "-", ".", "*"]
resletters = ["R", "M", "N", "E", "F",
	"I", "D", "L", "A", "Q",
	"G", "C", "W", "Y", "K",
	"P", "T", "S", "V", "H",
	"B", "Z", "X", "J", "-", ".", "*"]
resletterdict = OrderedDict(
	"ALA" => "A",
	"ARG" => "R",
	"ASN" => "N",
	"ASP" => "D",
	"ASX" => "B",
	"CYS" => "C",
	"GLN" => "Q",
	"GLU" => "E",
	"GLX" => "Z",
	"GLY" => "G",
	"HIS" => "H",
	"ILE" => "I",
	"LEU" => "L",
	"LYS" => "K",
	"MET" => "M",
	"PHE" => "F",
	"PRO" => "P",
	"SER" => "S",
	"THR" => "T",
	"TRP" => "W",
	"TYR" => "Y",
	"VAL" => "V",
	"XAA" => "X",
	"XLE" => "J",
	"*" => "*",
	"-" => "-",
	"." => ".",
	"A" => "ALA",
	"R" => "ARG",
	"N" => "ASN",
	"D" => "ASP",
	"B" => "ASX",
	"C" => "CYS",
	"Q" => "GLN",
	"E" => "GLU",
	"Z" => "GLX",
	"G" => "GLY",
	"H" => "HIS",
	"I" => "ILE",
	"L" => "LEU",
	"K" => "LYS",
	"M" => "MET",
	"F" => "PHE",
	"P" => "PRO",
	"S" => "SER",
	"T" => "THR",
	"W" => "TRP",
	"Y" => "TYR",
	"V" => "VAL",
	"X" => "XAA",
	"J" => "XLE"
)
heavyresbonds = Dict(
                "ARG" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD"],["CD","NE"],["NE","CZ"],["CZ","NH1"],["CZ","NH2"],["C","OXT"]],
                "MET" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","SD"],["SD","CE"],["C","OXT"]],
                "ASN" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","OD1"],["CG","ND2"],["C","OXT"]],
                "GLU" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD"],["CD","OE1"],["CD","OE2"],["C","OXT"]],
                "PHE" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD1"],["CG","CD1"],["CD1","CE1"],["CD2","CE2"],["CE1","CZ"],["CE2","CZ"],["C","OXT"]],
                "ILE" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG1"],
						["CB","CG2"],["CG1","CD1"],["C","OXT"]],
                "ASP" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","OD1"],["CG","OD2"],["C","OXT"]],
                "LEU" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD1"],["CG","CD2"],["C","OXT"]],
                "ALA" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["C","OXT"]],
                "GLN" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD"],["CD","OE1"],["CD","NE2"],["C","OXT"]],
                "GLY" => [["C","O"],["C","CA"],["CA","N"],["C","OXT"]],
                "CYS" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","SG"],["C","OXT"]],
                "TRP" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD1"],["CG","CD2"],["CD1","NE1"],["CE2","NE1"],["CD2","CE2"],["CD2","CE3"],
						["CE2","CZ2"],["CE3","CZ3"],["CZ2","CH2"],["CZ3","CH2"],["C","OXT"]],
                "TYR" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD1"],["CG","CD2"],["CD1","CE1"],["CD2","CE2"],["CE1","CZ"],
						["CE2","CZ"],["CZ","OH"],["C","OXT"]],
                "LYS" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","CD"],["CD","CE"],["CE","NZ"],["C","OXT"]],
                "PRO" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["N","CD"],["CB","CG"],["CG","CD"],["C","OXT"]],
                "THR" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","OG1"],["CB","CG2"],["C","OXT"]],
                "SER" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","OG"],["C","OXT"]],
				"VAL" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG1"],["CB","CG2"],["C","OXT"]],
				"NAG" => [],
				"HOH" => [],
                "HIS" => [["C","O"],["C","CA"],["CA","N"],["CA","CB"],["CB","CG"],
						["CG","ND1"],["CG","CD2"],["ND1","CE1"],["CD2","NE2"],["NE2","CE1"],["C","OXT"]],
				" ZN" => []
)
hresbonds = Dict(
                "ARG" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CA","HA2"],["CB","HB3"],["CB","HB2"],
						["CG","HG2"],["CG","HG3"],["CD","HD2"],["CD","HD3"],["NE","HE"],
						["NH1","HH11"],["NH1","HH12"],["NH2","HH21"],["NH2","HH22"]],
                "MET" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CA","HA2"],["CB","HB3"],["CB","HB2"],
						["CG","HG1"],["CG","HG2"],["CE","HE1"],["CE","HE2"],["CE","HE3"]],
                "ASN" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CA","HA2"],["CB","HB3"],["CB","HB2"],
						["ND2","HD21"],["ND2","HD22"]],
                "GLU" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CA","HA2"],["CB","HB3"],["CB","HB2"],
						["CG","HG3"],["CG","HG2"]],
                "PHE" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CD1","HD1"],
						["CD2","HD2"],["CE1","HE1"],["CE2","HE2"],["CZ","HZ"]],
                "ILE" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB"],["CG1","HG13"],["CG1","HG12"],
						["CG2","HG21"],["CG2","HG22"],["CG2","HG23"],["CD1","HD11"],["CD1","HD12"],["CD1","HD13"]],
                "ASP" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["OD2","HD2"]],
                "LEU" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CG","HG"],
						["CG1","HG1"],["CD1","HD11"],["CD1","HD12"],["CD1","HD13"],
						["CD2","HD21"],["CD2","HD22"],["CD2","HD23"]],
                "ALA" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB1"],["CB","HB2"],["CB","HB3"]],
                "GLN" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CG","HG3"],
						["CG","HG2"],["NE2","HE21"],["NE2","HE22"]],
                "GLY" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA2"],["CA","HA3"]],
                "CYS" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["SG","HG"]],
                "TRP" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CD1","HD1"],
						["NE1","HE1"],["CE3","HE3"],["CZ2","HZ2"],["CH2","HH2"],["CZ3","HZ3"]],
                "TYR" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CD1","HD1"],
						["CD2","HD2"],["CE1","HE1"],["CE2","HE2"],["OH","HH"]],
                "LYS" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["CG","HG3"],
						["CG","HG2"],["CD","HD3"],["CD","HD2"],["CE","HE3"],["CE","HE2"],["NZ","HZ1"],["NZ","HZ2"],["NZ","HZ3"]],
                "PRO" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["N","HN"],["CA","HA"],["CB","HB2"],["CB","HB3"],
						["CG","HG2"],["CG","HG3"],["CD","HD2"],["CD","HD3"]],
                "THR" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB"],["OG1","HG1"],["CG2","HG21"],
						["CG2","HG22"],["CG2","HG23"]],
                "SER" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["OG","HG"]],
                "VAL" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB"],["CG1","HG11"],["CG1","HG12"],
						["CG1","HG13"],["CG2","HG21"],["CG2","HG22"],["CG2","HG23"]],
                "HIS" => [["N","H1"],["N","H2"],["N","H3"],["N","H"],["CA","HA"],["CB","HB3"],["CB","HB2"],["ND1","HD1"],
						["CD2","HD2"],["CE1","HE1"]],
				"NAG" => [],
				"HOH" => [],
				" ZN" => []
)


text:
module TestExamples2DEuler

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler" begin
  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [9.321181253186009e-7, 1.4181210743438511e-6, 1.4181210743487851e-6, 4.824553091276693e-6],
      linf = [9.577246529612893e-6, 1.1707525976012434e-5, 1.1707525976456523e-5, 4.8869615580926506e-5])
  end

  @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fv.jl"),
      l2   = [0.026440292358506527, 0.013245905852168414, 0.013245905852168479, 0.03912520302609374],
      linf = [0.042130817806361964, 0.022685499230187034, 0.022685499230187922, 0.06999771202145322])
  end

  @trixi_testset "elixir_euler_density_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0010600778457964775, 0.00010600778457634275, 0.00021201556915872665, 2.650194614399671e-5],
      linf = [0.006614198043413566, 0.0006614198043973507, 0.001322839608837334, 0.000165354951256802],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.259440511766445e-6, 2.318888155713922e-6, 2.3188881557894307e-6, 6.3327863238858925e-6],
      linf = [1.498738264560373e-5, 1.9182011928187137e-5, 1.918201192685487e-5, 6.0526717141407005e-5])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
      linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
  end

  @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481471610306124, 0.027694280613944234, 0.027697905866996532, 0.12932052501462554],
      linf = [0.31052098400669004, 0.3481295959664616, 0.34807152194137336, 1.1044947556170719],
      maxiters = 10,
      surface_flux = flux_kennedy_gruber,
      volume_flux = flux_kennedy_gruber)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481122603050542, 0.027662840593087695, 0.027665658732350273, 0.12927455860656786],
      linf = [0.3110089578739834, 0.34888111987218107, 0.3488278669826813, 1.1056349046774305],
      maxiters = 10,
      surface_flux = flux_chandrashekar,
      volume_flux = flux_chandrashekar)
  end

  @trixi_testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.05380629130119074, 0.04696798008325309, 0.04697067787841479, 0.19687382235494968],
      linf = [0.18527440131928286, 0.2404798030563736, 0.23269573860381076, 0.6874012187446894])
  end

  @trixi_testset "elixir_euler_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
      l2   = [0.14170569763947993, 0.11647068900798814, 0.11647072556898294, 0.3391989213659599],
      linf = [1.6544204510794196, 1.35194638484646, 1.3519463848472744, 1.831228461662809],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"),
      l2   = [0.4758794741390833, 0.21045415565179362, 0.21045325630191866, 0.7022517958549878],
      linf = [1.710832148442441, 0.9711663578827681, 0.9703787873632452, 2.9619758810532653],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
      l2   = [0.472445774440313, 0.2090782039442978, 0.20885558673697927, 0.700569533591275],
      linf = [1.7066492792835155, 0.9856122336679919, 0.9784316656930644, 2.9372978989672873],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl with mortars" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
      l2   = [0.016486406327766923, 0.03097329879894433, 0.03101012918167401, 0.15157175775429868],
      linf = [0.27688647744873407, 0.5653724536715139, 0.565695523611447, 2.513047611639946],
      refinement_patches=(
          (type="box", coordinates_min=(-0.25, -0.25), coordinates_max=(0.25, 0.25)),
          (type="box", coordinates_min=(-0.125, -0.125), coordinates_max=(0.125, 0.125)),),
      initial_refinement_level = 4,
      maxiters = 5)
  end

  @trixi_testset "elixir_euler_blast_wave_neuralnetwork_cnn.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_cnn.jl"),
      l2   = [0.4795795496408325, 0.2125148972465021, 0.21311260934645868, 0.7033388737692883],
      linf = [1.8295385992182336, 0.9687795218482794, 0.9616033072376108, 2.9513245978047133],
      initial_refinement_level = 4,
      maxiters = 50)
  end

  @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_pure_fv.jl"),
      l2   = [0.39957047631960346, 0.21006912294983154, 0.21006903549932, 0.6280328163981136],
      linf = [2.20417889887697, 1.5487238480003327, 1.5486788679247812, 2.4656795949035857],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
      l2   = [0.6835576416907511, 0.2839963955262972, 0.28399565983676, 0.7229447806293277],
      linf = [3.0969614882801393, 1.7967947300740248, 1.7967508302506658, 3.040149575567518],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.4866953770742574, 0.1673477470091984, 0.16734774700934, 0.6184367248923149],
      linf = [2.6724832723962053, 1.2916089288910635, 1.2916089289001427, 6.474699399394252],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_neuralnetwork_perssonperaire.jl"),
      l2   = [0.0845430093623868, 0.09271459184623232, 0.09271459184623232, 0.4377291875101709],
      linf = [1.3608553480069898, 1.6822884847136004, 1.6822884847135997, 4.220147414536653],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_positivity.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"),
      l2   = [0.48862067511841695, 0.16787541578869494, 0.16787541578869422, 0.6184319933114926],
      linf = [2.6766520821013002, 1.2910938760258996, 1.2910938760258899, 6.473385481404865],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_blob_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_mortar.jl"),
      l2   = [0.22271619518391986, 0.6284824759323494, 0.24864213447943648, 2.9591811489995474],
      linf = [9.15245400430106, 24.96562810334389, 10.388109127032374, 101.20581544156934],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blob_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
      l2   = [0.2086261501910662, 1.2118352377894666, 0.10255333189606497, 5.296238138639236],
      linf = [14.829071984498198, 74.12967742435727, 6.863554388300223, 303.58813147491134],
      tspan = (0.0, 0.12))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      l2   = [0.055691508271624536, 0.032986009333751655, 0.05224390923711999, 0.08009536362771563],
      linf = [0.24043622527087494, 0.1660878796929941, 0.12355946691711608, 0.2694290787257758],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr.jl"),
      l2   = [0.05569452733654995, 0.033107109983417926, 0.05223609622852158, 0.08007777597488817],
      linf = [0.2535807803900303, 0.17397028249895308, 0.12321616095649354, 0.269046666668995],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr_neuralnetwork_perssonperaire.jl"),
      l2   = [0.0009823702998067061, 0.004943231496200673, 0.0048604522073091815, 0.00496983530893294],
      linf = [0.00855717053383187, 0.02087422420794427, 0.017121993783086185, 0.02720703869972585],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_colliding_flow.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"),
      l2   = [0.007237139090503349, 0.044887582765386916, 1.0453570959003603e-6, 0.6627307840935432],
      linf = [0.19437260992446315, 0.5554343646648533, 5.943891455255412e-5, 15.188919846360125],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_colliding_flow_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow_amr.jl"),
      l2   = [0.006768801432802192, 0.032184992228603666, 6.923887797276484e-7, 0.6784222932398366],
      linf = [0.2508663007713608, 0.4097017076529792, 0.0003528986458217968, 22.435474993016918],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_astro_jet_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_astro_jet_amr.jl"),
      l2   = [0.011338365293662804, 10.09743543555765, 0.00392429463200361, 4031.7811487690506],
      linf = [3.3178633141984193, 2993.6445033486402, 8.031723414357423, 1.1918867260293828e6],
      tspan = (0.0, 1.0e-7))
  end

  @trixi_testset "elixir_euler_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [3.53375983916925e-6, 0.0032123259330577325, 0.00321232443824996, 0.004547280616310348],
      linf = [7.719164482999918e-5, 0.030543222729985442, 0.0304822911023237, 0.042888536761282126])
  end

  @trixi_testset "elixir_euler_vortex_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar.jl"),
      l2   = [2.110390460364181e-6, 2.7230027429598542e-5, 3.657273339760332e-5, 8.735519818394382e-5],
      linf = [5.9743882399154735e-5, 0.000731856753784843, 0.0007915976735435315, 0.0022215051634404404])
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1104708260273183e-6, 2.7351867233787183e-5, 3.657480146175434e-5, 8.762994487073319e-5],
      linf = [5.9784376802540784e-5, 0.0007362412469449353, 0.0007903523992238703, 0.0022167215782289418])
  end

  @trixi_testset "elixir_euler_vortex_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_shockcapturing.jl"),
      l2   = [3.7412276700966986e-6, 5.4150680347525463e-5, 5.415287159571093e-5, 0.0001542834620109727],
      linf = [8.473507257800161e-5, 0.0009317864493174621, 0.0009371841830909666, 0.0030735931384739956])
  end

  @trixi_testset "elixir_euler_vortex_mortar_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_shockcapturing.jl"),
      l2   = [2.1106705938658844e-6, 2.73520503434512e-5, 3.6574899482028756e-5, 8.763675702381447e-5],
      linf = [5.9783989705297635e-5, 0.000736247574993687, 0.0007903506762051427, 0.002216704220305843])
  end

  @trixi_testset "elixir_euler_vortex_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_amr.jl"),
      l2   = [2.120552206480055e-6, 0.003281541473561042, 0.003280625257336616, 0.004645872821313438],
      linf = [4.500266027052113e-5, 0.031765399304366726, 0.03179340562764421, 0.04563622772500864])
  end

  @trixi_testset "elixir_euler_ec.jl with boundary_condition_slip_wall" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.06176845600430613, 0.05020624043492084, 0.05021389111189423, 0.22592682624517807],
      linf = [0.29347582879609024, 0.3108124923286465, 0.3107380389949771, 1.054035804988522],
      periodicity = false, boundary_conditions = boundary_condition_slip_wall)
  end
end

# Coverage test for all initial conditions
@testset "Compressible Euler: Tests for initial conditions" begin
  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_density_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.0034896883781420263, 0.0034896883781420245, 0.003489688378142025, 0.0034896883781420345],
      linf = [0.04763147976015558, 0.04763147976015536, 0.04763147976015558, 0.04763147976015514],
      maxiters = 1,
      initial_condition = initial_condition_density_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.00021747926167733352, 0.002201047574822531, 0.0022010475748225305, 0.010855330793503795],
      linf = [0.005450371783679131, 0.031279706118324424, 0.031279706118324535, 0.14847506763666196],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_pressure_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_density_pressure_pulse" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [0.0034736657511653182, 0.005490934309849779, 0.005490934309849781, 0.015625174846479427],
      linf = [0.04658987593455355, 0.07335595216951218, 0.07335595216951196, 0.21083755091992007],
      maxiters = 1,
      initial_condition = Trixi.initial_condition_density_pressure_pulse)
  end

  @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [1.1790213022362371e-16, 8.580657423476384e-17, 1.3082387431804115e-16, 1.6182739965672862e-15],
      linf = [3.3306690738754696e-16, 2.220446049250313e-16, 5.273559366969494e-16, 3.552713678800501e-15],
      maxiters = 1,
      initial_condition = initial_condition_constant)
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl one step" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0021196114178949396, 0.010703549234544042, 0.01070354923454404, 0.10719124037195142],
      linf = [0.11987270645890724, 0.7468615461136827, 0.7468615461136827, 3.910689155287799],
      maxiters=1)
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl one step with initial_condition_medium_sedov_blast_wave" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.00211844592093308, 0.01069763196354363, 0.01069763196354363, 0.1071519081170261],
      linf = [0.11981344588212117, 0.7464822093541509, 0.7464822093541509, 3.9091850198658853],
      maxiters=1, initial_condition=initial_condition_medium_sedov_blast_wave)
  end
end

end # module


text:
# """
#     isWaveguide(domain::Domain)
# """
# function isWaveguide(domain::Domain)
#     return domain.domain_type ∈ [:planar_waveguide, :planar_waveguide_background, :halfspace, :pc_waveguide, :pc_waveguide_background]
# end

#
# """
#     isBulkWaveguide(domain::Domain))
# """
# function isBulkWaveguide(domain::Domain)
#     return domain.domain_type ∈ [:bulk_planar_waveguide_x, :bulk_planar_waveguide_y, :bulk_pc_waveguide_x, :bulk_pc_waveguide_y]
# end

#
# """
#     isBackground(domain::Domain)
# """
# function isBackground(domain)
#     return domain.domain_type ∈ [:background, :planar_waveguide_background, :pc_waveguide_background]
# end

# """
#     isDefect(domain::Domain))
# """
# function isDefect(domain::Domain)
#     return domain.domain_type ∈ [:defect, :site_defect, :line_defect, :pc_waveguide]
# end

# # """
# #     isPC(dom::Domain)
# # """
# # function isPC(dom::Domain)
# #     return dom.domain_type ∈ [:pc, :pc_waveguide, :pc_waveguide_background]
# # end
#
#
# """
#     isPlanar(dom::Domain)
# """
# function isPlanar(dom::Domain)
#     return dom.domain_type ∈ [:planar_waveguide, :planar_waveguide_background, :bulk_planar_waveguide_x, :bulk_planar_waveguide_y]
# end

"""
    isHalfSpace(dom::Domain)
"""
function isHalfSpace(dom::Domain)
    return dom.domain_type ∈ [:halfspace_waveguide]
end

"""
    isDirichlet(bc)
"""
isDirichlet(bc::DirichletBC) = true
isDirichlet(bc) = false

"""
    isNeumann(bc)
"""
isNeumann(bc::NeumannBC) = true
isNeumann(bc) = false

"""
    isOpen(bc)
"""
isMatched(bc::MatchedBC) = true
isMatched(bc) = false

"""
    isPeriodic(bc)
"""
isFloquet(bc::FloquetBC) = true
isFloquet(bc) = false

"""
    isPMLout(bl)
"""
isPML(bl::PML) = true
isPML(bl) = false

"""
    isPMLin(bl)
"""
iscPML(bl::cPML) = true
iscPML(bl) = false

"""
    isNone(bl)
"""

isnoBL(bl::noBL) = true
isnoBL(bl) = false


text:
ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  inits1 = [
    Dict("alpha" => 0,"beta" => 0,"tau" => 1),
    Dict("alpha" => 1,"beta" => 2,"tau" => 1),
    Dict("alpha" => 3,"beta" => 3,"tau" => 2),
    Dict("alpha" => 5,"beta" => 2,"tau" => 5)
  ]

  function test(init::Array{Dict{String, Int64},1})
    res = map((x)->convert(Dict{String, Any}, x), init)
  end

  function test(init::Array{Dict{String, Float64},1})
    res = map((x)->convert(Dict{String, Any}, x), init)
  end

  function test(init::Array{Dict{String, Number},1})
    res = map((x)->convert(Dict{String, Any}, x), init)
  end

  function test(init::Array{Dict{String, Any},1})
    res = map((x)->convert(Dict{String, Any}, x), init)
  end

  line = "
  model {
    for (i in 1:n) {
          mu[i] <- alpha + beta*(x[i] - x.bar);
          y[i]   ~ dnorm(mu[i],tau);
    }
    x.bar   <- mean(x[]);
    alpha    ~ dnorm(0.0,1.0E-4);
    beta     ~ dnorm(0.0,1.0E-4);
    tau      ~ dgamma(1.0E-3,1.0E-3);
    sigma   <- 1.0/sqrt(tau);
  }
  "

  data = Dict{String, Any}()
  data["x"] = [1, 2, 3, 4, 5]
  data["y"] = [1, 3, 3, 3, 5]
  data["n"] = 5

  inits = test(inits1)

  monitors = Dict{String, Bool}(
    "alpha" => true,
    "beta" => true,
    "tau" => true,
    "sigma" => true,
  )

  jagsmodel = Jagsmodel(name="line1", model=line, monitor=monitors,
    ncommands=1, nchains=4, adapt=1000, nsamples=10000, thin=1,
    deviance=true, dic=true, popt=true, pdir=ProjDir);

  println("\nJagsmodel that will be used:")
  jagsmodel |> display
  println("Input observed data dictionary:")
  data |> display
  println("\nInput initial values dictionary:")
  inits |> display

  cd(ProjDir)
  isdir("tmp") &&
    rm("tmp", recursive=true);

end


text:
# global state management
#
# The functionality in this file serves to create a coherent environment to perform
# computations in, with support for Julia constructs like tasks (and executing those on
# multiple threads), using a GPU of your choice (with the ability to reset that device, or
# use different devices on different threads).

export context, context!, device, device!, device_reset!, deviceid, stream, stream!, pool

## hooks

"""
    CUDA.attaskswitch(f::Function)

Register a function to be called after switching to or initializing a task on a thread.

Use this hook to invalidate thread-local state that depends on the current task.
"""
attaskswitch(f::Function) = (pushfirst!(task_hooks, f); nothing)
const task_hooks = []
_attaskswitch() = foreach(f->Base.invokelatest(f), task_hooks)

"""
    CUDA.atdeviceswitch(f::Function)

Register a function to be called after switching to or initializing a device on a thread.

Use this hook to invalidate thread-local state that depends on the current device. If that
state is also context dependent, be sure to query the context in your callback.
"""
atdeviceswitch(f::Function) = (pushfirst!(device_switch_hooks, f); nothing)
const device_switch_hooks = []
_atdeviceswitch() = foreach(f->Base.invokelatest(f), device_switch_hooks)

"""
    CUDA.atdevicereset(f::Function)

Register a function to be called after resetting devices. The function is passed one
argument: the device which has been reset.

Use this hook to invalidate global state that depends on the current device.
"""
atdevicereset(f::Function) = (pushfirst!(device_reset_hooks, f); nothing)
const device_reset_hooks = []
_atdevicereset(dev) = foreach(f->Base.invokelatest(f, dev), device_reset_hooks)

## initialization

# Many CUDA APIs, like the CUDA driver API used by CUDA.jl, use implicit thread-local state
# to determine, e.g., which device to use. With Julia however, code is grouped in tasks.
# Execution can switch between them, and tasks can be executing on (and in the future
# migrate between) different threads. To synchronize these two worlds, we try to detect task
# and thread switches, making sure CUDA's state is mirrored appropriately.

# If you need to maintain your own task-local state, subscribe to device and task switch
# events using [`CUDA.atdeviceswitch`](@ref) and [`CUDA.attaskswitch`](@ref) for proper
# invalidation. If your state is device-specific, but global (i.e. not task-bound), it
# suffices to index your state with the current [`deviceid()`](@ref) and invalidate that
# state when the device is reset by subscribing to [`CUDA.atdevicereset()`](@ref).

@inline function detect_state_changes()
    detect_task_switches()

    # make sure to initialize the CUDA context, which is responsible for triggering a
    # device switch event after the device has been reset.
    initialize_cuda_context()

    return
end

@inline function detect_task_switches()
    tid = Threads.threadid()

    # detect when a different task is now executing on a thread
    task = @inbounds thread_tasks[tid]
    if task === nothing || task.value === nothing || task.value::Task !== current_task()
        switched_tasks(tid, current_task())
    end

    return
end

const thread_tasks = Union{Nothing,WeakRef}[]
@noinline function switched_tasks(tid::Int, task::Task)
    thread_tasks[tid] = WeakRef(task)
    _attaskswitch()

    # switch contexts if task switched to was already bound to one
    ctx = get(task_local_storage(), :CuContext, nothing)
    if ctx !== nothing && isvalid(ctx)
        # NOTE: the context may be invalid if another task reset it (which we detect here
        #       since we can't touch other tasks' local state from `device_reset!`)
        context!(ctx)
    else
        thread_state[tid] = nothing  # trigger `initialize_thread` down the line
        # NOTE: actually deactivating the CUDA context would be more correct,
        #       but that confuses CUDA and leads to invalid contexts later on.
    end
end

# prepare for a CUDA API call, which requires a CUDA device to be selected.
# we can't do this eagerly because it consumes quite some memory.
@inline function initialize_cuda_context()
    tid = Threads.threadid()

    # initialize CUDA state when first executing on a thread
    state = @inbounds thread_state[tid]
    if state === nothing
        initialize_current_thread()
    end

    return
end

# the default device unitialized tasks will use, set when switching devices.
# this behavior differs from the CUDA Runtime, where device 0 is always used.
# this setting won't be used when switching tasks on a pre-initialized thread.
const default_device = Ref{Union{Nothing,CuDevice}}(nothing)

# CUDA uses thread-bound state, but calling CuCurrent* all the time is expensive,
# so we maintain our own thread-local copy keeping track of the current CUDA state.
const CuCurrentState = NamedTuple{(:ctx, :dev), Tuple{CuContext,CuDevice}}
const thread_state = Union{Nothing,CuCurrentState}[]
@noinline function initialize_current_thread()
    dev = something(default_device[], CuDevice(0))
    device!(dev)

    # NOTE: we can't be compatible with externally initialize contexts here (i.e., reuse
    #       the CuCurrentContext we don't know about) because of how device_reset! works:
    #       contexts that got reset remain bound, and are not discernible from regular ones.
end

## context-based API

"""
    context()::CuContext

Get or create a CUDA context for the current thread (as opposed to
`CuCurrentContext` which may return `nothing` if there is no context bound to the
current thread).
"""
@inline function context()
    detect_state_changes()
    tid = Threads.threadid()

    state = @inbounds thread_state[tid]::CuCurrentState

    if Base.JLOptions().debug_level >= 2
        @assert state.ctx == CuCurrentContext()
    end
    state.ctx
end

"""
    context!(ctx::CuContext)

Bind the current host thread to the context `ctx`.

Note that the contexts used with this call should be previously acquired by calling
[`context`](@ref), and not arbitrary contexts created by calling the `CuContext` constructor.
"""
function context!(ctx::CuContext)
    detect_task_switches()  # we got a CuDevice, so CUDA is initialized already
    tid = Threads.threadid()

    # update the thread-local state
    state = @inbounds thread_state[tid]
    if state === nothing || state.ctx != ctx
        activate(ctx)
        dev = CuCurrentDevice()
        thread_state[tid] = (;ctx=ctx, dev=dev)
        _atdeviceswitch()
    end

    # update the task-local state
    task_local_storage(:CuContext, ctx)

    return
end

"""
    context!(f, ctx)

Sets the active context for the duration of `f`.
"""
function context!(f::Function, ctx::CuContext)
    @assert isvalid(ctx)
    old_ctx = CuCurrentContext()
    if ctx != old_ctx
        context!(ctx)
    end
    try
        f()
    finally
        if ctx !== old_ctx && old_ctx !== nothing
            context!(old_ctx)
        end
    end
end

# macro version for maximal performance (avoiding closures)
macro context!(ctx_expr, expr)
    quote
        ctx = $(esc(ctx_expr))
        @assert isvalid(ctx)
        old_ctx = CuCurrentContext()
        if ctx != old_ctx
            context!(ctx)
        end
        try
            $(esc(expr))
        finally
            if ctx !== old_ctx && old_ctx !== nothing
                context!(old_ctx)
            end
        end
    end
end

## device-based API

"""
    device()::CuDevice

Get the CUDA device for the current thread, similar to how [`context()`](@ref) works
compared to [`CuCurrentContext()`](@ref).
"""
@inline function device()
    detect_state_changes()
    tid = Threads.threadid()

    state = @inbounds thread_state[tid]::CuCurrentState

    if Base.JLOptions().debug_level >= 2
        @assert state.dev == CuCurrentDevice()
    end
    state.dev
end

const __device_contexts = Union{Nothing,CuContext}[]
device_context(i) = @after_init(@inbounds __device_contexts[i])
device_context!(i, ctx) = @after_init(@inbounds __device_contexts[i] = ctx)

function context(dev::CuDevice)
    devidx = deviceid(dev)+1

    # querying the primary context for a device is expensive, so cache it
    ctx = device_context(devidx)
    if ctx !== nothing
        return ctx
    end

    if capability(dev) < v"5"
        @warn("""Your $(name(dev)) GPU does not meet the minimal required compute capability ($(capability(dev)) < 5.0).
                 Some functionality might not work. For a fully-supported set-up, please use an older version of CUDA.jl""",
              maxlog=1, _id=devidx)
    end

    # configure the primary context
    pctx = CuPrimaryContext(dev)
    ctx = CuContext(pctx)
    device_context!(devidx, ctx)
    return ctx
end

"""
    device!(dev::Integer)
    device!(dev::CuDevice)

Sets `dev` as the current active device for the calling host thread. Devices can be
specified by integer id, or as a `CuDevice` (slightly faster).

If your library or code needs to perform an action when the active device changes,
add a hook using [`CUDA.atdeviceswitch`](@ref).
"""
function device!(dev::CuDevice, flags=nothing)
    detect_task_switches()  # we got a CuDevice, so CUDA is initialized already
    tid = Threads.threadid()

    # configure the primary context flags
    if flags !== nothing
        devidx = deviceid(dev)+1
        if device_context(devidx) !== nothing
            error("Cannot set flags for an active device. Do so before calling any CUDA function, or reset the device first.")
        end
        cuDevicePrimaryCtxSetFlags(dev, flags)
    end

    # make this device the new default
    default_device[] = dev

    # bail out if switching to the current device
    state = @inbounds thread_state[tid]
    if state !== nothing && state.dev == dev
        return
    end

    # actually switch contexts
    ctx = context(dev)
    context!(ctx)
end

"""
    device!(f, dev)

Sets the active device for the duration of `f`.

Note that this call is intended for temporarily switching devices, and does not change the
default device used to initialize new threads or tasks.
"""
function device!(f::Function, dev::CuDevice)
    ctx = context(dev)
    context!(f, ctx)
end

macro device!(dev_expr, expr)
    quote
        dev = $(esc(dev_expr))
        ctx = context(dev)
        @context! ctx $(esc(expr))
    end
end

"""
    device_reset!(dev::CuDevice=device())

Reset the CUDA state associated with a device. This call with release the underlying
context, at which point any objects allocated in that context will be invalidated.

If your library or code needs to perform an action when the active context changes,
add a hook using [`CUDA.atdevicereset`](@ref). Resetting the device will also cause
the [`CUDA.atdeviceswitch`](@ref) hook to fire when `initialize_cuda_context` is called,
so it is generally not needed to subscribe to the reset hook specifically.

!!! warning

    This function is broken on CUDA 11.2 when using the CUDA memory pool (the default).
    If you need to reset the device, use another memory pool by setting the
    `JULIA_CUDA_MEMORY_POOL` environment variable to, e.g., `binned` before importing
    this package.
"""
function device_reset!(dev::CuDevice=device())
    if async_alloc[] # NVIDIA bug #3240770
        @error """Due to a bug in CUDA, resetting the device is not possible on CUDA 11.2 when using the stream-ordered memory allocator.

                  If you are calling this function to free memory, that may not be required anymore
                  as the stream-ordered memory allocator releases memory much more eagerly.
                  If you do need this functionality, switch to another memory pool by setting
                  the `JULIA_CUDA_MEMORY_POOL` environment variable to, e.g., `binned`."""
        return
    end

    # unconditionally reset the primary context (don't just release it),
    # as there might be users outside of CUDA.jl
    pctx = CuPrimaryContext(dev)
    unsafe_reset!(pctx)

    # wipe the thread-local state for all threads using this device
    for (tid, state) in enumerate(thread_state)
        if state !== nothing && state.dev == dev
            thread_state[tid] = nothing
        end
    end

    # wipe the device-specific state
    devidx = deviceid(dev)+1
    device_context!(devidx, nothing)

    _atdevicereset(dev)

    return
end

## integer device-based API

device!(dev::Integer, flags=nothing) = device!(CuDevice(dev), flags)
device!(f::Function, dev::Integer) = device!(f, CuDevice(dev))

"""
    deviceid()::Int
    deviceid(dev::CuDevice)::Int

Get the ID number of the current device of execution. This is a 0-indexed number,
corresponding to the device ID as known to CUDA.
"""
deviceid(dev::CuDevice=device()) = Int(convert(CUdevice, dev))

## helpers

# helper struct to maintain state per device
# - make it possible to index directly with CuDevice (without converting to integer index)
# - initialize function to fill state based on a constructor function
# - automatically wiping state on device reset
struct PerDevice{T,F} <: AbstractVector{T}
    inner::Vector{T}
    ctor::F

    PerDevice{T,F}(ctor::F) where {T,F<:Function} = new(Vector{T}(), ctor)
end

PerDevice{T}(ctor::F) where {T,F} = PerDevice{T,F}(ctor)

function initialize!(x::PerDevice, n::Integer)
    @assert isempty(x.inner)

    resize!(x.inner, n)
    for i in 0:n-1
        x[i] = x.ctor(CuDevice(i))
    end

    atdevicereset() do dev
        x[dev] = x.ctor(dev)
    end

    return
end

# 0-based indexing for using CUDA device identifiers
Base.getindex(x::PerDevice, devidx::Integer) = x.inner[devidx+1]
Base.setindex!(x::PerDevice, val, devidx::Integer) = (x.inner[devidx+1] = val; )

# indexing using CuDevice objecs
Base.getindex(x::PerDevice, dev::CuDevice) = x[deviceid(dev)]
Base.setindex!(x::PerDevice, val, dev::CuDevice) = (x[deviceid(dev)] = val; )

Base.length(x::PerDevice) = length(x.inner)
Base.size(x::PerDevice) = size(x.inner)

function Base.show(io::IO, mime::MIME"text/plain", x::PerDevice{T}) where {T}
    print(io, "PerDevice{$T} with $(length(x)) entries")
end

## math mode

@enum MathMode begin
    # use prescribed precision and standardized arithmetic for all calculations.
    # this may serialize operations, and reduce performance.
    PEDANTIC_MATH

    # use at least the required precision, and allow reordering operations for performance.
    DEFAULT_MATH

    # additionally allow downcasting operations for better use of hardware resources.
    # whenever possible the `precision` flag passed to `math_mode!` will be used
    # to constrain those downcasts.
    FAST_MATH
end

# math mode and precision are sticky (once set on a task, inherit to newly created tasks)
const default_math_mode = Ref{Union{Nothing,MathMode}}(nothing)
const default_math_precision = Ref{Union{Nothing,Symbol}}(nothing)

function math_mode!(mode::MathMode; precision=nothing)
    # make sure we initialize first, or recursion might overwrite the math mode
    ctx = context()

    tls = task_local_storage()
    tls[(:CUDA, :math_mode)] = mode
    default_math_mode[] = mode
    if precision !== nothing
        tls[(:CUDA, :math_precision)] = precision
        default_math_precision[] = precision
    end

    # reapply the CUBLAS math mode if it had been set already
    cublas_handle = get(tls, (:CUBLAS, ctx), nothing)
    if cublas_handle !== nothing
        CUBLAS.math_mode!(cublas_handle, mode)
    end

    return
end

math_mode() =
    get!(task_local_storage(), (:CUDA, :math_mode)) do
        something(default_math_mode[],
                  Base.JLOptions().fast_math==1 ? FAST_MATH : DEFAULT_MATH)
    end
math_precision() =
    get!(task_local_storage(), (:CUDA, :math_precision)) do
        something(default_math_precision[], :TensorFloat32)
    end

## streams

# thread cache for task-local streams
const thread_streams = Vector{Union{Nothing,CuStream}}()

"""
    stream()

Get the CUDA stream that should be used as the default one for the currently executing task.
"""
@inline function stream()
    detect_state_changes()
    tid = Threads.threadid()

    if @inbounds thread_streams[tid] === nothing
        ctx = context()
        thread_streams[tid] = get!(task_local_storage(), (:CuStream, ctx)) do
            stream = CuStream(; flags=STREAM_NON_BLOCKING)

            t = current_task()
            tptr = pointer_from_objref(current_task())
            tptrstr = string(convert(UInt, tptr), base=16, pad=Sys.WORD_SIZE>>2)
            NVTX.nvtxNameCuStreamA(stream, "Task(0x$tptrstr)")

            stream
        end
    end
    something(@inbounds thread_streams[tid])
end

function set_library_streams(s)
    CUBLAS.set_stream(s)
    CUSPARSE.set_stream(s)
    CUSOLVER.set_stream(s)
    CURAND.set_stream(s)
    CUFFT.set_stream(s)

    CUDNN.set_stream(s)
    CUTENSOR.set_stream(s)
end

function stream!(s::CuStream)
    # task switch detected by context()
    tid = Threads.threadid()

    ctx = context()
    task_local_storage((:CuStream, ctx), s)

    # update the thread cache
    @inbounds thread_streams[tid] = s

    set_library_streams(s)
end

function stream!(f::Function, s::CuStream)
    # task switch detected by stream()
    tid = Threads.threadid()

    # NOTE: we can't read `thread_streams` directly here, or could end up with `nothing`,
    #       and we need a valid stream to fall back to and reset the library handles with.
    old_s = stream()
    try
        return task_local_storage(:CuStream, s) do
            thread_streams[tid] = s
            set_library_streams(s)
            f()
        end
    finally
        thread_streams[tid] = old_s
        set_library_streams(old_s)
    end
end

## memory pools

const memory_pools = Dict{CuContext,CuMemoryPool}()

function pool()
    if CUDA.version() < v"11.2"
        return nothing
    end

    return get!(memory_pools, context()) do
        CuMemoryPool(device())
    end
end


text:
#ConventionalApp: Launcher file generation
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
#const DEFAULT_JULIACMD="julia" #Do not use: Will not exist on most systems
const BASHSCRIPT_PATH="\$SCRIPTDIR"
const RUNSCRIPT_JULIA="run.jl"
#const RUNSCRIPT_BASH="run"

#==String generators
===============================================================================#
function gen_juliacmd()
	#Base.julia_cmd() #TODO: Investigate use of julia_cmd() instead.
	juliacmd = joinpath(Sys.BINDIR, Base.julia_exename())
	return juliacmd
end

function gen_runfile(::Type{String}, proj::Project)
	app_path = proj.root
	app_name = basename(app_path)
	runfile="""
using ConventionalApp

proj = Project(@__DIR__)
setup_env(proj)
@include_startup(proj)

#show(activeproject)

Main.include("src/$app_name.jl")
Main.$app_name.run_app()
"""
end

function gen_bashfile(::Type{String}, juliacmd::String, projpath::String)
	bashfile="""
#!/bin/bash
SCRIPTPATH=`readlink -f \$0` #This file here
SCRIPTDIR=`dirname \$SCRIPTPATH`

PROJPATH="$projpath"
JULIASCRIPT="\$PROJPATH/$RUNSCRIPT_JULIA"
JULIACMD="$juliacmd"
exec \$JULIACMD -i --color=yes --startup-file=no --project="\$PROJPATH" "\$JULIASCRIPT" "\$@"
"""
	return bashfile
end

function gen_winshortcut(::Type{String}, juliacmd::String, projpath::String)
	info="""
Target:
   \"$juliacmd\" --project=\"$projpath\" -L \"$RUNSCRIPT_JULIA\"

Start in:
    \"$projpath\"
"""
	return info
end

#==File generators (implementation)
===============================================================================#
function _gen_runfile(proj::Project)
	filepath = joinpath(proj.root, RUNSCRIPT_JULIA)
	data = gen_runfile(String, proj)
	open(filepath, "w") do io
		print(io, data)
	end
	return
end

function _gen_bashfile(proj::Project, bashfile::String)
	inprojroot = false
	projpath = inprojroot ? BASHSCRIPT_PATH : proj.root

	data = gen_bashfile(String, gen_juliacmd(), projpath)
	open(bashfile, "w") do io
		print(io, data)
	end
	return
end
_gen_bashfile(proj::Project, bashfile::Nothing) = _gen_bashfile(proj, basename(proj))

function _gen_winshortcut(proj::Project, bashfile::String)
	data = gen_winshortcut(String, gen_juliacmd(), proj.root)
	println()
	@info("To create a Windows shortcut, use the following values:")
	println(data)
	return
end
_gen_winshortcut(proj::Project, bashfile::Nothing) = _gen_winshortcut(proj, basename(proj) * ".lnk")

#==Top-level interface
===============================================================================#
gen_runfile(proj::Project) = _gen_runfile(proj)
gen_runfile(m::Module) = _gen_runfile(Project(m))

gen_bashfile(proj::Project, bashfile=nothing) = _gen_bashfile(proj, bashfile)
gen_bashfile(m::Module, bashfile=nothing) = _gen_bashfile(Project(m), bashfile)

gen_winshortcut(proj::Project, bashfile=nothing) = _gen_winshortcut(proj, bashfile)
gen_winshortcut(m::Module, bashfile=nothing) = _gen_winshortcut(Project(m), bashfile)

#Last line

max_stars_repo_name: ma-laforge/JuliaApplications

max_issues_repo_name: ma-laforge/JuliaApplications

max_forks_repo_name: ma-laforge/JuliaApplications


text:
# Error
# =====

"""
Decoding exception.
"""
struct DecodeError <: Exception
    msg::String
end

function Base.showerror(io::IO, error::DecodeError)
    print(io, "DecodeError: $(error.msg)")
end


text:
function my_f(x,y)
  2x + 3y
end

function derivative_of_my_f(x,y)
  ForwardDiff.derivative(x -> my_f(x,y), x)
end


text:
using Jive

ignores = Set()

for db in ["hive", "jdbc", "odbc"]
    push!(ignores, joinpath("adapters", db))
end

# LibPQ v0.11.1
#     ERROR: LoadError: ArgumentError: cannot convert NULL to string
push!(ignores, "adapters/postgresql/copy_test")

# skip options
push!(ignores, "adapters/postgresql/options.jl")
push!(ignores, "adapters/mysql/options.jl")

runtests(@__DIR__, skip=collect(ignores))


text:
#===== Projectile model =====
We test the line plot!(x->x, x->-24*x+375, 0., 50., line=2., color="red", linestyle=:solid, legend=:none)
=#
A = [0. 0.5 0. 0. ; 0. 0. 0. 0. ; 0. 0. 0. 0.7 ; 0. 0. 0. 0.]
X0 = Singleton([0.,5.,100.,0])
U = Singleton([0.,0.,0.,-9.81])
S = ContinuousSystem(A, X0, U)
time_horizon = 20.0
precision = 1e-4
initial_δ = 0.5
algorithm(N, δ) = solve(S, :mode => "check", :partition=>[1:2, 3:4],
                           :vars => [1, 3], :plot_vars => [1, 3], :δ => δ,
                           :T => time_horizon,
                           :property=>LinearConstraintProperty([24., 0., 1, 0],  375.)).satisfied

Properties.tune_δ(algorithm, time_horizon, precision, initial_δ)


text:
"""
Created on Wed 27 May 2020

@author: Yoann Pradat

Plot metrics as produced by the `rank_by_consensus` function.
"""

using  DataFrames
import PyPlot
const  plt = PyPlot

function plot_metric!(ax::PyCall.PyObject, df_metrics::DataFrame; axvline::Bool, variable::Symbol, color::String, label::String, legend::String, title::String)
    #### plot
    ax.plot(
        df_metrics.rank,
        df_metrics[:,variable],
        color  = color,
        marker = "o",
        label  = legend
    )

    ax.axvline(
        df_metrics[argmax(df_metrics[:,variable]),:rank],
        ymin  = 0,
        ymax  = 1,
        color = color,
        ls    = "--",
        lw    = 2
    )

    ax.text(
        x          = -0.1,
        y          = 1.1,
        s          = label,
        transform  = ax.transAxes,
        fontsize   = 15,
        fontweight = "bold"
    )

    #### axes
    ax.set_xlabel(
        xlabel     = "rank",
        fontsize   = 15,
        fontweight = "bold"
    )
    ax.set_ylabel(
        ylabel     = "$(variable)",
        fontsize   = 15,
        fontweight = "bold"
    )
    ax.spines["top"].set_visible(false)
    ax.spines["right"].set_visible(false)

    #### title
    ax.set_title(
        label=title,
        fontsize=18,
        fontweight="bold"
    )
end

"""
    plot_metrics_consensus(df_metrics, filepath)

Plot metrics as produced by the `rank_by_consensus` function.
"""
function plot_metrics_consensus(df_metrics::DataFrame, filepath::String)
    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(14,10))

    plot_metric!(
        ax[1,1],
        df_metrics,
        axvline     = true,
        variable    = :cophenetic,
        color       = "purple",
        label       = "A",
        legend      = "",
        title       = "Cophenetic correlation"
    )

    plot_metric!(
        ax[1,2],
        df_metrics,
        axvline     = true,
        variable    = :dispersion,
        color       = "purple",
        label       = "B",
        legend      = "",
        title       = "Dispersion"
    )

    plot_metric!(
        ax[2,1],
        df_metrics,
        axvline     = true,
        variable    = :silhouette,
        color       = "purple",
        label       = "C",
        legend      = "",
        title       = "Silhouette"
    )

    plot_metric!(
        ax[2,2],
        df_metrics,
        axvline     = true,
        variable    = :sparse_W,
        color       = "green",
        label       = "D",
        legend      = "basis",
        title       = "Sparseness"
    )

    plot_metric!(
        ax[2,2],
        df_metrics,
        axvline     = true,
        variable    = :sparse_H,
        color       = "red",
        label       = "D",
        legend      = "coefs",
        title       = "Sparseness"
    )

    ax[2,2].set_ylabel("sparsity index")
    ax[2,2].legend(loc="best", fontsize=15, frameon=false)

    plt.subplots_adjust(hspace=0.4, wspace=0.2)
    plt.savefig(filepath, bbox_inches="tight")
end

"""
    plot_metrics_stability(df_metrics, filepath)

Plot metrics as produced by the `rank_by_stability` function.
"""
function plot_metrics_stability(df_metrics::DataFrame, filepath::String)
    fig, ax1 = plt.subplots()

    variable1 = :stability
    color     = "red"

    p1 = ax1.plot(
        df_metrics.rank,
        df_metrics[:,variable1],
        color      = color,
        marker     = "o",
        markersize = 10,
        linestyle  = ":",
    )

    #### axes
    ax1.set_xlabel(
        xlabel     = "rank",
        fontsize   = 12,
        fontweight = "bold"
    )
    ax1.set_ylabel(
        ylabel     = "",
        fontsize   = 15,
        fontweight = "bold"
    )

    ax1.tick_params(axis="y", labelcolor=color)
    ax1.set_xticks(df_metrics.rank)

    #### new ax with twin x axis
    ax2 = ax1.twinx()

    variable2 = :fidelity
    color     = "dodgerblue"

    p2 = ax2.plot(
        df_metrics.rank,
        df_metrics[:,variable2],
        color      = color,
        marker     = "o",
        markersize = 10,
        linestyle  = ":",
    )

    #### axes
    ax2.set_xlabel(
        xlabel     = "rank",
        fontsize   = 12,
        fontweight = "bold"
    )
    ax2.set_ylabel(
        ylabel     = "",
        fontsize   = 15,
        fontweight = "bold"
    )

    ax2.tick_params(axis="y", labelcolor=color)

    #### legend + title
    plt.legend([p1[1], p2[1]], ["$(variable1)", "$(variable2)"], loc="best", fontsize=12)

    plt.title(
        label="Rank by stability",
        fontsize=15,
        fontweight="bold"
    )

    fig.tight_layout()  # otherwise the right y-label is slightly clipped
    plt.savefig(filepath, bbox_inches="tight")
end


text:
#=――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
    NOTE
    The algorithm and rules used for unpluralizing English words was originated
    by the Inflector .NET project, and is copyright Andrew Peters, Scott
    Kirkland and other contributors. It is made available under the MIT
    license. Modifications may have been made to the original algorithm.

    The Inflector .NET is hosted on GitHub at the following URL
    <https://github.com/srkirkland/Inflector>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――=#

export singularize

const SINGULARIZE_RULES = [
    r"moves$"               => s"move",
    r"monies$"              => s"money",
    r"matrices$"            => s"matrix",
    r"([io])rides$"         => s"\1ris",
    r"podes$"               => s"pus",
    r"(qu|wh)izzes$"        => s"\1iz",
    r"zes"                  => s"z",
    r"(vert|ind)ices$"      => s"\1ex",
    r"(cris|ax|test)es$"    => s"\1is",
    r"shoes$"               => s"shoe",
    r"oes$"                 => s"o",
    r"([ml])ice$"           => s"\1ouse",
    r"([iay])nges$"         => s"\1nx",
    r"(x|ch|ss|sh)es$"      => s"\1",
    r"movies$"              => s"movie",
    r"([^aeiouy]|qu)ies$"   => s"\1y",
    r"([lr])ves$"           => s"\1f",
    r"([ht])ives$"          => s"\1ive",
    r"([^f])ves$"           => s"\1fe",
    r"[us]ses$"             => s"us",
    r"ises$"                => s"is",
    r"ses$"                 => s"sis",
    r"s$"                   => s"",
    r"mina$"                => s"men",
    r"([aegmosu])ata$"      => s"\1a",
    r"([ti])a$"             => s"\1um",
    r"a$"                   => s"on",
    r"geese$"               => s"goose",
    r"ae$"                  => s"a",
    r"i$"                   => s"us",
    r"im$"                  => s"",
    r"men$"                 => s"man",
    r"feet$"                => s"foot",
    r"([ae])ux$"            => s"\1u"
]

const IGNORE_SUFFIXES = [" General", "-in-law"]

"""
    singularize(word)

Unpluralize a plural noun (given in canonical capitalization) using heuristics
and lists of exceptions. If the given `word` is not a plural noun, the result
may be unpredictable.

```jldoctest
julia> using EnglishText

julia> singularize("foxes")
"fox"

julia> singularize("data")
"datum"
```
"""
function singularize(s::String)
    # handle cases like Attorneys General, children-in-law
    for suffix in IGNORE_SUFFIXES
        if endswith(s, suffix)
            return singularize(stem(s, length(suffix))) * suffix
        end
    end

    # irregular words
    for (singular, plural) in IRREGULAR_CLS
        if plural == s
            return singular
        end
    end
    # words that do not inflect
    if isnoninflecting(s, true) || isnoninflecting(s, false)
        return s
    end
    # singular words ending in s
    if endswith(s, "ses")
        chopped = stem(s, 2)
        if chopped ∈ SINGLE_S
            return chopped
        end
    end

    # A18: -i to -o (cls)
    for w in A18
        if endswith(s, stem(w, 1, "i"))
            return stem(s, 1, "o")
        end
    end

    for (regex, subst) ∈ SINGULARIZE_RULES
        if occursin(regex, s)
            return replace(s, regex => subst)
        end
    end
    s
end
singularize(s::AbstractString) = singularize(String(s))


text:
using CImGui
using CImGui.CSyntax.CStatic

"""
    ShowExampleAppLayout(p_open::Ref{Bool})
Create a window with multiple child windows.
"""
function ShowExampleAppLayout(p_open::Ref{Bool})
    CImGui.SetNextWindowSize((500, 440), CImGui.ImGuiCond_FirstUseEver)
    if CImGui.Begin("Example: Layout", p_open, CImGui.ImGuiWindowFlags_MenuBar)
        if CImGui.BeginMenuBar()
            if CImGui.BeginMenu("File")
                CImGui.MenuItem("Close") && (p_open[] = false;)
                CImGui.EndMenu()
            end
            CImGui.EndMenuBar()
        end

        selected = @cstatic selected=1 begin
            # left
            CImGui.BeginChild("left pane", (150, 0), true)
            for i = 1:100
                CImGui.Selectable("MyObject $i", selected == i) && (selected = i;)
            end
            CImGui.EndChild()
            CImGui.SameLine()
        end

        # right
        CImGui.BeginGroup()
            # leave room for 1 line below us
            CImGui.BeginChild("item view", (0.0f0, -CImGui.GetFrameHeightWithSpacing()))
                CImGui.Text("MyObject: $selected")
                CImGui.Separator()
                if CImGui.BeginTabBar("##Tabs", CImGui.ImGuiTabBarFlags_None)
                    if CImGui.BeginTabItem("Description")
                        CImGui.TextWrapped("Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. ")
                        CImGui.EndTabItem()
                    end
                    if CImGui.BeginTabItem("Details")
                        CImGui.Text("ID: 0123456789")
                        CImGui.EndTabItem()
                    end
                    CImGui.EndTabBar()
                end
            CImGui.EndChild()
            if CImGui.Button("Revert")
                @info "Trigger Revert | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
            CImGui.SameLine()
            if CImGui.Button("Save")
                @info "Trigger Save | find me here: $(@__FILE__) at line $(@__LINE__)"
            end
        CImGui.EndGroup()
    end
    CImGui.End()
end


text:
using GnuTLS

if isempty(profile["key"])
    hmac(s1,s2,s3,s4) = ""
else
    signature_scheme = get(profile, "signature_scheme", "hmac-sha256")
    isempty(signature_scheme) && (signature_scheme = "hmac-sha256")
    signature_scheme = split(signature_scheme, "-")
    if signature_scheme[1] != "hmac" || length(signature_scheme) != 2
        error("unrecognized signature_scheme")
    end
    const hmacstate = HMACState(eval(symbol(uppercase(signature_scheme[2]))),
                                profile["key"])
    function hmac(s1,s2,s3,s4)
        update(hmacstate, s1)
        update(hmacstate, s2)
        update(hmacstate, s3)
        update(hmacstate, s4)
        hexdigest!(hmacstate)
    end

    if hmac("a","b","c","d") != hmac("a","b","c","d")
        error("GnuTLS HMAC is not working, probably because you have an older version of GnuTLS linked with Libgcrypt rather than Nettle.") # see GnuTLS.jl issue #5
    end
end


text:
using SwagUI
using Test

include("app.jl")


text:
#=
Companion file to the McEliece implementation to test the main functions of it.
It's not really a benchmark. Know that if you run the whole set of parameters,
it's going to take a long time, it isn't the most optimized implementation.

Author: Ivan Alejandro Moreno Soto
Last updated: 2019/July/15
=#

using DelimitedFiles

include("./McEliece.jl")
using .McEliece

"""
    testKeyGeneration(n, t, m, teststrings, publicDirectory, privateDirectory, saveKeysToFile=false)

Times key generation, encryption and decryption of the implementation. If you want, you can save the
keys to text files during the process.

n, t, and m are the parameters for the goppa code. teststrings is the list of strings to encrypt and
decrypt.
"""
function testKeyGeneration(n, t, m, teststrings, publicDirectory, privateDirectory, saveKeysToFile=false)
    publicKey, privateKey = generateKeys(n, t, m)

    if saveKeysToFile
        saveKeys(publicKey, privateKey, publicDirectory, privateDirectory)
    end

    for plaintext in teststrings
        println("Timing encryption of $(plaintext)")
        @time cyphertext = encrypt(plaintext, publicKey[1], t)

        println("Timing decryption of $(plaintext)")
        @time decryptedtext = decrypt(cyphertext,
            publicKey[1],
            privateKey[1],
            privateKey[2],
            privateKey[3],
            privateKey[4],
            m, t)

        println("decrypted == plaintext? $(decryptedtext == plaintext)")
    end
end

GPubLoc = "./keys/public_g.txt"
GPrivLoc = "./keys/private_g.txt"
tLoc = "./keys/public_t.txt"
mLoc = "./keys/public_m.txt"
SLoc = "./keys/private_s.txt"
PLoc = "./keys/private_p.txt"
supLoc = "./keys/private_support.txt"
goppaLoc = "./keys/private_polynomial.txt"

parameters = [
    [2048, 31, 11],
    # [2048, 51, 11],
    # [2048, 81, 11],
    # [4096, 41, 12],
    # [4096, 62, 12],
    # [4096, 101, 12],
]

directories = [
    [".\\keys\\2048-1\\public", ".\\keys\\2048-1\\private"],
    # [".\\keys\\2048-2\\public", ".\\keys\\2048-2\\private"],
    # [".\\keys\\2048-3\\public", ".\\keys\\2048-3\\private"],
    # [".\\keys\\4096-1\\public", ".\\keys\\4096-1\\private"],
    # [".\\keys\\4096-2\\public", ".\\keys\\4096-2\\private"],
    # [".\\keys\\4096-3\\public", ".\\keys\\4096-3\\private"],
]

teststrings = [
    "Hello!",
    "Where are we meeting?",
    "I was wondering when you would show up."
]

for ((n, t, m), d) in zip(parameters, directories)
    println("Timing generation of keys with n=$(n), t=$(t), m=$(m)")
    testKeyGeneration(n, t, m, teststrings, d[1], d[2])
end

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
#
# The types and constructors for derivative operators.
#
# The type of the operator objects is DerivativeOperator. This is a
# dumb data structure. It is initialized by the constructors for
# CenteredDifference and UpwindDifference, which call the calculate_weights
# routine in fornberg.jl to generate the stencils. The derivatives
# are computed by the mul! methods defined in
# derivative_operator_functions.jl.
#

index(i::Int, N::Int) = i + div(N, 2) + 1

struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2,S3<:SArray,T3,F} <: AbstractDerivativeOperator{T}
    derivative_order        :: Int
    approximation_order     :: Int
    dx                      :: T2
    len                     :: Int
    stencil_length          :: Int
    stencil_coefs           :: S1
    boundary_stencil_length :: Int
    boundary_point_count    :: Int
    low_boundary_coefs      :: S2
    high_boundary_coefs     :: S3
    offside                 :: Int
    coefficients            :: T3
    coeff_func              :: F
end

struct nonlinear_diffusion!{N} end
struct nonlinear_diffusion{N} end

function nonlinear_diffusion!{N}(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int,
    p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
    nknots::Int) where {T<:Real, N}
    #q is given by bc*u , u being the unknown function
    #p is given as some function of q , p being the diffusion coefficient

    @assert approx_order>1 "approximation_order must be greater than 1."
    if first_differential_order > 0
    du .= (CenteredDifference{N}(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference{N}(second_differential_order,approx_order,dx,nknots)*p)
    else
    du .= q[2:(nknots + 1)].*(CenteredDifference{N}(second_differential_order,approx_order,dx,nknots)*p)
    end

    for l = 1:(second_differential_order - 1)
    du .= du .+ binomial(second_differential_order,l)*(CenteredDifference{N}(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference{N}(second_differential_order - l,approx_order,dx,nknots)*p)
    end

    du .= du .+ (CenteredDifference{N}(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]

end

# An out of place workaround for the mutating version
function nonlinear_diffusion{N}(second_differential_order::Int, first_differential_order::Int, approx_order::Int,
    p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real},
    nknots::Int) where {T<:Real, N}

    du = similar(q,length(q) - 2)
    return nonlinear_diffusion!{N}(du,second_differential_order,first_differential_order,approx_order,p,q,dx,nknots)
end

struct CenteredDifference{N} end

function CenteredDifference{N}(derivative_order::Int,
                            approximation_order::Int, dx::T,
                            len::Int, coeff_func=nothing) where {T<:Real,N}
    @assert approximation_order>1 "approximation_order must be greater than 1."
    stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
    boundary_stencil_length = derivative_order + approximation_order
    dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)
    left_boundary_x         = 0:(boundary_stencil_length-1)
    right_boundary_x        = reverse(-boundary_stencil_length+1:0)

    boundary_point_count    = div(stencil_length,2) - 1 # -1 due to the ghost point
    # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
    deriv_spots             = (-div(stencil_length,2)+1) : -1  # unused
    L_boundary_deriv_spots  = left_boundary_x[2:div(stencil_length,2)]
    R_boundary_deriv_spots  = right_boundary_x[2:div(stencil_length,2)]

    stencil_coefs           = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x))
    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, left_boundary_x)) for x0 in L_boundary_deriv_spots]
    low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)

    # _high_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, reverse(right_boundary_x))) for x0 in R_boundary_deriv_spots]
    # high_boundary_coefs      = convert(SVector{boundary_point_count},reverse(_high_boundary_coefs))

    high_boundary_coefs      = convert(SVector{boundary_point_count},reverse(map(reverse, _low_boundary_coefs*(-1)^derivative_order)))

    offside = 0

    coefficients            = coeff_func isa Nothing ? nothing : fill!(Vector{T}(undef,len),0)

    DerivativeOperator{T,N,false,T,typeof(stencil_coefs),
        typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
        typeof(coeff_func)}(
        derivative_order, approximation_order, dx, len, stencil_length,
        stencil_coefs,
        boundary_stencil_length,
        boundary_point_count,
        low_boundary_coefs,
        high_boundary_coefs,offside,coefficients,coeff_func
        )
end

function generate_coordinates(i::Int, stencil_x, dummy_x, dx::AbstractVector{T}) where T <: Real
    len = length(stencil_x)
    stencil_x .= stencil_x.*zero(T)
    for idx in 1:div(len,2)
        shifted_idx1 = index(idx, len)
        shifted_idx2 = index(-idx, len)
        stencil_x[shifted_idx1] = stencil_x[shifted_idx1-1] + dx[i+idx-1]
        stencil_x[shifted_idx2] = stencil_x[shifted_idx2+1] - dx[i-idx]
    end
    return stencil_x
end

function CenteredDifference{N}(derivative_order::Int,
                            approximation_order::Int, dx::AbstractVector{T},
                            len::Int, coeff_func=nothing) where {T<:Real,N}

    stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
    boundary_stencil_length = derivative_order + approximation_order
    stencil_x               = zeros(T, stencil_length)
    boundary_point_count    = div(stencil_length,2) - 1# -1 due to the ghost point

    interior_x              = boundary_point_count+2:len+1-boundary_point_count
    dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)-1
    low_boundary_x          = [zero(T); cumsum(dx[1:boundary_stencil_length-1])]
    high_boundary_x         = cumsum(dx[end-boundary_stencil_length+1:end])
    # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
    deriv_spots             = (-div(stencil_length,2)+1) : -1

    stencil_coefs           = [convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), generate_coordinates(i, stencil_x, dummy_x, dx))) for i in interior_x]
    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
                                                                  calculate_weights(derivative_order, low_boundary_x[i+1], low_boundary_x)) for i in 1:boundary_point_count]
    low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T},
                                                                  calculate_weights(derivative_order, high_boundary_x[end-i], high_boundary_x)) for i in boundary_point_count:-1:1]
    high_boundary_coefs      = convert(SVector{boundary_point_count},_high_boundary_coefs)

    offside = 0

    coefficients            = coeff_func isa Nothing ? nothing : zeros(T,len)

    DerivativeOperator{T,N,false,typeof(dx),typeof(stencil_coefs),
        typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
        typeof(coeff_func)}(
        derivative_order, approximation_order, dx,
        len, stencil_length,
        stencil_coefs,
        boundary_stencil_length,
        boundary_point_count,
        low_boundary_coefs,
        high_boundary_coefs,offside,coefficients,coeff_func
        )
end

struct UpwindDifference{N} end

"""
```
UpwindDifference{N}(derivative_order, approximation_order, dx, len, coeff_func = nothing)
```
constructs a DerivativeOperator that automatically implements upwinding.

### Inputs
* `dx::T` or `dx::Vector{T}`: grid spacing
* `coeff_func`: function mapping index in the grid to coefficient at that grid location

### Examples
julia> drift = [1., 1., -1.]
julia> L1 = UpwindDifference(1, 1, 1., 3, i -> drift[i])
julia> L2 = UpwindDifference(1, 1, 1., 3, i -> 1.)
julia> Q = Neumann0BC(1, 1.)
julia> Array(L1 * Q)[1]
3×3 Array{Float64,2}:
 -1.0   1.0   0.0
  0.0  -1.0   1.0
  0.0   1.0  -1.0
julia> Array(L2 * Q)[1]
3×3 Array{Float64,2}:
 -1.0   1.0  0.0
  0.0  -1.0  1.0
  0.0   0.0  0.0

"""
function UpwindDifference{N}(derivative_order::Int,
                             approximation_order::Int, dx::T,
                             len::Int, coeff_func=nothing; offside::Int=0) where {T<:Real,N}

    @assert offside > -1 "Number of offside points should be non-negative"
    @assert offside <= div(derivative_order + approximation_order - 1,2) "Number of offside points should not exceed the primary wind points"

    stencil_length          = derivative_order + approximation_order
    boundary_stencil_length = derivative_order + approximation_order
    boundary_point_count    = boundary_stencil_length - 2 - offside

    # TODO: Clean up the implementation here so that it is more readable and easier to extend in the future
    dummy_x = (0.0 - offside) : stencil_length - 1.0 - offside
    stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, 0.0, dummy_x))

    low_boundary_x         = 0.0:(boundary_stencil_length-1)
    L_boundary_deriv_spots = 1.0:boundary_stencil_length - 2.0 - offside
    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots]
    low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)

    high_boundary_x         = 0.0:-1.0:-(boundary_stencil_length-1.0)
    R_boundary_deriv_spots = -1.0:-1.0:-(boundary_stencil_length-2.0)
    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots]
    high_boundary_coefs = convert(SVector{boundary_point_count + offside},_high_boundary_coefs)

    coefficients = zeros(T,len)
    if coeff_func != nothing
        compute_coeffs!(coeff_func, coefficients)
    end

    DerivativeOperator{T,N,true,T,typeof(stencil_coefs),
        typeof(low_boundary_coefs),typeof(high_boundary_coefs),Vector{T},
        typeof(coeff_func)}(
        derivative_order, approximation_order, dx, len, stencil_length,
        stencil_coefs,
        boundary_stencil_length,
        boundary_point_count,
        low_boundary_coefs,
        high_boundary_coefs,offside,coefficients,coeff_func
        )
end

# TODO implement the non-uniform grid
function UpwindDifference{N}(derivative_order::Int,
                          approximation_order::Int, dx::AbstractVector{T},
                          len::Int, coeff_func=nothing; offside::Int=0) where {T<:Real,N}

    @assert offside > -1 "Number of offside points should be non-negative"
    @assert offside <= div(derivative_order + approximation_order - 1,2) "Number of offside points should not exceed the primary wind points"

    stencil_length          = derivative_order + approximation_order
    boundary_stencil_length = derivative_order + approximation_order
    boundary_point_count    = boundary_stencil_length - 2 - offside

    # Compute Stencils
    # Compute grid from dx
    x = [0.0, cumsum(dx)...]

    # compute low_boundary_coefs: low_boundary_coefs[upwind = 1 downwind = 2, index of point]
    if offside == 0
        _upwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i+1:i+boundary_stencil_length])) for i in 1:boundary_point_count])
    else
        _upwind_coefs = SMatrix{1,boundary_point_count}(append!([convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, x[i+1], x[1:boundary_stencil_length])) for i in 1:offside-1],[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, x[i+1], x[i+1-offside:i+boundary_stencil_length-offside])) for i in offside:boundary_point_count]))
    end
    _downwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[1:boundary_stencil_length])) for i in 1:boundary_point_count])
    low_boundary_coefs = [_upwind_coefs ; _downwind_coefs]

    # compute stencil_coefs: low_boundary_coefs[upwind = 1 downwind = 2, index of point]
    _upwind_coefs = SMatrix{1,len - 2*boundary_point_count}([convert(SVector{stencil_length, T}, calculate_weights(derivative_order, x[i+1], x[i+1-offside:i+stencil_length-offside])) for i in boundary_point_count+1:len-boundary_point_count])
    _downwind_coefs = SMatrix{1,len - 2*boundary_point_count}([convert(SVector{stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2+offside:i+1+offside])) for i in boundary_point_count+1:len-boundary_point_count])
    stencil_coefs = [_upwind_coefs ; _downwind_coefs]

    # compute high_boundary_coefs: low_boundary_coefs[upwind = 1 downwind = 2, index of point]
    _upwind_coefs = SMatrix{1,boundary_point_count + offside}([convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, x[i+1], x[len-boundary_stencil_length+3:len+2])) for i in len-boundary_point_count+1-offside:len])
    if offside == 0
        _downwind_coefs = SMatrix{1,boundary_point_count}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2:i+1])) for i in len-boundary_point_count+1:len])
    elseif offside == 1
        _downwind_coefs = SMatrix{1,boundary_point_count + offside}([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2+offside:i+1+offside])) for i in len-boundary_point_count+1-offside:len-offside+1])
    else
        _downwind_coefs = SMatrix{1,boundary_point_count + offside}(append!([convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[i-stencil_length+2+offside:i+1+offside])) for i in len-boundary_point_count+1-offside:len-offside+1],[convert(SVector{boundary_stencil_length,T}, calculate_weights(derivative_order, x[i+1], x[len-stencil_length+3:len+2])) for i in len-offside+2:len]))
    end
    high_boundary_coefs = [_upwind_coefs ; _downwind_coefs]

    # Compute coefficients
    coefficients = zeros(T,len)
    if coeff_func != nothing
        compute_coeffs!(coeff_func, coefficients)
    end

    DerivativeOperator{T,N,true,typeof(dx),typeof(stencil_coefs),
        typeof(low_boundary_coefs),typeof(high_boundary_coefs),Vector{T},
        typeof(coeff_func)}(
        derivative_order, approximation_order, dx, len, stencil_length,
        stencil_coefs,
        boundary_stencil_length,
        boundary_point_count,
        low_boundary_coefs,
        high_boundary_coefs,offside,coefficients,coeff_func
        )
end

CenteredDifference(args...) = CenteredDifference{1}(args...)
UpwindDifference(args...;kwargs...) = UpwindDifference{1}(args...;kwargs...)
nonlinear_diffusion(args...) = nonlinear_diffusion{1}(args...)
nonlinear_diffusion!(args...) = nonlinear_diffusion!{1}(args...)
use_winding(A::DerivativeOperator{T,N,Wind}) where {T,N,Wind} = Wind
diff_axis(A::DerivativeOperator{T,N}) where {T,N} = N
function ==(A1::DerivativeOperator, A2::DerivativeOperator)
    return all([eval(:($A1.$name == $A2.$name)) for name in fieldnames(DerivativeOperator)])
end
function Laplacian(aor::Int, dxyz::Union{NTuple{N, T}, NTuple{N,AbstractVector{T}}}, s::NTuple{N,I}, coeff_func=nothing) where {T,N,I<:Int}
    return sum(CenteredDifference{i}(2, aor, dxyz[i], s[i], coeff_func) for i in 1:N)
end