text:
module GateMap
# =========================================================================== #
#                           Imported modules
# =========================================================================== #

using Base
using LinearAlgebra
using StaticArrays
import QXZoo.GateOps
import QXZoo.DefaultGates

# =========================================================================== #
#                           Exported functions
# =========================================================================== #

export create_gate_1q, create_gate_2q

# =========================================================================== #
#                       Gate caching and retrieval
# =========================================================================== #

"""Module gate cache

Since many different circuits may use the same gates, keeping a module-global
cache makes sense to avoid recreating them. Each subcircuit can have a subset of
the global cache's gates."""
gates = Dict{Union{GateOps.GateSymbol, GateOps.GateSymbolP, Function}, Function}()

# TODO
# const gates = let .... end

"""
    init_cache()

Initialise the cache for default gates used in circuit generation. Modules generating
gates not stored here should cache them for use by other circuits.

"""
function init_cache() #gates::Dict{Union{<:GateOps.AGateSymbol, Function}, Function})
    push!(gates, DefaultGates.GateSymbols.p00 => p00)
    push!(gates, DefaultGates.GateSymbols.p10 => p10)
    push!(gates, DefaultGates.GateSymbols.p01 => p01)
    push!(gates, DefaultGates.GateSymbols.p11 => p11)
    push!(gates, DefaultGates.GateSymbols.I => I)
    push!(gates, DefaultGates.GateSymbols.h => h)
    push!(gates, DefaultGates.GateSymbols.x => x)
    push!(gates, DefaultGates.GateSymbols.y => y)
    push!(gates, DefaultGates.GateSymbols.z => z)

    push!(gates, DefaultGates.GateSymbols.c_x => c_x)
    push!(gates, DefaultGates.GateSymbols.c_y => c_y)
    push!(gates, DefaultGates.GateSymbols.c_z => c_z)
    push!(gates, DefaultGates.GateSymbols.r_x => r_x)
    push!(gates, DefaultGates.GateSymbols.r_y => r_y)
    push!(gates, DefaultGates.GateSymbols.r_z => r_z)
    push!(gates, DefaultGates.GateSymbols.r_phase => r_phase)

    push!(gates, DefaultGates.GateSymbols.c_r_x => c_r_x)
    push!(gates, DefaultGates.GateSymbols.c_r_y => c_r_y)
    push!(gates, DefaultGates.GateSymbols.c_r_z => c_r_z)
    push!(gates, DefaultGates.GateSymbols.c_r_phase => c_r_phase)

    push!(gates, DefaultGates.GateSymbols.s => s )
    push!(gates, DefaultGates.GateSymbols.t => t )
    return ;
end

# =========================================================================== #
#                       Custom gate storage: 1 qubit
# =========================================================================== #

"""
    create_gate_1q(gate_label::String, gen_func::Function)

Creates a new user-defined gate, caches the generating function, and returns the
GateOps.GateSymbol key for use in a circuit.

# Examples
```julia-repl

```
"""
function create_gate_1q(gate_label::String, gen_func::Function)
    mat2x2 = nothing
    if applicable(gen_func, pi/2)
        mat2x2 = gen_func(pi/2)
    else
        mat2x2 = gen_func()
    end
    if size(mat2x2) != (2,2)
        error("Please ensure generated matrix is 2x2")
    end
    gs = GateOps.GateSymbol(Symbol(gate_label))
    cache_gate!(gs, gen_func)
    return gs
end

"""
    create_gate_1q(gate_label::String, mat::SArray{Tuple{2,2},Complex{Float64},2,4})

Creates a new user-defined gate, caches the generating matrix, and returns the
GateOps.GateSymbol key for use in a circuit.

# Examples
```julia-repl

```
"""
function create_gate_1q(gate_label::String, mat::SArray{Tuple{2,2},Complex{Float64},2,4})
    return create_gate_1q(gate_label, ()->mat)
end

"""
    create_gate_1q(gate_label::String, mat::Array{Complex{Float64},2})

Creates a new user-defined gate, caches the generating matrix, and returns the
GateOps.GateSymbol key for use in a circuit. Matrix is converted to 2x2 StaticArrays
SMatrix internally.

# Examples
```julia-repl

```
"""
function create_gate_1q(gate_label::String, mat::Array{<:Number,2})
    return create_gate_1q(gate_label, SArray{Tuple{2,2},Complex{Float64},2,4}(mat))
end

# =========================================================================== #
#                       Custom gate storage: 2 qubit
# =========================================================================== #

"""
    create_gate_2q(gate_label::String, gen_func::Function)

Creates a new user-defined gate, caches the generating function, and returns the
GateOps.GateSymbol key for use in a circuit.

# Examples
```julia-repl

```
"""
function create_gate_2q(gate_label::String, gen_func::Function)
    mat4x4 = nothing
    if applicable(gen_func, pi/2)
        mat4x4 = gen_func(pi/2)
    else
        mat4x4 = gen_func()
    end
    if size(mat4x4) != (4,4)
        error("Please ensure generated matrix is 4x4")
    end
    gs = GateOps.GateSymbol(Symbol(gate_label))
    cache_gate!(gs, gen_func)
    return gs
end

"""
    create_gate_2q(gate_label::String, mat::SArray{Tuple{2,2},Complex{Float64},2,4})

Creates a new user-defined gate, caches the generating matrix, and returns the
GateOps.GateSymbol key for use in a circuit.

# Examples
```julia-repl

```
"""
function create_gate_2q(gate_label::String, mat::SArray{Tuple{4,4},Complex{Float64},2,16})
    return create_gate_2q(gate_label, ()->mat)
end
"""
    create_gate_2q(gate_label::String, mat::Array{Complex{Float64},2})

Creates a new user-defined gate, caches the generating matrix, and returns the
GateOps.GateSymbol key for use in a circuit. Matrix is converted to 2x2 StaticArrays
SMatrix internally.

# Examples
```julia-repl

```
"""
function create_gate_2q(gate_label::String, mat::Array{<:Number,2})
    return create_gate_2q(gate_label, SArray{Tuple{4,4},Complex{Float64},2,16}(mat))
end

# =========================================================================== #

"""
    cache_gate!(label::GateOps.AGateSymbol, mat::Matrix{<:Number})

Adds a mapping between label=>mat for fast retrieval of gates in circuit generation.

# Examples
```julia-repl
julia> QXZoo.GateMap.cache_gate!(QXZoo.GateOps.GateSymbol(:mygate), ()->[ 1 0; 0 1])
```
"""
function cache_gate!(key::Union{GateOps.GateSymbol, GateOps.GateSymbolP, Function}, mat_func::Function)
    if ~haskey(gates, key)
        gates[key] = mat_func
    end
    return ;
end

"""
    replace_gate!(label::GateOps.AGateSymbol, mat::Matrix{<:Number})

Replaces existing gate mapping as defined with cache_gate! otherwise performs no-op.

# Examples
```julia-repl
julia> QXZoo.GateMap.replace_gate!(QXZoo.GateOps.GateSymbol(:mygate), ()->[ 1 0; 0 -1])
```
"""
function replace_gate!(key::Union{GateOps.GateSymbol, GateOps.GateSymbolP, Function}, mat_func::Function)
    if haskey(gates, key)
        gates[key] = mat_func
    end
    return ;
end

# =========================================================================== #
#                       Projection operators
# =========================================================================== #
"""
    p00()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the |0><0| projector gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.p00()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
```
"""
function p00()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [1.0+0.0im 0.0; 0.0 0.0]
end

"""
    p01()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the |0><1| projector gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.p01()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.0+0.0im  1.0+0.0im
 0.0+0.0im  0.0+0.0im
```
"""
function p01()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [0.0 1.0+0.0im; 0.0 0.0]
end

"""
    p10()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the |1><0| projector gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.p10()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.0+0.0im  0.0+0.0im
 1.0+0.0im  0.0+0.0im
```
"""
function p10()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [0.0 0.0; 1.0+0.0im 0.0]
end

"""
    p11()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the |1><1| projector gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.p11()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im
```
"""
function p11()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [0.0 0.0; 0.0 1.0+0.0im]
end

# =========================================================================== #
#                       Pauli operators
# =========================================================================== #
"""
    x()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the Pauli-X gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.x()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.0+0.0im  1.0+0.0im
 1.0+0.0im  0.0+0.0im
```
"""
function x()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [0.0 1.0+0.0im; 1.0+0.0im 0.0]
end

"""
    y()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the Pauli-Y gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.y()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.0+0.0im  -0.0-1.0im
 0.0+1.0im   0.0+0.0im
```
"""
function y()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [0.0 -1.0im; 1.0im 0.0]
end

"""
    z()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the Pauli-Z gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.z()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im   0.0+0.0im
 0.0+0.0im  -1.0+0.0im
```
"""
function z()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return @SMatrix [1.0+0.0im 0.0; 0.0 -1.0+0.0im]
end

# =========================================================================== #
#                       Additional operators
# =========================================================================== #

"""
    h()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the Hadamard gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.h()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.707107+0.0im   0.707107+0.0im
 0.707107+0.0im  -0.707107+0.0im
```
"""
function h()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return SMatrix{2,2}([1.0+0.0im 1.0+0.0im; 1.0+0.0im -1.0+0.0im].*(1/sqrt(2)))
end

"""
    I()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the 2x2 identity gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.I()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im
```
"""
function I()::Diagonal{Complex{Float64},SArray{Tuple{2},Complex{Float64},1,2}}
    return StaticArrays.SDiagonal(1.0+0im,1.0+0im)
end

"""
    s()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the S gate (sqrt(Z)).

# Examples
```julia-repl
julia> QXZoo.GateMap.s()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+1.0im
```
"""
function s()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return sqrt(z())
end

"""
    t()::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the T gate (sqrt(S)).

# Examples
```julia-repl
julia> QXZoo.GateMap.t()
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im       0.0+0.0im
 0.0+0.0im  0.707107+0.707107im
```
"""
function t()::SArray{Tuple{2,2},Complex{Float64},2,4}
    return sqrt(s())
end

# =========================================================================== #
#                       Rotation operators
# =========================================================================== #
"""
    r_x(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the parametric r_x gate: exp(-0.5im*θ*x())

# Examples
```julia-repl
julia> QXZoo.GateMap.r_x(pi/3)
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.866025+0.0im       0.0-0.5im
      0.0-0.5im  0.866025+0.0im
```
"""
function r_x(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}
    return exp(-1im*x()*θ/2)
end
"""
    r_y(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the parametric r_y gate: exp(-0.5im*θ*y())

# Examples
```julia-repl
julia> QXZoo.GateMap.r_y(pi/6)
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.965926+0.0im  -0.258819+0.0im
 0.258819+0.0im   0.965926+0.0im
```
"""
function r_y(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}
    return exp(-1im*y()*θ/2)
end
"""
    r_z(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the parametric r_z gate: exp(-0.5im*θ*z())

# Examples
```julia-repl
julia> QXZoo.GateMap.r_z(pi/2)
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 0.707107-0.707107im       0.0+0.0im
      0.0+0.0im       0.707107+0.707107im
```
"""
function r_z(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}
    return exp(-1im*z()*θ/2)
end
"""
    r_phase(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}

Function which generates the parametric phase shifting gate [1 0; 0 exp(1im*θ)].
Note, this is not the Phase gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.r_phase(pi/4)
2×2 StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4} with indices SOneTo(2)×SOneTo(2):
 1.0+0.0im       0.0+0.0im
 0.0+0.0im  0.707107+0.707107im
```
"""
function r_phase(θ::Number)::SArray{Tuple{2,2},Complex{Float64},2,4}
    return [1 0; 0 exp(1im*θ)]
end

# =========================================================================== #
#                       Controlled Pauli operators
# =========================================================================== #
"""
    c_x()::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled Pauli-X gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_x()
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im
```
"""
function c_x()::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), x())
end

"""
    c_y()::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled Pauli-X gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_y()
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0-1.0im
 0.0+0.0im  0.0+0.0im  0.0+1.0im  0.0+0.0im
```
"""
function c_y()::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), y())
end

"""
    c_z()::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled Pauli-X gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_z()
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im  0.0+0.0im   0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im   0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im   0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  -1.0+0.0im
```
"""
function c_z()::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), z())
end

# =========================================================================== #
#                       Controlled rotation operators
# =========================================================================== #
"""
    c_r_x(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled parametric r_x gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_r_x(pi/3)
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im       0.0+0.0im       0.0+0.0im
 0.0+0.0im  1.0+0.0im       0.0+0.0im       0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.866025+0.0im       0.0-0.5im
 0.0+0.0im  0.0+0.0im       0.0-0.5im  0.866025+0.0im
```
"""
function c_r_x(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), r_x(θ))
end

"""
    c_r_y(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled parametric r_y gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_r_y(pi/4)
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im       0.0+0.0im        0.0+0.0im
 0.0+0.0im  1.0+0.0im       0.0+0.0im        0.0+0.0im
 0.0+0.0im  0.0+0.0im   0.92388+0.0im  -0.382683+0.0im
 0.0+0.0im  0.0+0.0im  0.382683+0.0im    0.92388+0.0im
```
"""
function c_r_y(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), r_y(θ))
end
"""
    c_r_z(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled parametric r_y gate as a block diagonal 4x4 matrix.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_r_z(pi/6)
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im       0.0+0.0im            0.0+0.0im
 0.0+0.0im  1.0+0.0im       0.0+0.0im            0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.965926-0.258819im       0.0+0.0im
 0.0+0.0im  0.0+0.0im       0.0+0.0im       0.965926+0.258819im
```
"""
function c_r_z(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), r_z(θ))
end

"""
    c_r_phase(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}

Function which generates the controlled parametric phase shifting gate [1 0; 0 exp(1im*θ)] as a block diagonal 4x4 matrix.

Note, this is not the controlled Phase gate.

# Examples
```julia-repl
julia> QXZoo.GateMap.c_r_phase(pi/9)
4×4 StaticArrays.SArray{Tuple{4,4},Complex{Float64},2,16} with indices SOneTo(4)×SOneTo(4):
 1.0+0.0im  0.0+0.0im  0.0+0.0im       0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im       0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im       0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.939693+0.34202im
```
"""
function c_r_phase(θ::Number)::SArray{Tuple{4,4},Complex{Float64},2,16}
    return kron(p00(), I()) + kron(p11(), r_phase(θ))
end

# =========================================================================== #
"""
Allow access directly using a gate-call for ease-of-use
"""
function Base.getindex(gc::Dict{GateOps.AGateSymbol, Matrix{<:Number}}, gs::GateOps.AGateCall)
    return gc[gs.gate_label]
end

"""
    clear_cache()

Empty the stored gates from the cache.
"""
function clear_cache()
    empty!(gates)
end

end


text:
#Another auotmated testing package. First time trying it out

julia> using Checkers
julia> f(x) = x^2
julia> @test_formany -10<x<10, f(x) >= 0
Test Passed
  Expression: (:((-10 < x < 10, f(x) >= 0)), :(mode = test_formany))

julia> @test_forall x in -10:10, f(x) >= 0
Test Passed
  Expression: (x in -10:10, f(x) >= 0)
    @test_forall x in 0:10, f(x) >= 0

julia> @test_forall x in -10:10, f(x) > 0 #Should fail b/c f(0) = 0
Test Failed
  Expression: (x in -10:10, f(x) > 0)
ERROR: There was an error during testing

ext: jl


text:
####
#### Defines list of how-to-guides
####

extending_clima = Any[
    "Home" => "ExtendingCLIMA/Contributing/CONTRIBUTING.md",
    "Arrays" => Any[
    # "Home" => "ExtendingCLIMA/Arrays/index.md"
    ],
    "Common" => Any[
        "SurfaceFluxes" => "ExtendingCLIMA/Common/SurfaceFluxes.md",
        "MoistThermodynamics" => "ExtendingCLIMA/Common/MoistThermodynamics.md",
    ],
    "Atmos" => Any[
        "AtmosModel" => "ExtendingCLIMA/Atmos/AtmosModel.md",
        "Microphysics" => "ExtendingCLIMA/Atmos/Microphysics.md",
        "Turbulence" => "ExtendingCLIMA/Atmos/Model/turbulence.md",
        "Tracers" => "ExtendingCLIMA/Atmos/Model/tracers.md",
    ],
    "Ocean" => Any[
    # "Home" => "ExtendingCLIMA/Ocean/index.md"
    ],
    "Land" => Any[
    # "Home" => "ExtendingCLIMA/Land/index.md"
    ],
    "Diagnostics" => Any["Variables" => "ExtendingCLIMA/Diagnostics/DiagnosticVariables.md"],
    "Numerics" => Any[
        "Meshes" => Any[
        # "Home" => "ExtendingCLIMA/Numerics/Meshes/index.md",
        ],
        "DG methods" => Any["How to make a balance law" => "ExtendingCLIMA/Numerics/DGmethods/how_to_make_a_balance_law.md",],
        "ODE Solvers" => Any["Time-integration" => "ExtendingCLIMA/Numerics/ODESolvers/Timestepping.md",],
        "Linear Solvers" => Any["Iterative Solvers" => "ExtendingCLIMA/Numerics/LinearSolvers/IterativeSolvers.md",],
    ],
]


text:
module CoulombExchange

#LearnHatreeFock.jl Modules
using TypesBasis,TypesParticles

export getcoulexchange

@doc raw"""
function getcoulexchange()

"""
function getcoulexchange(density::Array,elecelecrepul::Array)

    nbasisfunc = size(density)[1];
    coulexchange = zeros(Float64,nbasisfunc,nbasisfunc);

    for l=1:nbasisfunc
        for m=1:nbasisfunc
            for n=1:nbasisfunc
                for o=1:nbasisfunc
                    diff = 2.00e0*elecelecrepul[l,m,n,o]-elecelecrepul[l,n,m,o];
                    coulexchange[l,m] += density[n,o]*(diff);
                end
            end #n=1:nbasisfunc
        end
    end #l=1:nbasisfunc

    return coulexchange
end #getcoulexchange

end #CoulombExchange


text:
## This downloads and populates the spec test files

using Mustache
using YAML
using Test

ghub = "https://raw.githubusercontent.com/mustache/spec/72233f3ffda9e33915fd3022d0a9ebbcce265acd/specs/{{:spec}}.yml"

specs = ["comments",
          "delimiters",
          "interpolation",
          "inverted",
          "partials",
          "sections"#,
          #"~lambdas"
          ]

D = Dict()
for spec in specs
    nm = Mustache.render(ghub, spec=spec)
    D[spec] = YAML.load_file(download(nm))
end

function write_spec_file(fname)
    io = open("spec_$fname.jl","w")
    println(io, """
using Mustache
using Test
""")

    x = D[fname]
    println(io, "@testset \" $fname \" begin\n")
    for (i,t) in enumerate(x["tests"])
        tpl = t["template"]
        tpl = replace(tpl, "\\" => "\\\\")
        tpl = replace(tpl, r"\"" => "\\\"")
        data = t["data"]
        desc = t["desc"]

        if haskey(t, "partials")
            for (k,v) in t["partials"]
                data[k] = v
            end
        end

        val = try
            Mustache.render(t["template"], data)
        catch err
            "Failed"
        end

        expected = t["expected"]
        expected_ = replace(expected, "\\" => "\\\\")
        expected_ = replace(expected_, r"\"" => "\\\"")
        expected_ = "\"\"\"" * expected_ * "\"\"\""

        println(io, "\n\t## $desc")
        print(io, "tpl = \"\"\"")
        print(io, tpl)
        println(io, "\"\"\"\n")
        println(io, "\t@test Mustache.render(tpl, $data) == $expected_")
        println("")
    end
    println(io, "end\n\n")
    close(io)
end

for spec in specs
    write_spec_file(spec)
end

# partials are different, as they refer to an external file
# tihs should clean up temp files, but we don't run as part of test suite
# test 7 fails, but I think that one is wrong
using Test
function test_partials()
    for spec in specs
        for (i,t) in enumerate(D[spec]["tests"])
            if haskey(t, "partials")
                println("""Test $spec / $i""")

                d = t["data"]
                partial = t["partials"]
                for (k,v) in partial
                    io = open(k, "w")
                    write(io, v)
                    close(io)
                end

                tpl = t["template"]
                expected = t["expected"]
                out =  Mustache.render(tpl, d)

                val = out == expected
                if val
                    @test val
                else
                    val = replace(out, r"\n"=>"") == replace(expected, r"\n"=>"")
                    if val
                        println("""$(t["desc"]): newline issue ...""")
                        @test val
                    else
                        println("""$(t["desc"]): FAILED:\n $out != $expected""")
                    end
                end

                for (k,v) in partial
                    rm(k)
                end
            end
        end
    end
end


text:
"""
Special matrix
"""

"""
eye(n, ::Type{Tv} = Float64)::Matrix{Tv}

Make an indentity matrix
"""

function eye(n, ::Type{Tv} = Float64)::Matrix{Tv} where {Tv}
    m = zeros(Tv, n,n)
    @inbounds for i = 1:n
        m[i,i] = Tv(1)
    end
    m
end


text:
import Zygote
import RestrictedBoltzmannMachines as RBMs
using Test: @test, @testset
using Random: bitrand
using Statistics: mean
using LinearAlgebra: norm

@testset "∂free_energy(rbm, v)" begin
    rbm = RBMs.BinaryRBM(randn(5,2), randn(4,3), randn(5,2,4,3))
    v = bitrand(size(rbm.visible)..., 7)
    gs = Zygote.gradient(rbm) do rbm
        mean(RBMs.free_energy(rbm, v))
    end
    ∂F = RBMs.∂free_energy(rbm, v)
    @test ∂F.visible.θ ≈ only(gs).visible.θ
    @test ∂F.hidden.θ ≈ only(gs).hidden.θ
    @test ∂F.w ≈ only(gs).w

    gnorms = RBMs.gradnorms(∂F)
    @test gnorms.visible.θ ≈ norm(∂F.visible.θ)
    @test gnorms.hidden.θ ≈ norm(∂F.hidden.θ)
    @test gnorms.w ≈ norm(∂F.w)

    gλ = RBMs.gradmult(∂F, 2.3)
    @test gλ.visible.θ ≈ ∂F.visible.θ * 2.3
    @test gλ.hidden.θ ≈ ∂F.hidden.θ * 2.3
    @test gλ.w ≈ ∂F.w * 2.3

    wts = rand(7)
    gs = Zygote.gradient(rbm) do rbm
        RBMs.wmean(RBMs.free_energy(rbm, v); wts)
    end
    ∂F = RBMs.∂free_energy(rbm, v; wts)
    @test ∂F.visible.θ ≈ only(gs).visible.θ
    @test ∂F.hidden.θ ≈ only(gs).hidden.θ
    @test ∂F.w ≈ only(gs).w
end

@testset "∂contrastive_divergence" begin
    rbm = RBMs.BinaryRBM(randn(5,2), randn(4,3), randn(5,2,4,3))
    vd = bitrand(size(RBMs.visible(rbm))..., 7)
    vm = bitrand(size(RBMs.visible(rbm))..., 7)
    gs = Zygote.gradient(rbm) do rbm
        RBMs.contrastive_divergence(rbm, vd, vm)
    end
    ∂ = RBMs.∂contrastive_divergence(rbm, vd, vm)
    @test ∂.visible.θ ≈ only(gs).visible.θ
    @test ∂.hidden.θ ≈ only(gs).hidden.θ
    @test ∂.w ≈ only(gs).w
end


text:
# This shows how to

# You will need to download the data from:
# https://www.seismologie.ifg.uni-kiel.de/en/research/research-data/mere2020model.
#
# And make sure that you are in the same directory as the data file

using DelimitedFiles, GeophysicalModelGenerator, GeoStats

# Load data:
data        =   readdlm("El-Sharkawy-etal-G3.2020_MeRE2020_Mediterranean.csv",'|',Float64,'\n', skipstart=23,header=false)
lat         =   data[:,1];
lon         =   data[:,2];
depth       =  -data[:,3];
Vs          =   data[:,4];

# Create 3D regular grid:
Depth_vec       =   unique(depth)
Lon,Lat,Depth   =   LonLatDepthGrid(-10:0.5:40,32:0.25:50,Depth_vec);

# Employ GeoStats to interpolate irregular data points to a regular grid
Cgrid = CartesianGrid((size(Lon,1),size(Lon,2)),(minimum(Lon),minimum(Lat)),(Lon[2,2,2]-Lon[1,1,1],Lat[2,2,2]-Lat[1,1,1]))
Vs_3D = zeros(size(Depth));
for iz=1:size(Depth,3)
    println("Depth = $(Depth[1,1,iz])")
    ind   = findall(x -> x==Depth[1,1,iz], depth)
    coord = PointSet([lon[ind]'; lat[ind]'])
    Geo   = georef((Vs=Vs[ind],), coord)
    P     = EstimationProblem(Geo, Cgrid, :Vs)
    S     = IDW(:Vs => (distance=Euclidean(),neighbors=2));
    sol   = solve(P, S)
    sol_Vs= values(sol).Vs
    Vs_2D = reshape(sol_Vs, size(domain(sol)))
    Vs_3D[:,:,iz] = Vs_2D;
end

# Save data to paraview:
Data_set    =   GeoData(Lon,Lat,Depth,(Vs_km_s=Vs_3D,))
Write_Paraview(Data_set, "MeRe_ElSharkawy")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "BerkeleyDB"
version = v"18.1.40"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://download.oracle.com/berkeley-db/db-$(version).tar.gz", "0cecb2ef0c67b166de93732769abdeba0555086d51de1090df325e18ee8da9c8")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/db-*/docs
# see https://stackoverflow.com/questions/64707079/berkeley-db-make-install-fails-on-linux
mkdir bdb-sql
mkdir gsg_db_server
cd ../dist

# update configure for powerpc
if [[ "${target}" == powerpc64le-* ]]; then
    sed -i s/elf64ppc/elf64lppc/ configure
fi

cd ../build_unix
../dist/configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("db_dump", :db_dump),
    LibraryProduct("libdb", :libdb),
    ExecutableProduct("db_load", :db_load),
    ExecutableProduct("db_stat", :db_stat),
    ExecutableProduct("db_log_verify", :db_log_verify),
    ExecutableProduct("db_archive", :db_archive),
    ExecutableProduct("db_deadlock", :db_deadlock),
    ExecutableProduct("db_hotbackup", :db_hotbackup),
    ExecutableProduct("db_tuner", :db_tuner),
    ExecutableProduct("db_checkpoint", :db_checkpoint),
    ExecutableProduct("db_recover", :db_recover),
    ExecutableProduct("db_replicate", :db_replicate),
    ExecutableProduct("db_printlog", :db_printlog),
    ExecutableProduct("db_verify", :db_verify),
    ExecutableProduct("db_upgrade", :db_upgrade),
    ExecutableProduct("db_convert", :db_convert)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# the data structures for storing variables histories

#=

Varchive should contain:
past, a dictionary indexed by cells,
and points to a VarState

VarState: a dict indexed by name of var,
 with keys being their hashes.
And, it should have ans.

And, Store: a dict from hash to variable.

=#

struct VarState
    vars::Dict
    ans::Any
end

struct Varchive
    past::Dict{Int,VarState}
    store::Dict
end

Varchive() = Varchive(Dict{Int,VarState}(), Dict())

"""
    put_state!(vx::Varchive, key, di::Dict, ansc)

Assuming that state indexed by key (=n) is packed into ijstate,
put it into the Varchive.
"""
function put_state!(vx::Varchive, key, di::Dict, ansc)
    vs = VarState(di, ansc)

    if haskey(vx.past, key)
        error("Key $(key) is already stored.")
    end
    vx.past[key] = vs

    for (name, pair) in di
        h, val = pair
        vs.vars[name] = h
        if !haskey(vx.store, h)
            vx.store[h] = val
        end
    end
end

#=
"""
    put_state_copied!(vx::Varchive, key, di::Dict, ansc)

Assuming that state indexed by key (=n) is packed into ijstate,
put it into the Varchive.
"""
function put_state_copied!(vx::Varchive, key, di::Dict, ansc)
    vs = VarState(di, ansc)

    if haskey(vx.past, key)
        error("Key $(key) is already stored.")
    end
    vx.past[key] = vs

    for (name, pair) in di
        h, copied = pair
        debug_mode && println(IJulia.orig_stdout[],"$name  : $copied")
        vs.vars[name] = h
        if !haskey(vx.store, h)
            vx.store[h] = copied
        end
    end
end
=#

"""
    vars(cell)
    vars(vx::Varchive, cell)

Returns a dictionary of the variables from cell `cell`.
If `vx` is omitted, returns from the default history.
This is how it is usually used.
"""
function vars(vx::Varchive, cell)
    if !haskey(vx.past, cell)
        error("Cell $(cell) was not saved.")
    end

    di = Dict()

    vs = vx.past[cell]

    for (name, h) in vs.vars
        di[name] = vx.store[h]
    end

    return di
end

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


text:
function _add_forces!(
    f::Tuple,
    p,
    idx::Integer,
    jdx::Integer,
    fp::AbstractFloat,
    rp::AbstractFloat,
)
    @inbounds for (k, c) in zip(p, f)
        c[idx] += (fp * k) / rp
        c[jdx] -= (fp * k) / rp
    end
end

function _compute_distance(
    x::SubArray,
    y::SubArray,
    z::SubArray,
    i::Integer,
    j::Integer,
    boxl::Real,
)
    xij = x[i] - x[j]
    yij = y[i] - y[j]
    zij = z[i] - z[j]

    xij -= boxl * round(xij / boxl)
    yij -= boxl * round(yij / boxl)
    zij -= boxl * round(zij / boxl)

    rij2 = xij * xij + yij * yij + zij * zij

    st_positions = @SVector [xij, yij, zij]

    return st_positions, rij2
end

function _compressibilityz(
    pos::StaticVector,
    fp::AbstractFloat,
    rij::AbstractFloat,
    zfactor::ZFactor,
)
    tmp_sum = @. (pos^2.0 * -fp) / rij
    zfactor.zval += sum(tmp_sum)
end

function _compute_energy!(
    positions::AbstractArray,
    forces::AbstractArray,
    params::NamedTuple,
    pot::PairwisePotential;
    rdfobj = nothing,
    zfactor = nothing,
)
    # Initialize necessary variables
    energy = zero(params.rc2)

    x = view(positions, :, 1)
    y = view(positions, :, 2)
    z = view(positions, :, 3)
    fx = view(forces, :, 1)
    fy = view(forces, :, 2)
    fz = view(forces, :, 3)

    for i = 1:params.N-1
        @inbounds @fastmath for j = (i+1):params.N
            (static_positions, rij2) = _compute_distance(x, y, z, i, j, params.boxl)

            if rij2 < params.rc2
                rij2 = √rij2
                u_pair, f_pair = apply!(pot, rij2)
                _add_forces!((fx, fy, fz), static_positions, i, j, f_pair, rij2)
                energy += u_pair
                if !isnothing(rdfobj)
                    simple_rdf!(rdfobj, rij2)
                end
                if !isnothing(zfactor)
                    _compressibilityz(static_positions, f_pair, rij2, zfactor)
                end
            end
        end
    end
    return energy / params.N
end

function energy_force!(
    positions::AbstractArray,
    forces::AbstractArray,
    params::NamedTuple,
    pot::PairwisePotential;
    gofr = nothing,
    zfactor = nothing,
)
    # Retrieve data from the system
    fill!(forces, zero(params.boxl))
    energy =
        _compute_energy!(positions, forces, params, pot; rdfobj = gofr, zfactor = zfactor)
end


text:
"""
    couple_terms(t1, t2)

Generate all possible coupling terms between `t1` and `t2`.  It is
assumed that `t1` and `t2` originate from non-equivalent electrons
(i.e. from _different_ subshells), since the vector model does not
predict correct term couplings for equivalent electrons; some of the
generated terms would violate the Pauli principle; cf. Cowan
p. 108–109.

# Examples

```jldoctest
julia> couple_terms(T"1Po", T"2Se")
1-element Array{Term,1}:
 ²Pᵒ

julia> couple_terms(T"3Po", T"2Se")
2-element Array{Term,1}:
 ²Pᵒ
 ⁴Pᵒ

julia> couple_terms(T"3Po", T"2De")
6-element Array{Term,1}:
 ²Pᵒ
 ²Dᵒ
 ²Fᵒ
 ⁴Pᵒ
 ⁴Dᵒ
 ⁴Fᵒ
```
"""
function couple_terms(t1::Term, t2::Term)
    L1 = t1.L
    L2 = t2.L
    S1 = t1.S
    S2 = t2.S
    p = t1.parity * t2.parity
    sort(vcat([[Term(L, S, p) for S in abs(S1-S2):(S1+S2)]
               for L in abs(L1-L2):(L1+L2)]...))
end

"""
    couple_terms(t1s, t2s)

Generate all coupling between all terms in `t1s` and all terms in
`t2s`.
"""
function couple_terms(t1s::Vector{T}, t2s::Vector{T}) where {T<:Union{Term,<:Real}}
    ts = map(t1s) do t1
        map(t2s) do t2
            couple_terms(t1, t2)
        end
    end
    sort(unique(vcat(vcat(ts...)...)))
end

"""
    final_terms(ts::Vector{<:Vector{<:Union{Term,Real}}})

Generate all possible final terms from the vector of vectors of
individual subshell terms by coupling from left to right.

# Examples

```jldoctest
julia> ts = [[T"1S", T"3S"], [T"2P", T"2D"]]
2-element Array{Array{Term,1},1}:
 [¹S, ³S]
 [²P, ²D]

julia> AtomicLevels.final_terms(ts)
4-element Array{Term,1}:
 ²P
 ²D
 ⁴P
 ⁴D
```
"""
final_terms(ts::Vector{<:Vector{<:T}}) where {T<:Union{Term,Real}} =
    foldl(couple_terms, ts)

couple_terms(J1::T, J2::T) where {T <: Union{Integer,HalfInteger}} =
    collect(abs(J1-J2):(J1+J2))
couple_terms(J1::Real, J2::Real) =
    couple_terms(convert(HalfInteger, J1), convert(HalfInteger, J2))

function intermediate_couplings(its::Vector{T}, t₀::T=zero(T)) where {T<:Union{Term,Integer,HalfInteger}}
    ts = Vector{Vector{T}}()
    for t in couple_terms(t₀, its[1])
        if length(its) == 1
            push!(ts, [t₀, t])
        else
            for ts′ in intermediate_couplings(its[2:end], t)
                push!(ts, vcat(t₀, ts′...))
            end
        end
    end
    ts
end

"""
    intermediate_couplings(its::Vector{IntermediateTerm,Integer,HalfInteger}, t₀ = T"1S")

Generate all intermediate coupling trees from the vector of
intermediate terms `its`, starting from the initial term `t₀`.

# Examples

```jldoctest
julia> intermediate_couplings([IntermediateTerm(T"2S", 1), IntermediateTerm(T"2D", 1)])
2-element Array{Array{Term,1},1}:
 [¹S, ²S, ¹D]
 [¹S, ²S, ³D]
```
"""
intermediate_couplings(its::Vector{<:IntermediateTerm{T}}, t₀::T=zero(T)) where T =
    intermediate_couplings(map(t -> t.term, its), t₀)

"""
    intermediate_couplings(J::Vector{<:Real}, j₀ = 0)

# Examples

```jldoctest
julia> intermediate_couplings([1//2, 3//2])
2-element Array{Array{HalfIntegers.Half{Int64},1},1}:
 [0, 1/2, 1]
 [0, 1/2, 2]
```
"""
intermediate_couplings(J::Vector{IT}, j₀::T=zero(IT)) where {IT <: Real, T <: Real} =
    intermediate_couplings(convert.(HalfInteger, J), convert(HalfInteger, j₀))

export couple_terms, intermediate_couplings


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]
ey = [0.;1.;0.]
ez = [0.;0.;1.]

l1 = 1.0
x, y = .1, .1
b1 = Box(x, l1, y, l1, color = RGBA(1., 1., 0.))

vert11 = [0.;l1 / 2;0.]
vert12 = -vert11

# Initial orientation
phi1, phi2, phi3 = 0., 2pi / 3, 4pi / 3
q1, q2, q3 = Quaternion(RotZ(phi1)), Quaternion(RotZ(phi2)), Quaternion(RotZ(phi3))

# Links
origin = Origin{Float64}()
link1 = Body(b1)
link2 = deepcopy(link1, b1)
link3 = deepcopy(link1, b1)

# Constraints
joint0to123 = EqualityConstraint(Planar(origin, link1, zeros(3), vert12, ez), Planar(origin, link2, zeros(3), vert12, ez), Planar(origin, link3, zeros(3), vert12, ez))
joint1to23 = EqualityConstraint(Spherical(link1, link2, vert11, vert11), Spherical(link1, link3, vert11, vert11))

links = [link1; link2; link3]
constraints = [joint0to123;joint1to23]
shapes = [b1]

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(mech,origin,link1,p2 = vert11,Δq = q1)
setPosition!(mech,origin,link2,p2 = vert11,Δq = q2)
setPosition!(mech,origin,link3,p2 = vert11,Δq = q3)
setForce!(mech,link1,τ = [0;0;2.])
setForce!(mech,link2,τ = [0;0;2.])
setForce!(mech,link3,τ = [0;0;2.])

simulate!(mech,save = true)
visualize!(mech)


text:
"""
`module  JAC.RadiativeAuger`
    ... a submodel of JAC that contains all methods for computing radiative Auger and autoionization amplitudes and rates.
"""
module RadiativeAuger

    using Printf, ..AngularMomentum, ..AtomicState, ..AutoIonization, ..Basics, ..Defaults, ..ManyElectron,
                  ..Radial, ..PhotoEmission, ..TableStrings

    """
    `struct  RadiativeAuger.Settings`  ... defines a type for the settings in estimating radiative-Auger and autoionization rates.

        + multipoles              ::Array{EmMultipole}      ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}         ... Specifies the gauges to be included into the computations.
        + green                   ::Array{GreenChannel,1}   ... Precalculated and user-specified Green function of the ion.
        + NoEnergySharings        ::Int64                   ... Number of energy sharings that are used in the computations for each line.
        + maxKappa                ::Int64                   ... Maximum kappa value of partial waves to be included.
        + printBefore             ::Bool                    ... True, if all energies and lines are printed before their evaluation.
        + operator                ::AbstractEeInteraction   ... Auger operator that is to be used for evaluating the Auger amplitudes:
                                                                allowed values are: CoulombInteraction(), BreitInteraction(), ...
        + lineSelection           ::LineSelection           ... Specifies the selected levels, if any.
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        green                     ::Array{GreenChannel,1}
        NoEnergySharings          ::Int64
        maxKappa                  ::Int64
        printBefore               ::Bool
        operator                  ::AbstractEeInteraction
        lineSelection             ::LineSelection
    end

    """
    `RadiativeAuger.Settings()`  ... constructor for the default values of RadiativeAuger line computations
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], GreenChannel[], 0, 3, false, CoulombInteraction(), LineSelection())
    end

    # `Base.show(io::IO, settings::RadiativeAuger.Settings)`  ... prepares a proper printout of the variable settings::RadiativeAuger.Settings.
    function Base.show(io::IO, settings::RadiativeAuger.Settings)
        println(io, "multipoles:                   $(settings.multipoles)  ")
        println(io, "gauges:                       $(settings.gauges)  ")
        println(io, "green:                         (settings.green)  ")
        println(io, "NoEnergySharings:             $(settings.NoEnergySharings)  ")
        println(io, "maxKappa:                     $(settings.maxKappa)  ")
        println(io, "printBefore:                  $(settings.printBefore)  ")
        println(io, "operator :                    $(settings.operator )  ")
        println(io, "lineSelection:                $(settings.lineSelection)  ")
    end

    """
    `struct  RadiativeAuger.ReducedChannel`
        ... defines a type for a RadiativeAuger (reduced) channel to help characterize a scattering (continuum) state of many
            electron-states with a single free electron.

        + multipole      ::EmMultipole          ... Multipole of the photon absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + omega          ::Float64              ... photon energy
        + epsilon        ::Float64              ... (free) electron energy
        + kappa          ::Int64                ... partial-wave of the free electron
        + phase          ::Float64              ... phase of the partial wave
        + tSymmetry      ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + amplitude      ::Complex{Float64}     ... Auger amplitude associated with the given channel.
    """
    struct  ReducedChannel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        omega            ::Float64
        epsilon          ::Float64
        kappa            ::Int64
        phase            ::Float64
        tSymmetry        ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Sharing`
        ... defines a type for a RadiativeAuger sharing to help characterize energy sharing between the emitted photon and
            the scattering (continuum) state with a single free electron.

        + omega            ::Float64         ... Energy of the emitted photon.
        + epsilon          ::Float64         ... Energy of the (outgoing free) electron.
        + weight           ::Float64         ... Gauss-Lengendre weight of this sharing for energy-integrated quantities.
        + differentialRate ::EmProperty      ... differential rate of this energy sharing.
        + channels         ::Array{RadiativeAuger.ReducedChannel,1}  ... List of RadiativeAuger channels of this line.
    """
    struct  Sharing
        omega              ::Float64
        epsilon            ::Float64
        weight             ::Float64
        differentialRate   ::EmProperty
        channels           ::Array{RadiativeAuger.ReducedChannel,1}
    end

    # RadiativeAuger line between an initial and final (bound-state) level
    """
    `struct  Line`  ... defines a type for a radiative Auger line that may include the definition of sublines and their corresponding amplitudes.

        + initialLevel   ::Level          ... initial-(state) level
        + finalLevel     ::Level          ... final-(state) level
        + totalRate      ::EmProperty     ... Total rate of this line.
        + sharings       ::Array{RadiativeAuger.Sharing,1}  ... List of RadiativeAuger sharings of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        totalRate        ::EmProperty
        sharings         ::Array{RadiativeAuger.Sharing,1}
    end

    # `Base.show(io::IO, line::RadiativeAuger.Line)`  ... prepares a proper printout of the variable line::RadiativeAuger.Line.
     function Base.show(io::IO, line::RadiativeAuger.Line)
        println(io, "initialLevel:           $(line.initialLevel)  ")
        println(io, "finalLevel:             $(line.finalLevel)  ")
        println(io, "totalRate:              $(line.totalRate)  ")
        println(io, "sharings:               $(line.sharings)  ")
    end

    """
    `RadiativeAuger.computeAmplitudesProperties(line::RadiativeAuger.Line, grid::Radial.Grid, settings::RadiativeAuger.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::RadiativeAuger.Line is returned for which the amplitudes
            and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(line::RadiativeAuger.Line, grid::Radial.Grid, settings::RadiativeAuger.Settings)
        newSharings = RadiativeAuger.Sharing[]
        for sharing in line.sharings
            newChannels = RadiativeAuger.ReducedChannel[]
            for ch in sharing.channels
                # Generate a continuum orbital
                phase     = 1.3
                amplitude = 1.0im
                push!( newChannels, ReducedChannel(ch.multipole, ch.gauge, ch.omega, ch.epsilon, ch.kappa, phase, ch.tSymmetry, amplitude) )
            end
            # Calculate the differential rate
            diffRate = EmProperty(-1., -1.)
            push!( newSharings, RadiativeAuger.Sharing( sharing.omega, sharing.epsilon, sharing.weight, diffRate, newChannels) )
        end
        # Calculate the totalRate
        totalRate = EmProperty(-1., -1.)
        line = RadiativeAuger.Line( line.initialLevel, line.finalLevel, totalRate, newSharings)
        return( line )
    end

    """
    `RadiativeAuger.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                 settings::RadiativeAuger.Settings; output=true)`
        ... to compute the radiative Auger transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{RadiativeAuger.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::RadiativeAuger.Settings; output=true)
        println("")
        printstyled("RadiativeAuger.computeLines(): The computation of radiative Auger rates starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = RadiativeAuger.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    RadiativeAuger.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = RadiativeAuger.Line[]
        for  line in lines
            newLine = RadiativeAuger.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        RadiativeAuger.displayTotalRates(stdout, lines, settings)
        RadiativeAuger.displayDifferentialRates(stdout, lines, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   RadiativeAuger.displayTotalRates(iostream, lines, settings)
                           RadiativeAuger.displayDifferentialRates(iostream, lines, settings)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `RadiativeAuger.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::RadiativeAuger.Settings)`
        ... to determine a list of RadiativeAuger.Line's for transitions between levels from the initial- and final-state multiplets,
            and  by taking into account the particular selections and settings for this computation; an Array{RadiativeAuger.Line,1} is
            returned. Apart from the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::RadiativeAuger.Settings)
        lines = RadiativeAuger.Line[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    energy   = iLevel.energy - fLevel.energy
                    channels = RadiativeAuger.determineSharingsAndChannels(fLevel, iLevel, energy, settings)
                    push!( lines, RadiativeAuger.Line(iLevel, fLevel, EmProperty(0., 0.,), channels) )
                end
            end
        end
        return( lines )
    end

    """
    `RadiativeAuger.determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::RadiativeAuger.Settings)`
        ... to determine a list of RadiativeAuger Sharing's and Channel's for a transitions from the initial to final level and by taking into
            account the particular settings of for this computation; an Array{RadiativeAuger.Sharing,1} is returned.
    """
    function determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::RadiativeAuger.Settings)
        sharings = RadiativeAuger.Sharing[];    eSharings = Basics.determineEnergySharings(energy, settings.NoEnergySharings)
        for  es in eSharings
            omega     = es[1];    epsilon = es[2];    weight = es[3]
            channels  = RadiativeAuger.ReducedChannel[];
            symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
            for  mp in settings.multipoles
                for  gauge in settings.gauges
                    symList = AngularMomentum.allowedMultipoleSymmetries(symi, mp)
                    for  symt in symList
                        kappaList = AngularMomentum.allowedKappaSymmetries(symt, symf)
                        for  kappa in kappaList
                            # Include further restrictions if appropriate
                            if     string(mp)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                                push!(channels, ReducedChannel(mp, Basics.Coulomb,   omega, epsilon, kappa, 0., symt, Complex(0.)) )
                            elseif string(mp)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                                push!(channels, ReducedChannel(mp, Basics.Babushkin, omega, epsilon, kappa, 0., symt, Complex(0.)) )
                            elseif string(mp)[1] == 'M'
                                push!(channels, ReducedChannel(mp, Basics.Magnetic,  omega, epsilon, kappa, 0., symt, Complex(0.)) )
                            end
                        end
                    end
                end
            end
            push!(sharings, RadiativeAuger.Sharing(omega, epsilon, weight, EmProperty(0., 0.), channels) )
        end
        return( sharings )
    end

    """
    `RadiativeAuger.displayDifferentialRates(stream::IO, lines::Array{RadiativeAuger.Line,1}, settings::RadiativeAuger.Settings)`
        ... to display all differential rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayDifferentialRates(stream::IO, lines::Array{RadiativeAuger.Line,1}, settings::RadiativeAuger.Settings)
        #
        # First, print lines and sharings
        nx = 130
        println(stream, " ")
        println(stream, "  Energy-differential rates of selected radiative-Auger lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(38, "Energies (all in " * TableStrings.inUnits("energy") * ")"; na=4);
        sb = sb * TableStrings.flushleft(38, "  i -- f          omega        epsilon"; na=4)
        sa = sa * TableStrings.center(14, "Weight"; na=0);                            sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(34, "Cou -- diff. rate -- Bab"; na=3)
        sb = sb * TableStrings.center(34, TableStrings.inUnits("rate") * "        " *
                                          TableStrings.inUnits("rate"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega))                    * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon))                  * "    "
                sb = sb * @sprintf("%.4e",                                              sharing.weight)                    * "       "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Coulomb))   * "   "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Babushkin)) * "   "
                println(stream,  sb )
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `RadiativeAuger.displayLines(lines::Array{RadiativeAuger.Line,1})`
        ... to display a list of lines, sharings and channels that have been selected due to the prior settings. A neat table
            of all selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{RadiativeAuger.Line,1})
        nx = 170
        println(" ")
        println("  Selected radiative-Auger lines:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(34, "Energy  " * TableStrings.inUnits("energy"); na=5);
        sb = sb * TableStrings.center(34, "i -- f        omega      e_Auger "; na=5)
        sa = sa * TableStrings.flushleft(57, "List of multipoles, gauges, kappas and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(57, "partial (multipole, gauge, total J^P)                  "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "  "
            #
            for  sharing  in  line.sharings
                sb =      @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega))   * "  "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon)) * "    "
                kappaMultipoleSymmetryList = Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry}[]
                for  i in 1:length(sharing.channels)
                    push!( kappaMultipoleSymmetryList, (sharing.channels[i].kappa, sharing.channels[i].multipole, sharing.channels[i].gauge,
                                                        sharing.channels[i].tSymmetry) )
                end
                wa = TableStrings.kappaMultipoleSymmetryTupels(85, kappaMultipoleSymmetryList)
                sc = sa * sb * wa[1];    println( sc )
                for  i = 2:length(wa)
                    sc = TableStrings.hBlank( length(sa*sb) ) * wa[i];    println( sc )
                end
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `RadiativeAuger.displayTotalRates(stream::IO, lines::Array{RadiativeAuger.Line,1}, settings::RadiativeAuger.Settings)`
        ... to display all total rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayTotalRates(stream::IO, lines::Array{RadiativeAuger.Line,1}, settings::RadiativeAuger.Settings)
        #
        # First, print lines and sharings
        nx = 88
        println(stream, " ")
        println(stream, "  Total rates of selected radiative-Auger lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=2);                         sb = sb * TableStrings.hBlank(21)
        sa = sa * TableStrings.center(12, "i--Energy--f"; na=3)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(30, "Cou --   total rate   -- Bab"; na=4)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("rate") * "          " *
                                          TableStrings.inUnits("rate"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))                 * "      "
            #
            sb = sa * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Coulomb))   * "     "
            sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Babushkin)) * "   "
            println(stream,  sb )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
module WarpPlots

using TimeWarp
using Reexport
@reexport using Plots
export dtwplot
pyplot(); # preferred backend

"""
    dtwplot(seq1, seq2, [dist=SqEuclidean()])
    dtwplot(seq1, seq2, i1, i2)

Given two sequences, perform dynamic time warping and plot
the results. If alignment has already been computed, pass
the indices `i1` and `i2` to make the plot.
"""
function dtwplot(
        seq1::TimeWarp.Sequence,
        seq2::TimeWarp.Sequence,
        dist::SemiMetric=SqEuclidean()
    )
    D = dtw_cost_matrix(seq1, seq2, dist)
    cost,i1,i2 = TimeWarp.trackback(D)
    dtwplot(seq1, seq2, D, i1, i2)
end

@userplot DTWPlot

@recipe function f(h::DTWPlot)
    seq1, seq2, D, i1, i2 = h.args
    n1, n2 = length(seq1), length(seq2)

    # set up the subplots
    legend --> false
    link := :both
    grid --> false
    layout --> @layout [
        a  b{0.9w,0.9h}
        _  c
    ]

    left_margin --> 0mm
    bottom_margin --> 0mm
    top_margin --> 0mm
    right_margin --> 0mm
    clim --> (0,3*D[end,end])
    colorbar --> :right

    # heatmap
    @series begin
        seriestype := :heatmap
        formatter --> (z)->""
        c --> :viridis
        subplot := 2
        D
    end

    # the rest of the plots are paths
    seriestype := :path
    linecolor --> RGB(0,0,0)

    # main plot
    @series begin
        linewidth --> 3
        subplot := 2
        formatter --> (z)->""
        i1, i2
    end

    # left line plot
    @series begin
        subplot := 1
        seq2, 1:n2
    end

    # bottom line plot
    @series begin
        subplot := 3
        1:n1, seq1
    end
end

end


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs

@testset "trace" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    # TODO: change to macro
    mode1 = Mode(0.2, 1.0)
    mode2 = Mode(0.3, 2.0)
    Energy = [0.0, 200.0]
    mol1 = Molecule([mode1], 3, [2.0, 200.0])
    mol2 = Molecule([mode2], 3, [3.0, 300.0])
    aggCore = AggregateCore([mol1, mol2])
    aggCore.coupling[2, 3] = 50
    aggCore.coupling[3, 2] = 50
    agg = setupAggregate(aggCore)
    aggTools = agg.tools

    # not dependent on aggregate
    Ham_I = Matrix(rand(StableRNG(1), Float64, agg.tools.bSize, agg.tools.bSize))
    Ham_I = DenseOperator(agg.tools.basis, agg.tools.basis, Ham_I)
    Ham_0 = Matrix(rand(StableRNG(2), Float64, agg.tools.bSize, agg.tools.bSize))
    Ham_0 = DenseOperator(agg.tools.basis, agg.tools.basis, Ham_0)
    Ham = Ham_I + Ham_0

    basis = agg.tools.basis
    indicesLen = agg.tools.bSize
    indices = agg.tools.indices
    indicesMap = agg.tools.indicesMap
    FCFact = agg.tools.FCfactors
    FCProd = agg.tools.FCproduct

    data = Matrix(Hermitian(rand(StableRNG(0), ComplexF64, indicesLen, indicesLen)))
    rho0 = DenseOperator(basis, basis, data)
    normalize!(rho0)

    t = 0.0
    U_op = evolutionOperator(Ham, t)
    rho = U_op * rho0 * U_op'
    rho_traced_ref = ComplexF64[0.30545089951602383 + 0.0im 0.14375197756961597 + 0.21658981002618588im 0.09446278210921694 + 0.06885750587885092im; 0.14375197756961597 - 0.21658981002618588im 0.30379617012935617 + 1.5178830414797062e-18im 0.06403549575681955 + 0.0490347030584035im; 0.09446278210921695 - 0.06885750587885092im 0.06403549575681954 - 0.04903470305840349im 0.3601657072334817 - 1.1629184732212025e-18im]
    rho_traced = trace_bath_slow(rho, aggCore, aggTools)
    @test 1e-7 > D(rho_traced.data, rho_traced_ref)
    rho_traced = trace_bath_slow(rho.data, aggCore, aggTools)
    @test 1e-7 > D(rho_traced, rho_traced_ref)

    rho_traced = trace_bath(rho, aggCore, aggTools)
    @test 1e-7 > D(rho_traced.data, rho_traced_ref)
    rho_traced = trace_bath(rho.data, aggCore, aggTools)
    @test 1e-7 > D(rho_traced, rho_traced_ref)

    for a = 1:3, b = 1:3
        rho_traced_ab = trace_bath(rho, a, b, aggTools)
        @test 1e-7 > abs(rho_traced_ref[a, b] - rho_traced_ab)
    end

    rho_bath = get_rho_bath(rho0, aggCore, aggTools)
    rho_traced = trace_bath(rho_bath, aggCore, aggTools)
    rho_traced_ref = [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]
    @test 1e-10 > D(rho_traced_ref, rho_traced.data)

    rho_traced = trace_bath(rho0, aggCore, aggTools)
    rho0_ad = ad(rho_traced, rho_bath, aggCore, aggTools)
    @test 1e-7 > D(rho0, rho0_ad)

    t = 1.0
    U_op = evolutionOperator(Ham, t)
    rho = U_op * rho0 * U_op'
    rho_traced_ref = ComplexF64[0.2595039421087527 + 1.1102230246251565e-16im 0.6626972429778065 - 0.3034395850675509im 0.12521539202671672 - 0.4920188624519769im; 0.6626972429778069 + 0.3034395850675506im 0.21095069278599016 + 2.2730010622775303e-16im 0.29117337783386915 + 1.2392616289083571im; 0.12521539202671667 + 0.492018862451977im 0.29117337783386926 - 1.239261628908357im 1.3575702817649218 + 4.0047474713266894e-17im]
    rho_traced = trace_bath_slow(rho, aggCore, aggTools)
    @test 1e-7 > D(rho_traced.data, rho_traced_ref)
    rho_traced = trace_bath_slow(rho.data, aggCore, aggTools)
    @test 1e-7 > D(rho_traced, rho_traced_ref)

    rho_traced = trace_bath(rho, aggCore, aggTools)
    @test 1e-7 > D(rho_traced.data, rho_traced_ref)
    rho_traced = trace_bath(rho.data, aggCore, aggTools)
    @test 1e-7 > D(rho_traced, rho_traced_ref)

    for a = 1:3, b = 1:3
        rho_traced_ab = trace_bath(rho, a, b, aggTools)
        @test 1e-7 > abs(rho_traced_ref[a, b] - rho_traced_ab)

        rho_ab = take_el_part(rho.data, a, b, indicesMap)
        rho_traced_ab = trace_bath_part(
            rho_ab,
            a,
            b,
            aggTools
        )
        @test 1e-7 > abs(rho_traced_ref[a, b] - rho_traced_ab)
    end

    rho_bath = get_rho_bath(rho0, aggCore, aggTools)
    rho_traced = trace_bath(rho_bath, aggCore, aggTools)
    rho_traced_ref = [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]
    @test 1e-10 > D(rho_traced_ref, rho_traced.data)

    rho_traced = trace_bath(rho0, aggCore, aggTools)
    rho0_ad = ad(rho_traced, rho_bath, aggCore, aggTools)
    @test 1e-7 > D(rho0, rho0_ad)

    rho_bath = get_rho_bath(rho0, aggCore, aggTools)

    t = 0.0
    Ham_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)
    prod = Ham_II_t * Ham_I * rho_bath
    corr_ref = trace_bath(prod, aggCore, aggTools)
    corr = correlation_function(
        t,
        rho_bath,
        Ham_0,
        Ham_I,
        aggCore,
        aggTools
    )
    @test 1e-7 > D(corr, corr_ref)

    t = 1.0
    Ham_II_t = getInteractionHamIPicture(Ham_0, Ham_I, t)
    prod = Ham_II_t * Ham_I * rho_bath
    corr_ref = trace_bath(prod, aggCore, aggTools)
    corr = correlation_function(
        t,
        rho_bath,
        Ham_0,
        Ham_I,
        aggCore,
        aggTools
    )
    @test 1e-7 > D(corr, corr_ref)

end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
# Autogenerated wrapper script for tree_sitter_php_jll for x86_64-w64-mingw32-cxx03
export libtreesitter_php

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Relative path to `libtreesitter_php`
const libtreesitter_php_splitpath = ["bin", "libtreesitter_php.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libtreesitter_php_path = ""

# libtreesitter_php-specific global declaration
# This will be filled out by __init__()
libtreesitter_php_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libtreesitter_php = "libtreesitter_php.dll"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"tree_sitter_php")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global libtreesitter_php_path = normpath(joinpath(artifact_dir, libtreesitter_php_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libtreesitter_php_handle = dlopen(libtreesitter_php_path)
    push!(LIBPATH_list, dirname(libtreesitter_php_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
"""
File: Kmeans.jl
Authors: Reese Pathak, Stephen Boyd

K-means implementation in Julia

Parameters:
   list_of_vectors - The vectors to cluster, represented as an array of N n-vectors
                     (N column vectors of length n).
   k               - The number of clusters (defaults to 2)

Returns:
   representatives - The cluster representatives represented as an array of k
                     column vectors of length n.
   labels          - An array of N cluster assignments in the same order as data,
                     so that the 1st data point's cluster assignment is labels[i].
                     The labels themselves index into the representatives list,
                     i.e. representatives[labels[1]] returns the first data
                     point's representative.
   losses          - An array containing the value of J at each iteration.
"""

module K_means

export k_means

function k_means(list_of_vectors; k=2, epsilon=1e-6)
  losses = Float64[] # array of losses
  # initialize representatives and assignments
  representatives, assignments = initialize_representatives(list_of_vectors, k)
  while has_not_converged(losses, epsilon)
    assignments = partition_data(list_of_vectors, representatives)
    representatives = update_representatives(list_of_vectors, assignments,
                                              k, representatives)
    push!(losses, loss(list_of_vectors, representatives, assignments))
   end
  return representatives, assignments, losses
end

function has_not_converged(losses, epsilon)
  return !(length(losses) >= 2  && abs(losses[end] - losses[end-1]) <= epsilon)
end
# Assigns each data point to the closest representative
function partition_data(data, representatives)
  labels = Int[]
  for i in 1:length(data)
    # calculate the euclidean distance from the ith data point to each representative
    representative_distances = pairwise_distance(data[i], representatives)
    # find the representative closest to the ith data point (argmin)
    distance, representative_index = findmin(representative_distances)
    push!(labels, representative_index)
  end
  return labels
end

# Calculates each representative as the mean of it's data points
function update_representatives(data, labels, k, old_representatives)
  representatives = []
  for i = 1:k
    # get the data points assigned to the representative
    representative_pts = data[labels .== i]
    # if the representative has assigned data points, set it to be their mean
    if length(representative_pts) > 0
      push!(representatives, mean(representative_pts))
    else # otherwise, leave it alone
      push!(representatives, old_representatives[i])
    end
  end
  return representatives
end

# Calculates the loss (J) of the given clustering.
function loss(data, representatives, assignments)
  N = length(data)
  J = 0.0
  for i = 1:N
    J += norm(data[i] - representatives[assignments[i]])^2
  end
  return J*(1/N)
end

function initialize_representatives(data, k)
  N = length(data)
  n = length(data[1])
  labels = rand(1:k,N)
  representatives = update_representatives(data, labels, k, [zeros(n) for i=1:k])
  return representatives, labels
end

# Calculates the Eudclidean distance between vector and each vectors in 'vectors',
# returning an array.
function pairwise_distance(vector, list_of_vectors)
  return Float64[norm(vector - list_of_vectors[i])
                  for i=1:length(list_of_vectors)]
end

end

max_stars_repo_name: reesepathak/mma-functions

max_issues_repo_name: reesepathak/mma-functions

max_forks_repo_name: reesepathak/mma-functions


text:
__precompile__()

module PotentialFlow

using Reexport

include("Properties.jl")
include("Utils.jl")

include("Elements.jl")
include("RigidBodyMotions.jl")
include("Motions.jl")

@reexport using .Elements
@reexport using .Motions
@reexport using .RigidBodyMotions

include("TimeMarching.jl")
@reexport using .TimeMarching

#== Some Built-in Potential Flow Elements ==#

export Vortex, Source, Sheets, Plates, Plate, Freestreams, Freestream,
       Doublets, Doublet, Bodies, Corners, Corner, Wedge

include("elements/Points.jl")
include("elements/Blobs.jl")
include("elements/Sheets.jl")

include("elements/Vortex.jl")
include("elements/Source.jl")

include("elements/Doublets.jl")
include("elements/Freestream.jl")

include("elements/Plates.jl")
include("elements/Bodies.jl")

include("elements/Corners.jl")

@reexport using SchwarzChristoffel

import .Plates: Plate
import .Bodies: ConformalBody

import .Corners: Corner, Wedge

#== Plot Recipes ==#

include("plot_recipes.jl")

end


text:
function SparseGradient(e::E) where E
    indexer = Dict{Int,Int}()
    d = Gradient(e,indexer)
    sparsity = Vector{Int}(undef,length(indexer))
    set_sparsity!(sparsity,indexer)
    return (d,sparsity)
end
SparseJacobian(s::Sink{Field}) = SparseJacobian(inner(s))
function SparseJacobian(f::Field1{G,I}) where {G,I}
    indexer = Dict{Tuple{Int,Int},Int}()
    jac = Jacobian(f,indexer)
    sparsity = Vector{Tuple{Int,Int}}(undef,length(indexer))
    set_sparsity!(sparsity,indexer)
    return (jac,sparsity)
end
SparseJacobian(::FieldNull) = (FIELD_NULL,Tuple{Int,Int}[])

function SparseHessian(e::E,d::G) where {E,G}
    indexer = Dict{Tuple{Int,Int},Int}()
    hess = Hessian(e,d,indexer)
    sparsity = Vector{Tuple{Int,Int}}(undef,length(indexer))
    set_sparsity!(sparsity,indexer)
    return (hess,sparsity)
end
function SparseLagrangianHessian(obj,grad,con,jac)
    indexer = Dict{Tuple{Int,Int},Int}()
    laghess = LagrangianHessian(obj,grad,con,jac,indexer)
    sparsity = Vector{Tuple{Int,Int}}(undef,length(indexer))
    set_sparsity!(sparsity,indexer)
    return (laghess,sparsity)
end

function set_sparsity!(sparsity,indexer)
    for (p,k) in indexer
        sparsity[k] = p
    end
end
function set_indexer!(indexer,index)
    if haskey(indexer,index)
        return indexer[index]
    else
        return indexer[index] = length(indexer)+1
    end
end
function set_indexer!(indexer,index1,index2)
    if haskey(indexer,(index1,index2))
        return indexer[index1,index2]
    else
        return indexer[index1,index2] = length(indexer)+1
    end
end


text:
using Documenter, DiscreteChoice

makedocs()


text:
using PartialWaveFunctions
using Test

@testset "PartialWaveFunctions.jl" begin
    # Clebsch Gordan coefficient
    @test clebschgordan_doublearg(3,  3, 1,  1, 2*2, 2*2) == 1
    @test clebschgordan_doublearg(3,  3, 1, -1, 2*1, 2*1) ≈ sqrt(3)/2
    @test clebschgordan_doublearg(3, -1, 1,  1, 2*1, 2*0) ≈ -sqrt(2)/2

    # Wigner Function
    @test wignerD(3, 2, 1, π, 0, -π) ≈ -sqrt(10)/8 + 0im
    @test wignerd(2, 0, 0, 1/sqrt(2)) ≈ (3/2-1)/2
    @test wignerd(1, 1, 0, 0.3) ≈ -sqrt(1-0.3^2)/sqrt(2)
    @test wignerd_doublearg(1, 1, 1, 0) ≈ 1/sqrt(2)

    #
    @test sum(kronecker(i,j) for i in 1:5, j in 1:5) == 5
end


text:
"""
    density1d(x::AbstractVector, grid::AbstractVector; normalize::Bool=True)
Compute number density on a line. Normalized such that int n dx =1.
"""
function density1d(x::AbstractVector, grid::AbstractVector; normalize::Bool=true)
    bins = length(grid) -1
    # initialize array to store number
    n = zeros(bins)
    # in a cell [a,b], particles that has a<=x<b is counted.
    # for the last cell, x==b is also counted.
    for i in 1:bins
        for position in x
            if grid[i] <=position <grid[i+1]
                n[i] += 1
            end
        end
    end
    # last cell: add particles on the right boundary
    n[end] += sum(x .== grid[end])
    # use bin centers as representative locations.
    dx = diff(grid)
    bin_centers = grid[1:end-1] + dx/2
    if normalize
        n ./= dx *length(x)
    end
    return bin_centers, n
end

"""
    density1d(x::AbstractVector, xmin::Real, xmax::Real; bins::Int=100)
Compute number density on a line. Normalized such that int n dx =1.
"""
function density1d(x::AbstractVector, xmin::Real, xmax::Real; bins::Int=100, normalize::Bool=true)
    xmin < xmax ||throw(ArgumentError("xmin>=xmax"))
    grid = LinRange(xmin, xmax, bins+1)
    return density1d(x, grid, normalize=normalize)
end

"""
    density1d(x::AbstractVector; bins::Int=100, normalize::Bool=true)
Compute number density on a line.
"""
function density1d(x::AbstractVector; bins::Int=100, normalize::Bool=true)
    xmin = minimum(x)
    xmax = maximum(x)
    return density1d(x, xmin, xmax, bins=bins, normalize=normalize)
end


text:
"""
Partial out variables in a Dataframe

### Arguments
* `df::AbstractDataFrame`
* `model::Model`: A `Model` created using `@model`. See `@model`.
* `add_mean::Bool`: Should the initial mean added to the returned variable?
* `method::Symbol`: A symbol for the method. Default is :lsmr (akin to conjugate gradient descent). Other choices are :qr and :cholesky (factorization methods)
* `maxiter::Integer`: Maximum number of iterations
* `tol::Real`: Tolerance

### Returns
* `::DataFrame`: a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.
* `::Vector{Int}`: a vector of iterations for each column
* `::Vector{Bool}`: a vector of success for each column

### Details
`partial_out` returns the residuals of a set of variables after regressing them on a set of regressors. The syntax is similar to `reg` - but it accepts multiple dependent variables. It returns a dataframe with as many columns as there are dependent variables and as many rows as the original dataframe.
The regression model is estimated on only the rows where *none* of the dependent variables is missing. With the option `add_mean = true`, the mean of the initial variable is added to the residuals.

### Examples
```julia
using  RDatasets, DataFrames, FixedEffectModels, Gadfly
df = dataset("datasets", "iris")
df[:SpeciesC] =  categorical(df[:Species])
result = partial_out(df, @model(SepalWidth + SepalLength ~ 1, fe = SpeciesC), add_mean = true)
plot(layer(result[1], x="SepalWidth", y="SepalLength", Stat.binmean(n=10), Geom.point),
   layer(result[1], x="SepalWidth", y="SepalLength", Geom.smooth(method=:lm)))
```
"""
function partial_out(df::AbstractDataFrame, m::Model; kwargs...)
    partial_out(df, m.f; m.dict..., kwargs...)
end

function partial_out(df::AbstractDataFrame, f::Formula;
    fe::Union{Symbol, Expr, Nothing} = nothing,
    weights::Union{Symbol, Expr, Nothing} = nothing,
    add_mean = false,
    maxiter::Integer = 10000, tol::Real = 1e-8,
    method::Symbol = :lsmr)
    feformula = fe
    weightvar = weights

    rf = deepcopy(f)
    (has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose_iv!(rf)
    has_absorb = feformula != nothing
    if has_iv
        error("partial_out does not support instrumental variables")
    end
    rt = Terms(rf)
    has_weights = (weightvar != nothing)
    xf = @eval(@formula($nothing ~ $(rf.rhs)))
    xt = Terms(xf)

    # create a dataframe without missing values & negative weights
    vars = allvars(rf)
    absorb_vars = allvars(feformula)
    all_vars = vcat(vars, absorb_vars)
    all_vars = unique(convert(Vector{Symbol}, all_vars))
    esample = completecases(df[all_vars])
    if has_weights
        esample .&= isnaorneg(df[weightvar])
    end

    # initialize iterations & converged
    iterations = Int[]
    convergeds = Bool[]

    # Build fixedeffects, an array of AbtractFixedEffects
    if has_absorb
        fes, ids = parse_fixedeffect(df, Terms(@eval(@formula(nothing ~ $(feformula)))))
    end
    nobs = sum(esample)
    (nobs > 0) || error("sample is empty")
    # Compute weight vector
    sqrtw = get_weights(df, esample, weightvar)
    if has_absorb
        # in case there is any intercept fe, remove the intercept
        if any([isa(fe.interaction, Ones) for fe in fes])
            xt.intercept = false
        end
        fes = FixedEffect[_subset(fe, esample) for fe in fes]
        pfe = FixedEffectMatrix(fes, sqrtw, Val{method})
    else
        pfe = nothing
    end

    # Compute residualized Y
    yf = @eval(@formula($nothing ~ $(rf.lhs)))
    yt = Terms(yf)
    yt.intercept = false
    mfY = ModelFrame2(yt, df, esample)
    Y = ModelMatrix(mfY).m
    Y .= Y .* sqrtw
    if add_mean
        m = mean(Y, dims = 1)
    end
    if has_absorb
        Y, b, c = solve_residuals!(Y, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)
    end

    # Compute residualized X
    xvars = allvars(xf)
    if length(xvars) > 0 || xt.intercept
        if length(xvars) > 0
            mf = ModelFrame2(xt, df, esample)
            X = ModelMatrix(mf).m
        else
            X = fill(one(Float64), (length(esample), 1))
        end
        X .= X .* sqrtw
        if has_absorb
            X, b, c = solve_residuals!(X, pfe; maxiter = maxiter, tol = tol)
            append!(iterations, b)
            append!(convergeds, c)
        end
    end

    # Compute residuals
    if length(xvars) > 0 || xt.intercept
        residuals = Y .- X * (X \ Y)
    else
        residuals = Y
    end

    # rescale residuals
    if add_mean
        residuals .= residuals .+ m
    end
    residuals .= residuals ./ sqrtw

    # Return a dataframe
    yvars = convert(Vector{Symbol}, Symbol.(yt.eterms))
    out = DataFrame()
    j = 0
    for y in yvars
        j += 1
        out[y] = Vector{Union{Float64, Missing}}(missing, size(df, 1))
        out[esample, y] = residuals[:, j]
    end
    return out, iterations, convergeds
end


text:
# using CSV, DataFrames

# struct WShape

#     d::Float64
#     tw::Float64
#     bf::Float64
#     tf::Float64
#     kdes::Float64
#     k1::Float64
#     A::Float64
#     Ix::Float64
#     Iy::Float64
#     J::Float64
#     Cw::Float64
#     Zx::Float64
#     Zy::Float64
#     Wno::Float64

# end

# function AISC(shape_name)

#     filename = string(@__DIR__, "/aisc-shapes-database-v15.0.csv")

#     data = CSV.read(filename, DataFrame, header=true)

#     shape_row = findfirst(==(shape_name), data.AISC_Manual_Label)

#     section_type = data[shape_row, :Type]

#     if section_type == "W"

#         d = parse(Float64, data[shape_row, :d])
#         tw = parse(Float64, data[shape_row, :tw])
#         bf = parse(Float64, data[shape_row, :bf])
#         tf = parse(Float64, data[shape_row, :tf])
#         kdes = parse(Float64, data[shape_row, :kdes])

#         #get k1 from AISC table, it is in fraction format
#         k1 = data[shape_row, :k1]
#         index = findfirst("/", k1)

#         whole = Int(data[shape_row, :k1][1] - '0')

#         if isempty(index) == false
#             top_fraction = parse(Float64, data[shape_row, :k1][index[1]-2:index[1]-1])
#             bottom_fraction = parse(Float64, data[shape_row, :k1][index[1]+1:index[1]+2])
#         else
#             top_fraction = 0.0
#             bottom_fraction = 0.0
#         end

#         k1 = whole + top_fraction/bottom_fraction

#         A = data[shape_row, :A]
#         Ix = data[shape_row, :Ix]
#         Iy = data[shape_row, :Iy]
#         J = parse(Float64, data[shape_row, :J])
#         Cw = parse(Float64, data[shape_row, :Cw])
#         Zx = data[shape_row, :Zx]
#         Zy = data[shape_row, :Zy]
#         Wno = parse(Float64, data[shape_row, :Wno])

#         section = WShape(d, tw, bf, tf, kdes,k1, A, Ix, Iy, J, Cw, Zx, Zy, Wno)

#     end

#     return section

# end

using StructuresKit

shape_name = "W14X90"

section = CrossSection.AISC(shape_name)


text:
using Test
using Unitful
@testset "Test kinetic functions" begin
#Arrhenius testing
@testset "Test Arrhenius function" begin
arr = Arrhenius(A=1e6,n=0.5,Ea=41.84*1000.0)
Tlist = 200.0:200.0:2000.0
@test [arr(T=T) for T in Tlist] ≈ [1.6721e-4, 6.8770e1, 5.5803e3, 5.2448e4, 2.0632e5, 5.2285e5, 1.0281e6, 1.7225e6, 2.5912e6, 3.6123e6] rtol=1e-3
end;

#PDepArrhenius testing
@testset "Test pressure dependent Arrhenius function" begin
arr1 = Arrhenius(A=1e6,n=1.0,Ea=10.0*1000.0)
arr2 = Arrhenius(A=1e12,n=1.0,Ea=20.0*1000.0)
Ps = [0.1,10.0].*1e5
parr = PdepArrhenius(Ps=Ps,arrs=[arr1,arr2])
Tlist = 300.0:100.0:1500.0
Plist = [1e4,1e6,1e5]

parr.arrs[1](T=800.0)
parr(T=800.0,P=1.0e5)

for q = 1:3
    P = Plist[q]
    if q == 1
        for T in Tlist
            @test parr(T=T,P=P) ≈ parr.arrs[1](T=T) rtol=1e-5
        end
    elseif q == 2
        for T in Tlist
            @test parr(T=T,P=P) ≈ parr.arrs[2](T=T) rtol=1e-5
        end
    else
        for T in Tlist
            @test parr(T=T,P=P) ≈ sqrt(parr.arrs[1](T=T)*parr.arrs[2](T=T))  rtol=1e-5
        end
    end
end

#MultiArrhenius testing
marr = MultiArrhenius(arrs=[arr1,arr2])
@test marr(T=900.0) ≈ arr1(T=900.0)+arr2(T=900.0) rtol=1e-6
end;

#MultiPdepArrhenius testing
@testset "Test multiple pressure dependent Arrhenius function" begin
arr1 = Arrhenius(A=upreferred((9.3e-16*Na)u"cm^3/(mol*s)").val,n=0.0,Ea=upreferred((4740.0*R*0.001)u"kJ/mol").val)
arr2 = Arrhenius(A=upreferred((9.3e-14*Na)u"cm^3/(mol*s)").val,n=0.0,Ea=upreferred((4740.0*R*0.001)u"kJ/mol").val)
arr3 = Arrhenius(A=upreferred((1.4e-11*Na)u"cm^3/(mol*s)").val,n=0.0,Ea=upreferred((11200.0*R*0.001)u"kJ/mol").val)
arr4 = Arrhenius(A=upreferred((1.4e-9*Na)u"cm^3/(mol*s)").val,n=0.0,Ea=upreferred((11200.0*R*0.001)u"kJ/mol").val)
Ps = [0.1,10.0].*1e5
parr1 = PdepArrhenius(Ps=Ps,arrs=[arr1,arr2])
parr2 = PdepArrhenius(Ps=Ps,arrs=[arr3,arr4])
mparr = MultiPdepArrhenius(parrs=[parr1,parr2])

Tlist = [200.0,400.0,600.0,800.0,1000.0,1200.0,1400.0,1600.0,1800.0,2000.0]
Plist = [1e4,1e5,1e6]

kexplist = [
            [2.85400e-08 4.00384e-03 2.73563e-01 8.50699e+00 1.20181e+02 7.56312e+02 2.84724e+03 7.71702e+03 1.67743e+04 3.12290e+04];
            [2.85400e-07 4.00384e-02 2.73563e+00 8.50699e+01 1.20181e+03 7.56312e+03 2.84724e+04 7.71702e+04 1.67743e+05 3.12290e+05];
            [2.85400e-06 4.00384e-01 2.73563e+01 8.50699e+02 1.20181e+04 7.56312e+04 2.84724e+05 7.71702e+05 1.67743e+06 3.12290e+06];
        ]
for i = 1:length(Tlist)
    for j = 1:length(Plist)
        kexp = kexplist[j,i]
        kact = mparr(T=Tlist[i],P=Plist[j])
        @test kact ≈ kexp rtol=1e-3
    end
end

#ThirdBody testing
tb = ThirdBody(arr=arr1,efficiencies=Dict(5=>2.3))
@test tb(T=900.0,C=5.0) ≈ tb.arr(T=900.0)*5.0 rtol=1e-6
end;

#Lindemann testing
@testset "Test Lindemman function" begin
arrhigh = Arrhenius(A=upreferred(1.39e16u"cm^3/(mol*s)").val,n=-0.534,Ea=upreferred(2.243u"kJ/mol").val)
arrlow = Arrhenius(A=upreferred(2.62e33u"cm^6/(mol^2*s)").val,n=-4.76,Ea=upreferred(10.21u"kJ/mol").val)
efficiencies = Dict(2=>2.3)
lnd = Lindemann(arrhigh=arrhigh,arrlow=arrlow,efficiencies=efficiencies)
Tlist = [300.0,500.0,1000.0,1500.0]
Plist = [1e4,1e5,1e6]
Kexp = [
            [1.38023e+08 2.45661e+08 2.66439e+08];
            [6.09146e+07 2.12349e+08 2.82604e+08];
            [4.75671e+06 4.09594e+07 1.71441e+08];
            [7.03616e+05 6.85062e+06 5.42111e+07];
        ]
for i in 1:length(Tlist)
    for j in 1:length(Plist)
        Kact = lnd(T=Tlist[i],C=Plist[j]/(R*Tlist[i]))
        @test Kact ≈ Kexp[i,j] rtol=1e-3
    end
end
end;

#Troe testing
@testset "Test Troe function" begin
arrhigh = Arrhenius(A=upreferred(1.39e16u"cm^3/(mol*s)").val,n=-0.534,Ea=upreferred(2.243u"kJ/mol").val)
arrlow = Arrhenius(A=upreferred(2.62e33u"cm^6/(mol^2*s)").val,n=-4.76,Ea=upreferred(10.21u"kJ/mol").val)
efficiencies = Dict(2=>2.3)
troe = Troe(arrhigh=arrhigh,arrlow=arrlow,a=0.783,T3=74.0,T1=2941.0,T2=6964.0,efficiencies=efficiencies)
Kexp = [
            [1.00866e+08 2.03759e+08 2.55190e+08];
            [4.74623e+07 1.41629e+08 2.47597e+08];
            [3.97397e+06 2.89521e+07 9.57569e+07];
            [5.91277e+05 5.14013e+06 3.12239e+07];
        ]
Tlist = [300.0,500.0,1000.0,1500.0]
Plist = [1e4,1e5,1e6]
for i in 1:length(Tlist)
    for j in 1:length(Plist)
        Kact = troe(T=Tlist[i],C=Plist[j]/(R*Tlist[i]))
        @test Kact ≈ Kexp[i,j] rtol=1e-3
    end
end
end;

#Chebyshev testing
@testset "Test Chebyshev function" begin
coefs = [[ 5.67723e+00  7.29281e-01 -1.19840e-01  8.82175e-03];
       [-1.02669e+00  8.53639e-01 -3.23485e-02 -2.73670e-02];
       [-4.47011e-01  2.44144e-01  5.59122e-02 -1.01723e-02];
       [-1.28261e-01  1.11596e-02  2.81176e-02  6.04353e-03];
       [-1.17034e-02 -2.35646e-02  6.10090e-04  4.01309e-03];
       [ 1.55433e-02 -1.36846e-02 -4.63048e-03 -2.61353e-04]]

ch = Chebyshev(coefs=coefs,Tmin=300.0,Tmax=2000.0,Pmin=0.01e5,Pmax=100.0e5)
Tlist = [300.0,500.0,1000.0,1500.0]
Plist = [1e4,1e5,1e6]
Kexp2 = [[2291002.84234755 1101983.59484935 43791.93027973 5201.44365151];
       [2584518.40955842 2040371.05643461 236480.90466183 41012.26130665];
       [2572037.53399365 2574282.31472186 857727.08665108 250401.25475447];]

for i = 1:length(Tlist)
    for j = 1:length(Plist)
        Kact = ch(T=Tlist[i],P=Plist[j])
        @test Kact ≈ Kexp2[j,i] rtol=1e-3
    end
end
end;

end;


text:
################################################################################
# VehicleState
################################################################################

mutable struct VehicleState{T}
    x::T # x-position
    y::T # y-position
    θ::T # heading
    v::T # velocity
end

function VehicleState()
    x = 0.0
    y = 0.0
    θ = 0.0
    v = 0.0
    return VehicleState{typeof(x)}(x,y,θ,v)
end

import Base.==
import Base.+
import Base.*
import Base.copy

function (==)(x1::VehicleState{T}, x2::VehicleState{T}) where {T}
	out = true
	for name in fieldnames(VehicleState)
		out &= getfield(x1, name) == getfield(x2, name)
	end
	return out
end

function (+)(x1::VehicleState{T}, x2::VehicleState{T}) where {T}
	return VehicleState{T}(x1.x+x2.x, x1.y+x2.y, x1.θ+x2.θ, x1.v+x2.v)
end

function (*)(x1::VehicleState{T}, α::T) where {T}
	return VehicleState{T}(x1.x*α, x1.y*α, x1.θ*α, x1.v*α)
end

function (*)(α::T, x1::VehicleState{T}) where {T}
	return VehicleState{T}(x1.x*α, x1.y*α, x1.θ*α, x1.v*α)
end

function copy(x::VehicleState{T}) where {T}
	return VehicleState{T}(x.x, x.y, x.θ, x.v)
end

################################################################################
# standardize
################################################################################

function standardize(model::DoubleIntegratorGame{n,m,p}, x::SVector{n,T}) where {n,m,p,T}
    return [standardize(model, x[model.pz[i]]) for i=1:p]
end

function standardize(model::UnicycleGame{n,m,p}, x::SVector{n,T}) where {n,m,p,T}
    return [standardize(model, x[model.pz[i]]) for i=1:p]
end

function standardize(model::BicycleGame{n,m,p}, x::SVector{n,T}) where {n,m,p,T}
    return [standardize(model, x[model.pz[i]]) for i=1:p]
end

function standardize(model::DoubleIntegratorGame{n,m,p}, x::SVector{ni,T}) where {n,m,p,ni,T}
    vs = VehicleState()
    vs.x = x[1]
    vs.y = x[2]
    vs.θ = atan(x[4],x[3])
    vs.v = sqrt(x[3]^2+x[4]^2)
    return vs
end

function standardize(model::UnicycleGame{n,m,p}, x::SVector{ni,T}) where {n,m,p,ni,T}
    vs = VehicleState()
    vs.x = x[1]
    vs.y = x[2]
    vs.θ = x[3]
    vs.v = x[4]
    return vs
end

function standardize(model::BicycleGame{n,m,p}, x::SVector{ni,T}) where {n,m,p,ni,T}
    # x = [x,y,v,ψ]
    vs = VehicleState()
    vs.x = x[1]
    vs.y = x[2]
    vs.θ = x[4] # ψ
    vs.v = x[3]
    return vs
end

################################################################################
# specialize
################################################################################

function specialize(model::AbstractGameModel, vs::AbstractVector{VehicleState})
    p = length(vs)
    return [specialize(model, vs[i]) for i=1:p]
end

function specialize(model::UnicycleGame{n,m,p}, x::SVector{n,T}) where {n,m,p,T}
    return [specialize(model, x[model.pz[i]]) for i=1:p]
end

function specialize(model::BicycleGame{n,m,p}, x::SVector{n,T}) where {n,m,p,T}
    return [specialize(model, x[model.pz[i]]) for i=1:p]
end

function specialize(model::DoubleIntegratorGame{n,m,p}, vs::VehicleState{T}) where {n,m,p,T}
    x = @SVector [vs.x, vs.y, vs.v*cos(vs.θ), vs.v*sin(vs.θ)]
    return x
end

function specialize(model::UnicycleGame{n,m,p}, vs::VehicleState{T}) where {n,m,p,T}
    x = @SVector [vs.x, vs.y, vs.θ, vs.v]
    return x
end

function specialize(model::BicycleGame{n,m,p}, vs::VehicleState{T}) where {n,m,p,T}
    x = @SVector [vs.x, vs.y, vs.v, vs.θ]
    return x
end

@generated function specialize(model::DoubleIntegratorGame{n,m,p}, vs::AbstractVector{VehicleState{T}}) where {n,m,p,T}
	x = [:(vs[$i].x) for i=1:p]
	y = [:(vs[$i].y) for i=1:p]
	vx = [:(vs[$i].v*cos(vs[$i].θ)) for i=1:p]
	vy = [:(vs[$i].v*sin(vs[$i].θ)) for i=1:p]
	return :(SVector{$n}($(x...), $(y...), $(vx...), $(vy...)))
end

@generated function specialize(model::UnicycleGame{n,m,p}, vs::AbstractVector{VehicleState{T}}) where {n,m,p,T}
	x = [:(vs[$i].x) for i=1:p]
	y = [:(vs[$i].y) for i=1:p]
	θ = [:(vs[$i].θ) for i=1:p]
	v = [:(vs[$i].v) for i=1:p]
	return :(SVector{$n}($(x...), $(y...), $(θ...), $(v...)))
end

@generated function specialize(model::BicycleGame{n,m,p}, vs::AbstractVector{VehicleState{T}}) where {n,m,p,T}
	x = [:(vs[$i].x) for i=1:p]
	y = [:(vs[$i].y) for i=1:p]
	v = [:(vs[$i].v) for i=1:p]
	ψ = [:(vs[$i].θ) for i=1:p]
	return :(SVector{$n}($(x...), $(y...), $(v...), $(ψ...)))
end


text:
"""
  Middleware

Http request-handling middleware.
"""
module Middleware

using HTTP: Response, Request

using ..Resources
using ..HttpErrors
using ..HttpErrors: HttpError
using ..HttpMethods

export affix_newline, authorize, handle_errors, log_access

# affix_newline

function _affix_newline(handler, req::Request)
  res = handler(req)
  body = String(res.body)
  if !(endswith(body, "\n"))
    res.body = codeunits(body)
    return res
  end
  return res
end

"""
  affix_newline(handler)

Middleware to affix a newline to the end of a response if it doesn't already have one.
"""
affix_newline(handler) = req -> _affix_newline(handler, req)

# handle_errors

function _handle_errors(handler, req::Request)::Response
  try
    return handler(req)
  catch err
    @error exception=(err, catch_backtrace())
    if isa(err, HttpError)
      msg = sprint(showerror, err)
      return Response(err.status, msg)
    else
      return Response(500, "500 Internal Server Error")
    end
  end
end

"""
  handle_errors(handler)

Wraps an api handler with functionality to automatically handle and translate
application-layer errors into HTTP responses with the appropriate status code.

Returns an api handler, i.e. intended with type
`(Request -> Response) -> (Request -> Response)`.
"""
handle_errors(handler) = req -> _handle_errors(handler, req)

# authorize

function _authorize(handler, req::Request)::Response
  @warn "`authorize` not implemented yet"
  return handler(req)
end

"""
  authorize(handler)

Wraps an api handler with functionality to restrict access to a resource.

Returns an api handler, i.e. intended with type
`(Request -> Response) -> (Request -> Response)`.
"""
authorize(handler) = req -> _authorize(handler, req)

# log_access

function _log_access(handler, req::Request)::Response
  @warn "`log_access` not implemented yet"
  return handler(req)
end

"""
  log_access(handler)

Wraps an api handler with access-logging functionality.

Returns an api handler, i.e. intended with type
`(Request -> Response) -> (Request -> Response)`.
"""
log_access(handler) = req -> _log_access(handler, req)

"""
  Routers

Router middleware.
"""
module Routers

using HTTP: Request
using HTTP.URIs: URI, splitpath

using ..Resources
using ..Resources: Resource, path
using ..HttpErrors: method_not_allowed, not_found
using ..HttpMethods: HttpMethod, http_method

# Router
# Because HTTP.@register sucks

const RoutePair = Pair{Resource, HttpMethod}

struct Router
  base_path::AbstractString
  routes::Vector{RoutePair}
end

# route

function _match(template::AbstractString, path::AbstractString)
  vtemplate = splitpath(template)
  vpath = splitpath(path)
  if length(vtemplate) != length(vpath)
    return false
  end
  return all(
    x -> first(x) == "*" ? true : first(x) == last(x),
    zip(vtemplate, vpath)
  )
end

"""
    route(router::Router, request::HTTP.Request)

Route a `request` through a `router`.
"""
function route(router::Router, req::Request)
  method = http_method(req.method)
  req_path = req.target |> URI |> x -> x.path

  potential_routes = filter(y -> _match(y |> first |> path, req_path), router.routes)
  if isempty(potential_routes)
    not_found()
  end

  idx = findfirst(x -> x |> last == method, potential_routes)
  if isnothing(idx)
    method_not_allowed("$method not allowed for $path")
  end

  @inbounds resource = potential_routes[idx] |> first

  return Resources.handle(method, resource, req)
end

"""
  route(router::Router)

Route requests through a `router`. Throws 404s and 405s when necessary.
"""
route(router::Router) = req -> route(router, req)

end # module Routers

end # module Middleware