text:
"""
    struct Density{M,B}
        μ::M
        base::B
    end

For measures μ and ν with μ≪ν, the density of μ with respect to ν (also called
the Radon-Nikodym derivative dμ/dν) is a function f defined on the support of ν
with the property that for any measurable a ⊂ supp(ν), μ(a) = ∫ₐ f dν.

Because this function is often difficult to express in closed form, there are
many different ways of computing it. We therefore provide a formal
representation to allow comptuational flexibilty.
"""
struct Density{M,B}
    μ::M
    base::B
end

"""
    struct DensityMeasure{F,B} <: AbstractMeasure{X}
        density :: F
        base    :: B
    end

A `DensityMeasure` is a measure defined by a density with respect to some other
"base" measure
"""
struct DensityMeasure{X,F,B} <: AbstractMeasure{X}
    density :: F
    base    :: B
end

# function density(μ::M, ν::M) where {M}
#     if  μ==ν
#         return () -> 1.0
#     end
# end

density(μ, base=basemeasure(μ)) = Density(μ, base)
logdensity(μ, base=basemeasure(μ)) = LogDensity(μ, base)

density(μ::Dists.Distribution{Dists.Univariate,Dists.Continuous}, x::Real) = pdf(μ,x)
logdensity(μ::Dists.Distribution{Dists.Univariate,Dists.Continuous}, x::Real) = logpdf(μ,x)

density(μ::AbstractMeasure{X}, x::X) where {X} = density(μ, baseMeasure(μ))(x)

logdensity(μ::AbstractMeasure{X}, y::Y) where {X, Y <: X} = logdensity(μ, baseMeasure(μ))(x)


text:
using KeywordCalls, BenchmarkTools

letters = Symbol.('a':'z')

for n in 1:26
    fkeys = Tuple(letters[1:n])

    @eval begin
        f(nt::NamedTuple{$fkeys}) = sum(values(nt))
        $(KeywordCalls._kwcall(:(f($(fkeys...)))))
    end
end

function runbenchmark()
    times = Matrix{Float64}(undef, 26,2)
    for n in 1:26
        fkeys = Tuple(letters[1:n])
        rkeys = reverse(fkeys)

        nt = NamedTuple{fkeys}(1:n)
        rnt = NamedTuple{rkeys}(n:-1:1)

        times[n,1] = @belapsed($f($nt))
        times[n,2] = @belapsed($f($rnt))
    end
    return times
end

times = runbenchmark()

using Plots

plt = plot(times[:,1] .* 1e9, legend=:topleft, label="Preferred order", dpi=200)
plot!(plt, times[:,2] .* 1e9, label="Reverse order")
xlabel!(plt, "Number of keys")
ylabel!(plt, "Time (ns)")

ext: jl


text:
# Functionality for reading a mesh from COMSOL
# Currently 2D, only 1 object, only domain elements, currently only triangles

# 21 lines before the coordinates start
# this is not robust
const LINES_BEFORE_CORDS = 22
const LINES_BETWEEN_CORDS_ELEMENT = 10

function read_mphtxt(filename)
    mesh = GeoMesh()
    open(filename, "r") do f
        for i in 1:LINES_BEFORE_CORDS
            readline(f)
        end

        read_coordinates(f, mesh)

        for i in 1:LINES_BETWEEN_CORDS_ELEMENT
            readline(f)
        end

        read_elements(f, mesh)
    end
    return mesh
end

function read_coordinates(f, mesh)
    n_nodes = 1
    while true
        line = strip(readline(f))
        # The list of coordinates is broken by an empty line
        line == "" && break
        coords = split(line, " ")
        coords = [parse(Float64, c) for c in coords]
        push!(mesh, GeoNode2(n_nodes, Point2(coords[1], coords[2])))
        n_nodes += 1
    end
end

function read_elements(f, mesh)
    n_eles = 1
    while true
        line = strip(readline(f))
        # The list of elements is broken by an empty line
        line == "" && break
        vs = split(line, " ")
        vs = [parse(Int, v)+1 for v in vs]
        push!(mesh, GeoTrig(n_eles, Vertex3(vs[1], vs[2], vs[3])))
        n_eles += 1
    end
end


text:
export rec2coeff, showpoly

title_color =   :underline
name_color  =   :light_blue

function show(io::IO,m::Measure)
    printstyled("Measure dλ(t)=w(t)dt\n";color=title_color)
    print("name:\t\t")
    printstyled("$(m.name)\n"; color=name_color)
    println("w(t):\t\t$(m.w)")
    println("dom:\t\t($(m.dom[1]), $(m.dom[2]))")
    print("symmetric:\t")
    m.symmetric ? col=:green : col=:red
    printstyled("$(m.symmetric)\n";color=col)
    if typeof(m.pars)!=Dict{Any,Any}
        print("pars:")
        for (i,d) in enumerate(m.pars)
            printstyled("\t\t$d\n")
        end
    end
end

function show(io::IO,m::MultiMeasure)
    p = length(m.name)
    printstyled("$p-variate measure dλ(t)= ∏_{i=1}^$p w_i(t)dt_i\n";color=title_color)
    print("name:")
    [ printstyled("\t\t$(m.name[i])\n"; color=name_color) for i = 1:p]
    println("w(t):\t\t∏_{i=1}^$p w_i(t)")
    print("dom:")
    [ printstyled("\t\t($(m.dom[i][1]), $(m.dom[i][2]))\n") for i=1:p ]
    print("symmetric:")
    for i=1:p
        m.symmetric[i] ? col=:green : col=:red
        i==1 ? printstyled("\t$(m.symmetric[i])\n";color=col) : printstyled("\t\t$(m.symmetric[i])\n";color=col)
    end
    if length(findall(x->typeof(x)!=Dict{Any,Any},m.pars))>=1 print("pars:")
        for i=1:p
            print("\t\t")
            if typeof(m.pars[i])!=Dict{Any,Any}
                for (i,d) in enumerate(m.pars[i])
                    i==1 ? print("$d") : print(", $d")
                end
                print("\n")
            else
                print("∅\n")
            end
        end
    end
end

function show(io::IO,op::OrthoPoly;showmeasure::Bool=true)
    printstyled("\nUnivariate orthogonal polynomials\n";color=title_color)
    print("name:"); printstyled("\t\t$(op.name)\n";color=name_color)
    print("degree:"); printstyled("\t\t$(op.deg)\n")
    print("#coeffs:"); printstyled("\t$(length(op.α))\n")
    length(op.α)<=7 ? n=length(op.α) : n=7
    print("α ="); printstyled("\t\t$(op.α[1:n])")
    n<length(op.α) ? print("...\n") : print("\n")
    print("β ="); printstyled("\t\t$(op.β[1:n])")
    n<length(op.α) ? print("...\n") : print("\n")
    if showmeasure
        print("\n")
        show(io,op.meas)
    end
end

function show(io::IO,q::Quad;showmeasure::Bool=true)
    printstyled("\nQuadrature rule\n";color=title_color)
    print("name:"); printstyled("\t\t$(q.name)\n";color=name_color)
    print("N:"); printstyled("\t\t$(q.Nquad)\n")
    q.Nquad<=7 ? n=q.Nquad : n=7
    print("nodes"); printstyled("\t\t$(q.nodes[1:n])")
    n<q.Nquad ? print("...\n") : print("\n")
    print("weights"); printstyled("\t\t$(q.weights[1:n])")
    n<q.Nquad ? print("...\n") : print("\n")
    if showmeasure
        print("\n")
        show(io,q.meas)
    end
end

function show(io::IO,opq::OrthoPolyQ;showmeasure::Bool=false)
    show(io,opq.op)
    show(io,opq.quad;showmeasure=showmeasure)
end

function show(io::IO,mop::MultiOrthoPoly;showmeasure::Bool=false)
    p = length(mop.name)
    printstyled("\n$p-variate orthogonal polynomials\n";color=title_color)
    print("name:");
    [ printstyled("\t\t$(mop.name[i])\n";color=name_color) for i=1:p ]
    print("deg:"); printstyled("\t\t$(mop.deg)\n")
    print("dim:"); printstyled("\t\t$(mop.dim)\n")
    print("ind:");
    mop.dim<=7 ? n=mop.dim : n=7
    # printstyled("\t$(mop.ind[1,:])\n")
    [ printstyled("\t\t$(mop.ind[i,:])\n") for i=1:n ]
    n<mop.dim ? print("\t\t...\n") : print("\n")

    print("\n")
    showmeasure ? show(io,mop.meas) : ()
end

function show(io::IO,t::Tensor)
    printstyled("\n$(t.dim)-dimensional tensor\n";color=title_color)
    print("dim:"); printstyled("\t\t$(t.dim)\n")
    print("nonzeros:"); printstyled("\t$(length(t.T.nzind))\n")
    show(io,t.op;showmeasure=false)
end

"""
```
rec2coeff(deg::Int,a::Vector{Float64},b::Vector{Float64})
rec2coeff(a,b) = rec2coeff(length(a),a,b)
```
Get the coefficients of the orthogonal polynomial of degree up to `deg` specified by its
recurrence coefficients `(a,b)`. The function returns the values ``c_i^{(k)}`` from
```math
p_k (t) = t^d + \\sum_{i=0}^{k-1} c_i t^i,
```
where ``k`` runs from `1` to `deg`.

The call `rec2coeff(a,b)` outputs all possible recurrence coefficients given `(a,b)`.
"""
function rec2coeff(deg::Int,a::Vector{Float64},b::Vector{Float64})
    @assert deg >= 1 "degree must be positive (you asked for degree = $deg)"
    @assert length(a) == length(b) && length(a) >= deg "incorrect number of recurrence coefficients"
    c = Vector{Vector{Float64}}(undef,deg)
    @inbounds c[1] = [ -a[1] ]
    deg == 1 && return c
    @inbounds c[2] = [ a[1]*a[2] - b[2], -a[1] - a[2] ]
    deg == 2 && return c
    for k in 2:deg-1
        c[k+1] = [ -a[k+1]*c[k][1] - b[k+1]*c[k-1][1] ]
        for i = 1:k-2
            push!(c[k+1],c[k][i] - a[k+1]*c[k][i+1] - b[k+1]*c[k-1][i+1])
        end
        push!(c[k+1],c[k][k-1] - a[k+1]*c[k][k] - b[k+1])
        push!(c[k+1],c[k][k] - a[k+1])
    end
    return c
end
rec2coeff(α,β) = rec2coeff(length(α),α,β)

"""
```
showpoly(coeffs::Vector{Float64};sym::String,digits::Integer)
showpoly(coeffs::Vector{Vector{Float64}};sym::String,digits::Integer)
```
Show the monic polynomial with coefficients `coeffs` in a human readable way.
They keyword `sym` sets the name of the variable, and `digits` controls the number of shown digits.
```@repl
julia> showpoly([1.2, 2.3, 3.4456])
x^3 + 3.45x^2 + 2.3x + 1.2
julia> showpoly([1.2, 2.3, 3.4456], sym="t", digits=2)
t^3 + 3.45t^2 + 2.3t + 1.2
```

Tailored to types from `PolyChaos.jl`
```
showpoly(d::Integer,α::Vector{Float64},β::Vector{Float64};upto::Bool,sym::String,digits::Integer)
showpoly(d::Integer,op::OrthoPoly;upto::Bool=true,sym::String,digits::Integer)
showpoly(d::Integer,opq::OrthoPolyQ;upto::Bool=true,sym::String,digits::Integer)
```
Show the monic orthogonal polynomials with recurrence coefficients `(α,β)` up to degree `d`.
Setting the keyword `upto` to false prints the monic polynomial of degree equal to `d`.
```@repl
julia> op = OrthoPoly("gaussian",5);
julia> showpoly(3,op)
1
x
x^2 - 1.0
x^3 - 3.0x

julia> showpoly(3,op; upto=false)
x^3 - 3.0x

julia> showpoly(3,op; sym="t")
1
t
t^2 - 1.0
t^3 - 3.0t
```

The following function calls show all orthogonal polynomials given `(α,β)`.
```
showpoly(α::Vector{Float64},β::Vector{Float64};sym::String,digits::Integer)
showpoly(op::Union{OrthoPoly,OrthoPolyQ};sym::String,digits::Integer)
```

```@repl
julia> showpoly(op; sym="y")
1
y
y^2 - 1.0
y^3 - 3.0y
y^4 - 6.0y^2 + 3.0
y^5 - 10.0y^3 + 15.0y
y^6 - 15.0y^4 + 45.0y^2 - 15.0
```

[Thanks @pfitzseb for providing this functionality.](https://discourse.julialang.org/t/how-to-define-verbose-output-for-a-polynomial/21317/5)

"""
function showpoly(coeffs::Vector{Float64};sym::String="x",digits::Integer=2)
    io = Base.stdout
    length(coeffs) > 1 ? print(io, sym*"^", length(coeffs)) : print(io, sym)
    for (i, c) in enumerate(reverse(coeffs))
        abs(round(c, digits=digits)) == 0. && continue
        ex = length(coeffs) - i
        print(io, ' ', c > 0 ? '+' : '-', ' ')
        print(io, abs(round(c, digits=digits)))
        ex > 0 && print(io, sym)
        ex > 1 && print(io, '^', ex)
    end
    print(io, '\n')
end

function showpoly(coeffs::Vector{Vector{Float64}};sym::String="x",digits::Integer=2)
    print("1\n")
    map(c->showpoly(c,sym=sym,digits=digits),coeffs)
    print()
end

function showpoly(d::Integer,α::Vector{Float64},β::Vector{Float64};upto::Bool=true,sym::String="x",digits::Integer=2)
    @assert d >= 0 "degree has to be non-negative."
    d == 0 && return print("1")
    cfs = upto ? rec2coeff(d,α,β) : rec2coeff(d,α,β)[end]
    showpoly(cfs,sym=sym,digits=digits)
end

showpoly(d::Integer,op::OrthoPoly;upto::Bool=true,sym::String="x",digits::Integer=2) = showpoly(d,op.α,op.β;upto=upto,sym=sym,digits=digits)
showpoly(d::Integer,opq::OrthoPolyQ;upto::Bool=true,sym::String="x",digits::Integer=2) = showpoly(d,opq.op;upto=upto,sym=sym,digits=digits)

showpoly(α::Vector{Float64},β::Vector{Float64};sym::String="x",digits::Integer=2) = showpoly(length(α),α,β;sym=sym,digits=digits)
showpoly(op::Union{OrthoPoly,OrthoPolyQ};sym::String="x",digits::Integer=2) = showpoly(size(coeffs(op),1),op;sym=sym,digits=digits)


text:
immutable Node{T}
    isleaf::Bool
    index::Int
    radius::T
    inside::Node{T}
    outside::Node{T}

    Node{T}() where T = new(true, zero(Int))
    Node{T}(k::Int) where T = new(true, k)
    Node{T}(k::Int, d::T, l::Node{T}, r::Node{T}) where T = new(false, k, d, l, r)

    function Node{T}(metric::Function, indices::Vector{Int}, i::Int, j::Int) where T
        j <  i && return Node{T}()
        @inbounds j == i && return Node{T}(indices[i])

        # Choose a random element as the vantage point
        r = rand(i:j)
        @inbounds vp = indices[r]
        @inbounds indices[i], indices[r] = indices[r], indices[i]
        i += 1

        # Find the index of the median
        m = (i + j) >> 1

        d = quickselect!(T, metric, indices, vp, i, j, m)::T

        inside = Node{T}(metric, indices, i, m)
        outside = Node{T}(metric, indices, m + 1, j)
        Node{T}(vp, d, inside, outside)
    end
end

immutable VPTree
    num_points::Int
    root::Node
    metric::Function
end

function partition!(metric::Function, indices::Vector{Int}, vp::Int, d, i::Int, j::Int)
    while true
        @inbounds while metric(vp, indices[i]) <= d
            i == j && return j
            i += 1
        end

        @inbounds while metric(vp, indices[j]) > d
            i == j && return j - 1
            j -= 1
        end

        @inbounds indices[i], indices[j] = indices[j], indices[i]
    end
end

function quickselect!{T}(::Type{T}, metric::Function, indices::Vector{Int}, vp::Int, i::Int, j::Int, k::Int)
    while true
        # Choose a random element as the pivot
        r = rand(i:j)
        @inbounds d = metric(vp, indices[r])::T
        n = partition!(metric, indices, vp, d, i, j)

        if n > k
            j = n
        elseif n < k
            i = n + 1
        else
            return d
        end
    end
end

"""
    VPTree(metric, indices)

Construct a VPTree by providing

1) a function metric(i::Int, j::Int) that computes the distance between point i and point j
2) a list of point indices for which the VP Tree is constructed
"""
function VPTree(metric::Function, indices::Vector{Int})
    n = length(indices)
    T = typeof(metric(first(indices), first(indices)))
    root = Node{T}(metric, indices, 1, n)
    VPTree(n, root, metric)
end

function distsq(points::Matrix, i1::Int, i2::Int)
    x = points[1, i1] - points[1, i2]
    y = points[2, i1] - points[2, i2]
    z = points[3, i1] - points[3, i2]

    x^2 + y^2 + z^2
end


text:
# For optimal performance, run this script with
#   julia --threads=1 --check-bounds=no compare_entropy_projection_variants.jl
# to let Julia use half of your `Sys.CPU_THREADS` cores. You can also specify
# the number of processes to use via
#   julia --threads=1 --procs=N --check-bounds=no compare_entropy_projection_variants.jl
# where `N` is the number of processes you want to use.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

using Distributed

# Add worker processes if none are present
if nprocs() == 1
  addprocs(Sys.CPU_THREADS ÷ 2)
end

# Package setup required on each processor
@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
end
@everywhere using Trixi
@everywhere using OrdinaryDiffEq
@everywhere using Setfield

# used to modify discretization matrices for approximation_type = Val{:staggered}()
@everywhere using Trixi.StartUpDG
@everywhere using LinearAlgebra: diagm, norm, I

# Setup on the driver process
using CSV, DataFrames

@everywhere begin
    function set_initial_condition(atwood_number)
        A = atwood_number
        rho1 = 0.5 * one(A) # recover original with A = 3/7
        rho2 = rho1 * (1 + A) / (1 - A)
        function initial_condition_khi(x, t, equations::CompressibleEulerEquations2D)
            # domain size is [-1,+1]^2

            # B is a discontinuous function with value 1 for -.5 <= x <= .5 and 0 elsewhere
            slope = 15
            B = 0.5 * (tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5))

            rho = rho1 + B * (rho2 - rho1)  # rho ∈ [rho_1, rho_2]
            v1 = B - 0.5                    # v1  ∈ [-.5, .5]
            v2 = 0.1 * sin(2 * pi * x[1])
            p = 1.0
            return prim2cons(SVector(rho, v1, v2, p), equations)
        end
        return initial_condition_khi
    end

    function build_solver(approximation_type::Val{:dgsem}, polydeg, cells_per_dimension, initial_condition, equations)
        dg = DGMulti(polydeg = polydeg, element_type = Quad(), approximation_type = SBP(),
                    surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
        mesh = DGMultiMesh(dg, cells_per_dimension=cells_per_dimension, periodicity=true)
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
        return semi
    end

    function build_solver(approximation_type::Val{:gauss}, polydeg, cells_per_dimension, initial_condition, equations)
        dg = DGMulti(polydeg = polydeg, element_type = Quad(), approximation_type = GaussSBP(),
                    surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
        mesh = DGMultiMesh(dg, cells_per_dimension=cells_per_dimension, periodicity=true)
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
        return semi
    end

    # hybrid DGSEM (Lobatto on volume, Gauss on face)
    function build_solver(approximation_type::Val{:hybrid}, polydeg, cells_per_dimension, initial_condition, equations)
        N = polydeg
        r1D, w1D = StartUpDG.gauss_lobatto_quad(0, 0, N)
        rq, sq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D))
        wq = .*(vec.(StartUpDG.NodesAndModes.meshgrid(w1D))...)

        # Clenshaw-Curtis quadrature
        function clenshaw_curtis(polydeg)
            n = polydeg + 1 # number of Chebyshev nodes
            r = sort(@. cos((2*(1:n)-1)/(2*n)*pi)) # first kind
            Vq_approx = StartUpDG.vandermonde(StartUpDG.Line(), length(r)-1, r)
            w = Vq_approx' \ [sqrt(2); zeros(size(Vq_approx, 2)-1)]
            return r, w
        end
        r1D, w1D = clenshaw_curtis(N)

        dg = DGMulti(polydeg = N, element_type = Quad(), approximation_type = Polynomial(),
                    surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
                    quad_rule_vol = (rq, sq, wq), quad_rule_face=(r1D, w1D))
        mesh = DGMultiMesh(dg, cells_per_dimension=cells_per_dimension, periodicity=true)
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)
    end

    function build_solver(approximation_type::Val{:staggered}, polydeg, cells_per_dimension, initial_condition, equations)
        N = polydeg

        # (N+1)-Lobatto/Gauss staggering. Uses entropy projection at Gauss nodes
        # but evaluates fluxes at degree (N+1) Lobatto volume and face nodes
        r1D, w1D = StartUpDG.gauss_lobatto_quad(0, 0, N)
        rq, sq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D))
        wq = .*(vec.(StartUpDG.NodesAndModes.meshgrid(w1D))...)
        dg = DGMulti(polydeg = N, element_type = Quad(), approximation_type = Polynomial(),
                    surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
                    volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
                    quad_rule_vol = (rq, sq, wq), quad_rule_face=(r1D, w1D))

        mesh = DGMultiMesh(dg, cells_per_dimension=cells_per_dimension, periodicity=true)
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg)

        # replace Lobatto volume quadrature with Gauss nodes
        r1D, w1D = StartUpDG.gauss_quad(0, 0, N)
        rq, sq = vec.(StartUpDG.NodesAndModes.meshgrid(r1D))
        wq = .*(vec.(StartUpDG.NodesAndModes.meshgrid(w1D))...)

        # keep face quadrature = Lobatto nodes
        r1D, w1D = StartUpDG.gauss_lobatto_quad(0, 0, N)

        # build lobatto RefElemData
        rd = RefElemData(Quad(), N, quad_rule_vol=(rq,sq,wq), quad_rule_face=(r1D,w1D))
        Vh = [I(rd.Np); rd.Vf] # map from nodes (Lobatto) to interpolation points (Lobatto)

        VhP = Vh * rd.Pq # rebuild VhP - rd.Pq maps from Gauss to Lobatto
        Ph = rd.M \ Matrix(Vh') # rebuild Ph - rd.M is now dense instead of full

        _, VhP, Ph = hybridized_SBP_operators(rd)

        # swap out quadrature interpolation, projection, lift operators
        # (anything involving the volume quadrature nodes or mass matrix).
        # These operators are not used in the actual flux differencing kernel
        dg.basis.Vq .= rd.Vq # maps from Lobatto nodes to Gauss nodes
        semi.cache.VhP .= VhP
        semi.cache.Ph .= Ph

        # new LIFT built from Gauss volume quadrature and Lobatto surface quadrature
        dg.basis.LIFT .= rd.LIFT

        # just used for AnalysisCallback; not strictly necessary
        mesh.md.wJq .= diagm(rd.wq) * (rd.Vq * mesh.md.J)

        return semi
    end

    function blow_up_time(approximation_type, polydeg, cells_per_dimension, Atwood_number;
                          abstol=1.0e-7, reltol=1.0e-7)

    	println("Running with approximation_type=$(approximation_type), polydeg=$polydeg, cells_per_dimension=$cells_per_dimension, Atwood=$Atwood_number")

        initial_condition = set_initial_condition(Atwood_number)
        equations = CompressibleEulerEquations2D(1.4)

        semi = build_solver(Val{approximation_type}(), polydeg, cells_per_dimension,
                            initial_condition, equations)

        tspan = (0.0, 10.0)
        ode = semidiscretize(semi, tspan)

        alive_callback = AliveCallback(alive_interval = 10)

        integrator = init(ode, RDPK3SpFSAL49(), abstol=abstol, reltol=reltol,
                          dt = 1e-3 * estimate_dt(semi.mesh, semi.solver),
                          save_everystep=false, callback=alive_callback)
        try
            solve!(integrator)
        catch e
            @warn "Blow-up" e
        end
        integrator.t
    end
end

polydegs = (3, 7)
mesh_resolutions = (32, 16)

approximation_types = (:dgsem, :gauss, :hybrid, :staggered)
atwood_numbers = range(0.1, 0.9, length=20)

parameters = [(approximation_type, polydeg, cells_per_dimension, A)
               for (polydeg, cells_per_dimension) in zip(polydegs, mesh_resolutions)
               for approximation_type in approximation_types for A in atwood_numbers]

blow_up_times = pmap(args -> blow_up_time(args...), parameters)

data = DataFrame(approximation_type = map(x -> x[1], parameters),
                 polydeg = map(x -> x[2], parameters),
                 cells_per_dimension = map(x -> x[3], parameters),
                 atwood_number = map(x -> x[4], parameters),
                 blow_up_time = blow_up_times)

CSV.write(joinpath(@__DIR__, "compare_entropy_projection_variants.txt"), data)

# extra code to plot snapshots of density for Atwood number 3/7
function plot_sol(approximation_type, polydeg, cells_per_dimension, Atwood_number=3/7;
                  abstol=1.0e-7, reltol=1.0e-7)

    println("Running with approximation_type=$(approximation_type), polydeg=$polydeg, cells_per_dimension=$cells_per_dimension, Atwood=$Atwood_number")

    initial_condition = set_initial_condition(Atwood_number)
    equations = CompressibleEulerEquations2D(1.4)

    semi = build_solver(Val{approximation_type}(), polydeg, cells_per_dimension,
                        initial_condition, equations)

    tspan = (0.0, 5.0)
    ode = semidiscretize(semi, tspan)

    alive_callback = AliveCallback(alive_interval = 10)

    integrator = init(ode, RDPK3SpFSAL49(), abstol=abstol, reltol=reltol,
                        dt = 1e-3 * estimate_dt(semi.mesh, semi.solver),
                        save_everystep=false, callback=alive_callback)
    sol = solve!(integrator)

    fig, ax, plt = plot(Trixi.PlotData2D(sol)["rho"], plot_mesh=false)
    Colorbar(fig[1,2], plt)
    ax.aspect=DataAspect(); Makie.tightlimits!(ax); hidedecorations!(ax); hidespines!(ax)
    set_theme!(Theme(fontsize=20))
    save("khi_density_$(string(approximation_type))_p$(polydeg)_$(cells_per_dimension)_cells_t5.png", fig.scene)

    return sol
end
# plot_sol(:gauss, 3, 64)
# plot_sol(:hybrid, 3, 64)
# plot_sol(:staggered, 3, 64)


text:
function α(x) #clearing error function. To be replaced with solution to issue 1
      sum(x)
  end

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
#!/usr/bin/env julia

"""
bam_record.jl - Houses functions that deal with BAM Record objects and properties

By: Shaun Adkins (sadkins@som.umaryland.edu)
"""

abstract type AbstractAlignment end

struct FragmentAlignment<:AbstractAlignment
    count_multiplier::Float32
end
struct ReadAlignment<:AbstractAlignment
    count_multiplier::Float32
end

FragmentAlignment() = FragmentAlignment(one(Float32))
ReadAlignment() = ReadAlignment(one(Float32)/2)

multiplier(f::T) where {T<:AbstractAlignment} = f.count_multiplier
#multiplier(f::ReadAlignment) = f.count_multiplier

is_reversestranded(strand_type::String) = return (strand_type == "reverse" ? true : false)

#isduplicate(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_DUP == 0x0400
#ismatereverse(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_MREVERSE == 0x0020
isproperpair(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_PROPER_PAIR == 0x0002
isread1(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_READ1 == 0x0040
isread2(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_READ2 == 0x0080
isreverse(record::BAM.Record) = BAM.flag(record) & SAM.FLAG_REVERSE == 0x0010

# TODO: Checking specific case where dovetailed reads where the two mates completely overlap each other
iscompleteoverlap(record::BAM.Record) = get_alignment_start_end(record::BAM.Record, true) == get_alignment_start_end(record::BAM.Record, false)

# Forward-stranded assay
## Positive Strand:
### R1 - forward (65), R2 - revcom (145)
## Negative Strand
### R1 - revcom (81), R2 - forward (129)
# Reverse-stranded assay (i.e. Illumina) have these strands flipped

### For non-properly paired reads
# Singletons will either have flag 0 or 16

# Forward strand flags
flag65(record::BAM.Record) = isread1(record) && !isreverse(record)
flag145(record::BAM.Record) = isread2(record) && isreverse(record)
# Reverse strand flags
flag81(record::BAM.Record) = isread1(record) && isreverse(record)
flag129(record::BAM.Record) = isread2(record) && !isreverse(record)

function assign_read_to_strand(record::BAM.Record, reversestrand::Bool=false)
    """Use the bitwise flags to assign the paired read to the correct strand."""

    positiveflags = Set{Bool}([flag145(record), flag65(record)])
    negativeflags = Set{Bool}([flag81(record), flag129(record)])
    if reversestrand
        positiveflags = Set{Bool}([flag81(record), flag129(record)])
        negativeflags = Set{Bool}([flag145(record), flag65(record)])
    end

    # Does record flag belong in the pos or neg set?
    any(positiveflags) && return '+'
    any(negativeflags) && return '-'

    # Read must be a singleton
    if isreverse(record)
        # If reverse-stranded, strand is flipped
        reversestrand && return '+'
        return '-'
    end
    reversestrand && return '-'
    return '+'
end

#function determine_record_type(record::BAM.Record, max_frag_size::UInt)
function canbefragment(record::BAM.Record, max_frag_size::UInt)
    """Determine if BAM Record is a read (R) or fragment (F)."""
    # True = Can validate as a fragment
    # False = Can validate as a read
    if validate_fragment(record) && is_templength_smaller_than_max_fragment_size(BAM.templength(record), max_frag_size)
        return true
    elseif validate_read(record, max_frag_size)
        return false
    end
    # For fragments, only one read is looked at.  The other one is essentially skipped to avoid overcounting.
    # Reads that also fail validation go here.
    return nothing
end

function get_alignment_coords_set(alignment::Interval)
    """Get the range of coordinates for this alignment, returned as a Set."""
    return Set{UInt}(leftposition(alignment) : rightposition(alignment))
end

function get_alignment_interval(record::BAM.Record, max_frag_size::UInt, reversestrand::Bool=false)
    """Return an alignment-based Interval for the current record."""
    isfragment = canbefragment(record, max_frag_size)
    isfragment === nothing && return nothing  # Skip alignments not being considered
    return Interval(BAM.refname(record),
                    get_alignment_start_end(record, isfragment),
                    assign_read_to_strand(record, reversestrand),
                    isfragment
    )
end

function get_alignment_interval(record::BAM.Record, max_frag_size::UInt, strand_type::String)
    """Return an alignment-based Interval for the current record."""
    return get_alignment_interval(record, max_frag_size, is_reversestranded(strand_type))
end

function get_alignment_start_end(record::BAM.Record, isfragment::Bool)
    """ Get the start and end coordinates of the alignment record."""
    return get_alignment_start_end(record, getalignmenttype(isfragment))
end

function get_alignment_start_end(record::BAM.Record, record_type::FragmentAlignment)
    """Get the start and end coordinates of the fragment.  Assuming the use of the negative template"""
    return BAM.nextposition(record):BAM.rightposition(record)
end

function get_alignment_start_end(record::BAM.Record, record_type::ReadAlignment)
    """Get the start and end coordinates of the read."""
    return BAM.position(record):BAM.rightposition(record)
end

function getstrand(interval::Interval, stranded::Bool)
    """Get strand of interval with respect to strandedness arguments."""
    strand = '+'
    if stranded
        strand = convert(Char, interval.strand)
    end
    return strand
end

function getalignmenttype(isfragment::Bool)
    return (isfragment ? FragmentAlignment() : ReadAlignment())
end

function ismultimapped(record::BAM.Record)
    """Test to see if alignment is multimapped across multiple regions of the genome."""
    try
        return haskey(record, "NH") && record["NH"] > 1
    catch
        # Ran into bug where some attributes are 'X' type which is not valid
        # Bioalignments BAM.Record.auxdata(record) throws LoadError for these
        return false
    end
end

function is_templength_negative(templength::Int64)
    """Check to see if the read template is going in the opposite version."""
    return templength < 0
end

function is_templength_smaller_than_max_fragment_size(templength::Int64, max_frag_size::UInt)
    """Check to see if fragment template length is smaller than specified maximum fragment size."""
    return abs(templength) <= max_frag_size
end

function validate_fragment(record::BAM.Record)
    """Ensure alignment record can be used to calculate fragment depth. (Only need negative template of fragment).  In the case of completely dovetailling reads, just take the first read."""
    if isproperpair(record)
        if iscompleteoverlap(record)
            isread1(record) && return true
        elseif is_templength_negative(BAM.templength(record))
            return true
        end
    end
    return false
end

function validate_read(record::BAM.Record, max_frag_size::UInt)
    """Ensure alignment record can be used to calculate read depth."""
    return !(isproperpair(record) && is_templength_smaller_than_max_fragment_size(BAM.templength(record), max_frag_size))
end

max_stars_repo_name: IGS/FADU

max_issues_repo_name: adkinsrs/FADU

max_forks_repo_name: adkinsrs/FADU


text:
using LightNLP
using LightNLP.Tokenization
using JLD2, FileIO

training = true

if training
    embedsfile = ".data/glove.6B.100d.h5"
    trainfile = ".data/eng.train.BIOES"
    testfile = ".data/eng.testb.BIOES"
    nepochs = 20
    learnrate = 0.0005
    batchsize = 10
    ner = NER.Decoder(embedsfile, trainfile, testfile, nepochs, learnrate, batchsize)
    save("ner.jld2", "ner", ner)
else
    ner = load("ner.jld2", "ner")
    testfile = ".data/eng.train.BIOES"
    decode(ner, testfile)
end


text:
function getDivGrad(n1,n2,n3)
	# D = getDivGrad(n1,n2,n3)
	# builds 3D divergence operator
	D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
	D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
	D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))

	Div = [D1 D2 D3]
	return Div*Div'
end

function ddx(n)
# generate 1D finite difference on staggered grid
	return d = spdiags(ones(n)*[-1 1],[0,1],n,n+1)
end

function spdiags(B,d,m,n)
# A = spdiags(B,d,m,n)
# creates a sparse matrix from its diagonals
	d = d[:]
	p = length(d)
	len = zeros(p+1,1)
	for k = 1:p
		len[k+1] = round(Int,len[k]+length(max(1,1-d[k]):min(m,n-d[k])))
	end
	a = zeros(round(Int,len[p+1]),3)
	for k = 1:p
		# Append new d[k]-th diagonal to compact form
		i = max(1,1-d[k]):min(m,n-d[k])
		a[(round(Int,len[k])+1):round(Int,len[k+1]),:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
	end

	A = sparse(round(Int,a[:,1]),round(Int,a[:,2]),a[:,3],m,n)
	return A
end


text:
using Continuables
using Test

@testset "utils" begin
  include("utils.jl")
end

# check that @Ref is working correctly
# ====================================

@testset "Ref" begin
  expr1_simple = @macroexpand @Ref begin
    a = Ref(1)
    f(a) = 2 * a
    a + f(a)
  end
  expr2_simple = quote
    a = Ref(1)
    f(a) = 2 * a
    a.x + f(a.x)
  end
  @test same_expr(expr1_simple, expr2_simple)

  # difficult case to parse:
  expr1_complex = @macroexpand @Ref begin
    a = Ref(Some(1))
    b = Ref("hi")
    f(a::Int = b; c = b, d = 4) = 2 * a
    a.value + f(a.value)
  end
  expr2_complex = quote
    a = Ref(Some(1))
    b = Ref("hi")
    f(a::Int = b.x; c = b.x, d = 4) = 2 * a
    a.x.value + f(a.x.value)
  end
  @test same_expr(expr1_complex, expr2_complex)
end

@testset "specializing functions for Continuable" begin
  mycont(x) = @cont cont(x)
  Base.length(::Continuable{<:innerfunctype(mycont(:examplesignature))}) = :justatest
  @test length(mycont(1)) == :justatest
  @test length(mycont("a")) == :justatest
  # other continuables are not affected
  anothercont(x) = @cont cont(x)
  @test length(anothercont(42)) == 1
end

@testset "constructing Continuables" begin

  cont1 = @cont cont(1)
  @test collect(cont1) == [1]
  @test collect(singleton(3)) == [3]

  @test collect(take(repeated(() -> 4), 3)) == [4, 4, 4]
  @test collect(repeated(() -> 4, 3)) == [4, 4, 4]

  @test collect(take(3, iterated(x -> x*x, 2))) == [2, 4, 16]
end

@testset "DropWhile TakeWhile iterables" begin
  @test length(dropwhile(x -> x<4, 1:10)) == 7
  @test Base.IteratorSize(typeof(dropwhile(x -> x<4, 1:10))) == Base.SizeUnknown()

  @test eltype(dropwhile(x -> x<4, 1:10)) == Int
  @test Base.IteratorEltype(typeof(dropwhile(x -> x<4, 1:10))) == Base.HasEltype()

  @test length(takewhile(x -> x<4, 1:10)) == 3
  @test Base.IteratorSize(typeof(takewhile(x -> x<4, 1:10))) == Base.SizeUnknown()

  @test eltype(takewhile(x -> x<4, 1:10)) == Int
  @test Base.IteratorEltype(typeof(takewhile(x -> x<4, 1:10))) == Base.HasEltype()
end

# Check Continuables standard Interface
# =====================================

@testset "standard interface" begin
  @test_throws ErrorException cont(1)

  @test collect(@i2c 2:4:10) == collect(2:4:10)
  @test collect(i2c(2:4:10), length(2:4:10)) == collect(2:4:10)

  @test collect(enumerate(@i2c 2:4:10)) == collect(enumerate(2:4:10))
  @test collect(filter(x -> x < 7, @i2c 2:4:10)) == collect(filter(x -> x < 7, 2:4:10))

  @test collect(map(x->x^2, @i2c 2:4:10)) == collect(map(x->x^2, 2:4:10))
  @test collect(take(cycle(@i2c 1:3), 11)) == collect(take(cycle(1:3), 11))
  @test collect(take(cycle(@i2c 1:3), 11)) == collect(take(cycle(1:3), 11))
  @test collect(cycle(i2c(1:3), 3)) == [1,2,3, 1,2,3, 1,2,3]

  @test collect(drop(i2c(1:10), 3)) == collect(drop(1:10, 3))
  @test collect(drop(3, i2c(1:10))) == collect(drop(1:10, 3))

  @test nth(i2c(0:10), 8) == (0:10)[8]
  @test_throws ErrorException nth(12, i2c(0:10))

  @test Base.IteratorSize(@i2c 1:4) == Base.SizeUnknown()
  @test length(@i2c 3:20) == 18

  @test Base.IteratorEltype(emptycontinuable) == Base.EltypeUnknown()
  @test Base.eltype(@i2c 1:4) == Any

  test_continuable_product = @i2c 1:10
  @test product(test_continuable_product) === test_continuable_product  # should be a no-op
  @test collect(product(i2c(1:10), i2c(1:3))) == [(i,j) for i in 1:10 for j in 1:3]
  @test collect(product(i2c(1:10), i2c(1:3), i2c(2:4))) == [(i,j,k) for i in 1:10 for j in 1:3 for k in 2:4]

  # however they both are not easily comparable... because transpose does not work on arrays of tuples right now...

  @test collect(partition(i2c(1:15), 4)) == collect(partition(1:15, 4))
  @test collect(partition(i2c(1:20), 3, 5)) == [
    Any[1, 2, 3],
    Any[6, 7, 8],
    Any[11, 12, 13],
    Any[16, 17, 18],
  ]
  @test collect(partition(i2c(1:10), 5, 2)) == Any[
    Any[1, 2, 3, 4, 5],
    Any[3, 4, 5, 6, 7],
    Any[5, 6, 7, 8, 9],
  ]

  @test collect(zip(i2c(1:10), i2c(4:13))) == collect(zip(1:10, 4:13))
  @test collect(zip(i2c(1:10), i2c(4:13), lazy = false)) == collect(zip(1:10, 4:13))

  @test any(x -> x == 4, @i2c 1:10) == any(x -> x == 4, 1:10)
  @test any(x -> x == 4, i2c(1:10), lazy=false) == any(x -> x == 4, 1:10)
  @test any(x -> x == 4, i2c(1:3)) == any(x -> x == 4, 1:3)
  @test any(x -> x == 4, i2c(1:3), lazy=false) == any(x -> x == 4, 1:3)

  @test all(x -> x <= 4, i2c(1:10)) == all(x -> x <= 4, 1:10)
  @test all(x -> x <= 4, i2c(1:10), lazy=false) == all(x -> x <= 4, 1:10)
  @test all(x -> x <= 4, i2c(1:3)) == all(x -> x <= 4, 1:3)
  @test all(x -> x <= 4, i2c(1:3), lazy=false) == all(x -> x <= 4, 1:3)

  @test reduce((acc, x) -> acc + x*x, @i2c 1:10) == reduce((acc, x) -> acc + x*x, 1:10)
  @test foldl((acc, x) -> acc + x*x, @i2c 1:10) == foldl((acc, x) -> acc + x*x, 1:10)
  @test reduce((acc, x) -> acc + x*x, i2c(1:10), init = 100) == reduce((acc, x) -> acc + x*x, 1:10, init = 100)
  @test reduce!(push!, i2c(1:10), init = []) == collect(1:10)

  @test sum(@i2c 1:10) == sum(1:10)
  @test prod(@i2c 1:10) == prod(1:10)

  @test collect(chain(@i2c(1:10), @i2c(3:5))) == collect(chain(1:10, 3:5))
  contcont = @cont begin
    for i in 1:10
      cont(@i2c i:10)
    end
  end
  @test collect(flatten(contcont)) == collect(flatten(i:10 for i in 1:10))

  @test collect(takewhile(x -> x <= 4, @i2c 1:10)) == collect(takewhile(x -> x <= 4, 1:10))
  @test collect(dropwhile(x -> x <= 4, @i2c 1:10)) == collect(dropwhile(x -> x <= 4, 1:10))

  @test groupby(x -> x % 4, @i2c 1:10) == groupby(x -> x % 4, 1:10)
  @test groupbyreduce(x -> x % 4, @i2c(1:10), (x, y) -> x + y, x -> x+5) == groupbyreduce(x -> x % 4, 1:10, (x, y) -> x + y, x -> x+5)
end


text:
using DDF

using DifferentialForms
using StaticArrays
using Test

@testset "Circumcentre D=$D N=$N" for D in 1:Dmax, N in 1:(D + 1)
    T = Rational{BigInt}
    xs = SVector{N}(rand(Form{D,1,T}) for n in 1:N)
    cc = circumcentre(xs)
    cc::Form{D,1,T}
    # Check that circumcentre has same distance from all vertices
    r2 = norm2(xs[1] - cc)
    @test all(==(r2), norm2(xs[n] - cc) for n in 1:N)
    # Check that circumcentre lies in hyperplane spanned by vertices
    if sum(frank.(xs)) > D
        @test true
    else
        subspace = ∧(xs...)
        if frank(cc) + frank(subspace) > D
            @test true
        else
            @test iszero(cc ∧ subspace)
        end
    end
end

@testset "Weighted circumcentre D=$D N=$N" for D in 1:Dmax, N in 1:(D + 1)
    T = Rational{BigInt}
    xs = SVector{N}(rand(Form{D,1,T}) for n in 1:N)
    ws = SVector{N}(rand(Form{D,0,T}) for n in 1:N)
    cc = circumcentre(xs, ws)
    cc::Form{D,1,T}
    # Check that circumcentre has correct weighted distance from all vertices
    r2 = norm2(xs[1] - cc) - ws[1][]
    @test all(==(r2), norm2(xs[n] - cc) - ws[n][] for n in 1:N)
    # Check that circumcentre lies in hyperplane spanned by vertices
    if sum(frank.(xs)) > D
        @test true
    else
        subspace = ∧(xs...)
        if frank(cc) + frank(subspace) > D
            @test true
        else
            @test iszero(cc ∧ subspace)
        end
    end
end

@testset "Volume D=$D N=$N" for D in 1:Dmax, N in 1:(D + 1)
    T = Rational{BigInt}

    x0 = zero(Form{D,1,T})
    xs = SVector{N}(x0, (x0 + unit(Form{D,1,T}, n) for n in 1:(N - 1))...)
    v = volume(xs)
    @test v == T(1) / factorial(N - 1)

    xs = SVector{N}(x0,
                    (x0 + unit(Form{D,1,T}, n - (N - 1) + D) for n in 1:(N - 1))...)
    v = volume(xs)
    @test v == T(1) / factorial(N - 1)

    if N ≥ 2
        xs = SVector{N}(xs[2], xs[1], xs[3:end]...)
        v = volume(xs)
        @test v == T(1) / factorial(N - 1)
    end

    if N ≥ 2
        xs = SVector{N}(xs[2:end]..., xs[1])
        v = volume(xs)
        @test v == T(1) / factorial(N - 1)
    end

    if N > 1
        xs = SVector{N}(x0,
                        (x0 + unit(Form{D,1,T}, n - (N - 1) + D) for n in
                                                                     1:(N - 1))...)
        xs = SVector{N}(xs[1:(end - 1)]..., 5 * xs[end])
        v = volume(xs)
        @test v == T(5) / factorial(N - 1)
    end

    x1 = rand(Form{D,1,T})
    xs = SVector{N}(x1,
                    (x1 + unit(Form{D,1,T}, n - (N - 1) + D) for n in 1:(N - 1))...)
    v = volume(xs)
    @test v == T(1) / factorial(N - 1)

    xs = 2 * xs
    v = volume(xs)
    @test v == T(2^(N - 1)) / factorial(N - 1)
end

#TODO @testset "Dual volume" begin end

@testset "Morton ordering D=$D" for D in 1:Dmax
    # nbits = 10 ÷ D
    nbits = 2 ÷ D
    imax = (1 << nbits) - 1
    found = falses(1 << (D * nbits))
    for i in
        CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> imax, D))
        u = SVector{D,UInt64}(i.I)
        m = morton(u)
        @test !found[m + 1]
        found[m + 1] = true
    end

    T = Float64
    xmin = SVector{D,T}(-1 for d in 1:D)
    xmax = SVector{D,T}(+1 for d in 1:D)
    dx = (xmax - xmin) / (imax .+ 1)
    ifound = UInt64[]
    for i in
        CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> imax, D))
        u = SVector{D,UInt64}(i.I)
        # x = xmin + u .* dx + (T(0.01) .+ T(0.98) * rand(SVector{D,T}) .* dx)
        x = xmin + (u + rand(SVector{D,T})) .* dx
        @assert all(xmin .≤ x .≤ xmax)
        m = morton(x, xmin, xmax)
        push!(ifound, m)
    end
    @test allunique(sort!(ifound))
end


text:
"""
    viral_load_infectivity_testpos.jl

Written by Carl Whitfield, University of Manchester, 2021
"""

#============== Options (uncomment and change here) ==============#
#VL_model = ke_model_no         #choice of viral load model (Ke et al is default)
#LFD_model = SC_data_2021         #choice of LFD sensitivity model (social care binned 2021 data is default)
#PCR_sens_max = 0.83            #max PCR sensitivity (Ferretti et al 2021)
#Inf_model = ke_inf_model_no    #infectivity model
#p_asymp = 0.5                  #asymptomatic fraction
#onset_opt = V0_model_opt        #IF using kissler VL model, option for how to generate onset time
#peak_inf_opt = marks_peakinf_opt   #IF using flat or linear inf model, what function is used for peak value
#PCR_TaT_scale = 1.0
#===============================================================#

include("definitions.jl")

function generate_VL_params(Asymp::Bool)
    if VL_model == kissler_model_no
        return generate_VL_params_Kissler(Asymp)
    elseif VL_model == ke_model_no
        return generate_ke_VL_params()
    end
end

"""
    generate_symp_time(tp::Float64)

Generate and return single time of symptom onset,
bounded by the times where viral load is at detectable levels

## Arguments:
`tp` = Time from infection to peak VL (days)

## Returns:
`Float64` = Symptom onset time (days)

## See also:
`generate_VL_params(Asymp::Bool)`
"""
function generate_symp_time(tp::Float64)
    #Truncate in narrow window, 2 days either side of VL peak
    Strunc = truncated(Gamma(symp_alpha,1.0/symp_beta),
                       max(0.0, tp - 2.0), tp + 2.0)
    return rand(Strunc)
end

"""
    generate_isolations(Ntot::Int, Pisol::Float64)

Randomly generate whether people obey symptomatic isolation or not

## Arguments:
`Ntot` = Number to generate (indices are 1:Ntot)
`Pisol` = Probability of isolation

## Returns:
`Array{Int,1}` = Array of indices of those who will isolate at symptom
                 onset (if symptomatic)
"""
function generate_isolations(Ntot::Int, Pisol::Float64)
    return randsubseq(1:Ntot, Pisol)
end

function infectivity(PVL::Float64, r::Float64, d::Float64, tp::Float64,
                     T::Int64)
    if Inf_model == ke_inf_model_no
        J,h = generate_ke_inf_params()
        inf = infectivity_Ke(PVL, r, d, tp, J, h, T)
    elseif Inf_model == flat_inf_model_no
        PInf = generate_peak_infectivity(PVL)
        #divide by 2 so area under inf curve same as linear case
        inf = infectivity_flat(PVL, r, d, tp, 0.5*PInf, T)
    elseif Inf_model == linear_inf_model_no
        PInf = generate_peak_infectivity(PVL)
        inf = infectivity_linear(PVL, r, d, tp, PInf, T)
    end

    return inf
end

"""
    build_viral_load_distribution!(v::Array{Float64,1}, inf::Array{Float64,1})

Generate a viral load and infectivity trajectory

## Arguments:
`sim` = Dict container with simulation info

` i` = Index of individual to build

## See also:
`build_viral_load_distributions!(sim::Dict)`
"""
function build_viral_load_distribution!(sim::Dict, index::Int64)
    Asymp = generate_asymptomatic()
    PVL, tp, r, d = generate_VL_params(Asymp)

    T = Int64(ceil(tp + log(10)*(PVL - VL_LOD_PCR)/d))
    if T < 1
        print("Strange parameter set generated: PVL = ", PVL,
              ", d = ", d, ", tp = ", tp, '\n')
        T = 1
    end
    v = zeros(T)
    i = 1:T                              #day index
    t = i .- 1                           #days since infection

    v[t .<= tp] = PVL .+  (r/log(10))*(t[t .<= tp] .- tp)
    v[t .> tp] = PVL .- (d/log(10))*(t[t .> tp] .- tp)

    ST = generate_symp_time(tp)
    inf = infectivity(PVL, r, d, tp, T)
    SD = Int64(round(ST))

    sim["symp_day"][index] = SD
    sim["VL_mag"][index] = PVL
    sim["asymptomatic"][index] = Asymp
    sim["VL_profiles"][index] = v
    sim["infection_profiles"][index] = inf
end

"""
    build_viral_load_distributions!(sim::Dict)

Generate multiple viral load and infectivity trajectories.

## Arguments:
`sim` = Dict container with the following keys:

    "VL_profiles"=>Array{Array{Float64,1},1} = Each entry i is
        a container for a generated viral load trajectory for
        individual i, interpolated to one value per day.

    "infection_profiles"=>Array{Array{Float64,1},1} = Each
        entry i is a container for a generated infectivity profile
        of individual i, interpolated to one value per day.

    "symp_day"=>Array{Int64,1} = Each entry i is the symptom
        onset time of individual i.

    "VL_mag"=>Array{Float64,1} = Each entry i is the peak viral
        load of individual i.

    "asymptomatic"=>Array{Float64,1}  = Each entry i is the
        asymptomatic status of individual i.

    All arrays must have the same size (number of trajectories
    to be generated), and the entries will be overwritten with
    the generated values.

## Returns:
none

## See also
`init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)`
"""
function build_viral_load_distributions!(sim::Dict)
    build_viral_load_distribution!.(Ref(sim), collect(1:sim["Ntot"]))
end

"""
    init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)

Randomly generate all viral load and infectivity data for a simulation

## Arguments:
`Ntot` = Number of individuals to simulate
`Pisol` = Probability of self-isolation at symptom onset

## Returns:
`Dict` = Container with generated parameters including:

        "Ntot" => `Int` = Number of individuals in simulation (indexed 1:Ntot)

        "asymptomatic" => `Array{Bool,1}` = Asymptomatic status for each individual

        "isolation_time" => `Array{Int64,1}` = Day of isolation for each individual
                                           (to be filled later)

         "will_isolate" => `Array{Bool,1}` = Indicator of whether an individual will
                                         isolate at symptom onset time

         "VL_mag" => `Array{Float64,1}` = Peak viral load for each individual

         "inf_mag" => `Array{Float64,1}` = Peak infectivity for each individual

         "infection_profiles" => `Array{Array{Float64,1},1}` = infectivity for each day
                            since infection, for each individual

         "VL_profiles" => `Array{Array{Float64,1},1}` = viral load for each day
                            since infection, for each individual

         "days_infectious" => `Array{Int64,1}` = Number of days spent with infectivity
                                               > 0 for each individual

          "symp_day" => `Array{Int64,1}` = Day at which individuals would isolate due to
                                         symptom onset

          "non_isolators" => `Array{Int64,1}` = Array of indices ofindividuals who would
                                          refuse to isolate if they developed symptoms
"""
function init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)
    #simulate Ntot people with baseline isolation probability Pisol
    sim = Dict("Ntot"=>Ntot, "asymptomatic"=>zeros(Bool, Ntot),
               "VL_mag"=>zeros(Float64,Ntot),
               "symp_day"=>zeros(Int64,Ntot),
               "isolation_time"=>zeros(Int64, Ntot),
               "will_isolate"=>zeros(Bool, Ntot),
               "inf_mag"=>zeros(Float64, Ntot),
               "infection_profiles"=>Array{Array{Float64,1},1}(undef,Ntot),
               "VL_profiles"=>Array{Array{Float64,1},1}(undef,Ntot),
               "days_infectious" => zeros(Int64,Ntot))
#     for i in 1:Ntot
#         push!(sim["infection_profiles"],zeros(0))
#         push!(sim["VL_profiles"],zeros(0))
#     end
    build_viral_load_distributions!(sim)
    sim["will_isolate"][generate_isolations(Ntot, Pisol)] .= true
    nr = 1:sim["Ntot"]
    sim["non_isolators"] = nr[(sim["will_isolate"] .== false)]
    sim["will_isolate"][sim["asymptomatic"]] .= false #asymptomatics don't self isolate, but are not "non isolators"

    return sim
end

"""
    plot_infection_profiles(sim::Dict)

Create a plot of all infectiousness profiles

## Arguments:
`sim` = Dict created by `init_VL_and_infectiousness`

## Returns:

## See also:
`init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)`
"""
function plot_infection_profiles(sim::Dict)
    if sim["Ntot"] > 0
        pp = Plots.plot(sim["infection_profiles"][1], xlabel = "Days", ylabel = "Infectiousness")
        for i in 2:sim["Ntot"]
            pp = Plots.plot!(sim["infection_profiles"][i])
        end
    end

    pout = Plots.display(pp)

    return pout
end

function plot_trajectories(A::Array{Array{Float64,1},1})
    Ntraj = length(A)
    tmax = max(length.(A)...)
    traj = zeros(tmax,Ntraj)
    for n in 1:Ntraj
        traj[1:length(A[n]),n] = A[n]
    end
    npick = collect(1:Ntraj)
    if Ntraj > 100
        npick = rand(npick,100)
    end

    plot(1:tmax, traj[:,npick],label=:none, color=:grey, alpha=0.3)
    plot!(1:tmax, median(traj,dims=2), style=:dashdot, color=:black,
          linewidth=3, label="median")
    p1 = plot!(1:tmax, mean(traj,dims=2), style=:dot, color=:red, linewidth=3,
        xlabel="Days since infection", label="mean")
    return p1
end


text:
const ImmutableContainer = Union{
    ImmutableDict,
    NamedTuple,
    Tuple,
}

push(xs, i1, i2, items...) =
    foldl(push, items, init=push(push(xs, i1), i2))

push(xs::AbstractVector, x) = vcat(xs, SingletonVector(x))
push(xs::AbstractSet, x) = union(xs, SingletonVector(x))
push(xs::AbstractDict, x::Pair) = merge(xs, SingletonDict(x[1], x[2]))

push(xs::Tuple, items...) = (xs..., items...)

push(xs::NamedTuple{names}, x::Pair{Symbol}) where {names} =
    NamedTuple{(names..., x.first)}((xs..., x.second))

push(xs::NamedTuple{names}, x::Pair{Val{name}}) where {names, name} =
    NamedTuple{(names..., name)}((xs..., x.second))

push(::NamedTuple, x) =
    error("`push(::NamedTuple, x::$(typeof(x)))` is not supported.\n",
          "Use `push(::NamedTuple, :NAME => x)` or ",
          "`push(::NamedTuple, Val(:NAME) => x)`.")

push(xs::ImmutableDict, x::Pair) = ImmutableDict(xs, x)

append(xs, ys) = _append(xs, ys)
_append(xs, ys) = append!(copy(xs), ys)
_append(xs, ys::Tuple) = push(xs, ys...)
_append(xs, ys::Pairs{Symbol, <:Any, <:Any, <:NamedTuple}) = push(xs, ys...)

append(xs::AbstractVector, ys::AbstractVector) = vcat(xs, ys)

append(xs::ImmutableContainer, ys) = push(xs, ys...)

append(xs::AbstractString, ys::AbstractString) = string(xs, ys)

function pushfirst(xs::AbstractVector, ys...)
    T = promote_type(eltype(xs), map(typeof, ys)...)
    zs = similar(xs, T, length(xs) + length(ys))
    copyto!(zs, length(ys) + 1, xs, 1, length(xs))
    copyto!(zs, 1, ys, 1, length(ys))
    return zs
end

pushfirst(xs::Tuple, ys...) = (ys..., xs...)

@inline splice(xs::Tuple, i::Integer) = _splice_del((), i - 1, xs...)
@inline _splice_del(pre, n, x, xs...) =
    if length(pre) < n
        _splice_del((pre..., x), n, xs...)
    else
        (pre..., xs...), x
    end

@inline splice(xs::Tuple, i::Integer, v) = _splice_insert((), i - 1, v, xs...)
@inline _splice_insert(pre, n, v, x, xs...) =
    if length(pre) < n
        _splice_insert((pre..., x), n, v, xs...)
    else
        (pre..., v, xs...), x
    end

pop(xs::Tuple) = xs[1:end-1], xs[end]
pop(xs::NamedTuple{names}) where {names} =
    NamedTuple{names[1:end-1]}(Tuple(xs)[1:end-1]), xs[end]

maybepop(xs::NamedTuple, key) = maybepop(xs, Val(key))
maybepop(xs::NamedTuple{names}, ::Val{key}) where {names, key} =
    if haskey(xs, key)
        Some((Base.structdiff(xs, NamedTuple{(key,)}), xs[key]))
    else
        nothing
    end

pop(xs::NamedTuple, key, default) =
    something(maybepop(xs, key), Some((xs, default)))

function pop(xs::NamedTuple, key)
    result = maybepop(xs, key)
    result === nothing && throw(KeyError(key))
    return something(result)
end

popfirst(xs::Tuple) = xs[2:end], xs[1]
popfirst(xs::NamedTuple{names}) where {names} =
    NamedTuple{names[2:end]}(Tuple(xs)[2:end]), xs[1]

delete(xs::NamedTuple, key) = something(maybepop(xs, key), (xs,))[1]

_empty(xs) = empty(xs)
if VERSION < v"1.3.0-DEV.533"
    _empty(xs::NamedTuple) = NamedTuple()
end

_setindex(xs, v, I...) = Base.setindex(xs, v, I...)
_setindex(xs::NamedTuple, value, name) = setproperty(xs, name, value)

function _setindex(xs::AbstractArray, v, I...)
    T = promote_type(eltype(xs), typeof(v))
    ys = similar(xs, T)
    if eltype(xs) !== Union{}
        copy!(ys, xs)
    end
    ys[I...] = v
    return ys
end

function _setindex(d0::AbstractDict, v, k)
    K = promote_type(keytype(d0), typeof(k))
    V = promote_type(valtype(d0), typeof(v))
    d = Dict{K, V}()
    copy!(d, d0)
    d[k] = v
    return d
end

setproperty(value, name, x) = setproperties(value, NamedTuple{(name,)}((x,)))

materialize(::Any, x) = Broadcast.materialize(x)


text:
## SEXPREC types as defined in ENV["R_INCLUDE_DIR"]/Rinternals.h

const NILSXP     =  0 # NULL in R
const SYMSXP     =  1 # R Symbol
const LISTSXP    =  2 # internal "pairs list", the R list type is 19
const CLOSXP     =  3 # closures
const ENVSXP     =  4 # environments
const PROMSXP    =  5 # promises: unevaluated closure arguments
const LANGSXP    =  6 # language constructs (special lists
const SPECIALSXP =  7 # special forms
const BUILTINSXP =  8 # builtin non-special forms
const CHARSXP    =  9 # "scalar" string type (internal only
const LGLSXP     = 10 # logical vectors
# 11 and 12 were factors and ordered factors in the 1990s
const INTSXP     = 13 # integer vectors
const REALSXP    = 14 # real variables
const CPLXSXP    = 15 # complex variables
const STRSXP     = 16 # string vectors
const DOTSXP     = 17 # dot-dot-dot object
const ANYSXP     = 18 # make "any" args work.
const VECSXP     = 19 # generic vectors
const EXPRSXP    = 20 # expressions vectors
const BCODESXP   = 21 # byte code
const EXTPTRSXP  = 22 # external pointer
const WEAKREFSXP = 23 # weak reference
const RAWSXP     = 24 # raw bytes
const S4SXP      = 25 # S4, non-vector


text:
using .CUDA
import .CUDA: CuDevice, CuContext, devices, attribute

export CuArrayDeviceProc

"Represents a single CUDA GPU device."
struct CuArrayDeviceProc <: Dagger.Processor
    owner::Int
    #ctx::CuContext
    device::Int
end
@gpuproc(CuArrayDeviceProc, CuArray)
#= FIXME: DtoD copies and CUDA IPC
function Dagger.move(from::CuArrayDeviceProc, to::CuArrayDeviceProc, x)
    if from === to
        return x
    else
        error("Not implemented")
    end
end
=#
function Dagger.execute!(proc::CuArrayDeviceProc, func, args...)
    fetch(Threads.@spawn begin
        task_local_storage(:processor, proc)
        CUDA.device!(proc.device)
        CUDA.@sync func(args...)
    end)
end
Base.show(io::IO, proc::CuArrayDeviceProc) =
    print(io, "CuArrayDeviceProc on worker $(proc.owner), device $(proc.device)")

processor(::Val{:CUDA}) = CuArrayDeviceProc
cancompute(::Val{:CUDA}) = CUDA.has_cuda()
kernel_backend(::CuArrayDeviceProc) = CUDADevice()

if CUDA.has_cuda()
    for dev in devices()
        Dagger.add_callback!(proc -> begin
            return CuArrayDeviceProc(Distributed.myid(), #=CuContext(dev),=# dev.handle)
        end)
    end
end


text:
export  setpml!

function  setpml! end

_s(d::Real,L::Real,lam ; m::Real=1,s_min::Real = 3.0) = -im * s_min * (d^m*lam)/L^(m+1)
s(::High,x,d,l,lam)::Complex{Float} = (_s(x-(l - d),d,lam))
s(::Low,x,d,l,lam)::Complex{Float}  = (_s((-x + d) , d, lam))

function get_S_comp(sim::Simulation{Dim},  dir::Direction{D}) where {D,Dim}
		   if D == 1 return sim.S_x end
		   if D == 2 return sim.S_y end
	       if D == 3 return sim.S_z end
end

function (sim::Simulation)(::typeof(setpml!),dir::Direction{D},side::Low, thickness::Number,gridtype::GridType= p̂) where {D}
 S = get_S_comp(sim,dir)
 sim(add!,S,x -> s(side,x[D],thickness, extent(sim.grid,dir),sim.λ₀) ,x -> x[D] <= thickness, gridtype)
end

function (sim::Simulation)(::typeof(setpml!),dir::Direction{D},side::High, thickness::Number,gridtype::GridType= p̂) where {D}
S = get_S_comp(sim,dir)
 sim(add!,S,x -> s(side,x[D],thickness, extent(sim.grid,dir),sim.λ₀) ,x -> x[D] >= extent(sim.grid,dir) - thickness, gridtype)
end

function (sim::Simulation)(::typeof(setpml!),dir::Direction{D}, thickness::Number,gridtype::GridType= p̂) where {D}
 (sim::Simulation)(setpml!,dir,HIGH, thickness,gridtype)
 (sim::Simulation)(setpml!,dir,LOW, thickness,gridtype)
end

function (sim::Simulation{1})(::typeof(setpml!), thickness::Number,gridtype::GridType= p̂)
sim(setpml!,x̂, thickness,gridtype)
end
function (sim::Simulation{2})(::typeof(setpml!), thickness::Number,gridtype::GridType= p̂)
sim(setpml!,x̂, thickness,gridtype)
sim(setpml!,ŷ, thickness,gridtype)
end
function (sim::Simulation{3})(::typeof(setpml!), thickness::Number,gridtype::GridType= p̂)
sim(setpml!,x̂, thickness,gridtype)
sim(setpml!,ŷ, thickness,gridtype)
sim(setpml!,ẑ, thickness,gridtype)
end


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>15.0, :alg=>TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000)),
    "X0"   => Hyperrectangle(low=[0.6353, 14.7, 0.5573, 0.017, 0.0], high=[0.6353, 14.7, 0.5573, 0.017, 0.0]),
    "info" => Dict(
        "name" => "powertrain_control",
        "dim" => 5,
        "linear" => false,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

n = 10
plot(x=rand(n), y=rand(n), color=8.^rand(n),
     Scale.color_log2(minvalue=1, maxvalue=8))


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Wind

struct with inner constructor:
```julia
function Wind(;
    min_kw = 0.0,
    max_kw = 1.0e9,
    installed_cost_per_kw = 0.0,
    om_cost_per_kw = 40.0,
    prod_factor_series = missing,
    size_class = "",
    wind_meters_per_sec = [],
    wind_direction_degrees = [],
    temperature_celsius = [],
    pressure_atmospheres = [],
    macrs_option_years = 5,
    macrs_bonus_pct = 0.0,
    macrs_itc_reduction = 0.5,
    federal_itc_pct = 0.26,
    federal_rebate_per_kw = 0.0,
    state_ibi_pct = 0.0,
    state_ibi_max = 1.0e10,
    state_rebate_per_kw = 0.0,
    state_rebate_max = 1.0e10,
    utility_ibi_pct = 0.0,
    utility_ibi_max = 1.0e10,
    utility_rebate_per_kw = 0.0,
    utility_rebate_max = 1.0e10,
    production_incentive_per_kwh = 0.0,
    production_incentive_max_benefit = 1.0e9,
    production_incentive_years = 1,
    production_incentive_max_kw = 1.0e9,
    can_net_meter = true,
    can_wholesale = true,
    can_export_beyond_nem_limit = true
)
```

`size_class` must be one of ["residential", "commercial", "medium", "large"]. If `size_class` is not provided then it is determined based on the average electric load.

If no `installed_cost_per_kw` is provided (or it is 0.0) then it is determined from:
```julia
size_class_to_installed_cost = Dict(
    "residential"=> 11950.0,
    "commercial"=> 7390.0,
    "medium"=> 4440.0,
    "large"=> 3450.0
)
```

The Federal Investment Tax Credit is adjusted based on the `size_class` as follows (if the default of 0.3 is not changed):
```julia
size_class_to_itc_incentives = Dict(
    "residential"=> 0.3,
    "commercial"=> 0.3,
    "medium"=> 0.12,
    "large"=> 0.12
)
```

If the `prod_factor_series` is not provided then NREL's System Advisor Model (SAM) is used to get the wind turbine
production factor.

Wind resource values are optional, i.e.
(`wind_meters_per_sec`, `wind_direction_degrees`, `temperature_celsius`, and `pressure_atmospheres`).
If not provided then the resource values are downloaded from NREL's Wind Toolkit.
These values are passed to SAM to get the turbine production factor.

"""
struct Wind <: AbstractTech
    min_kw::Float64
    max_kw::Float64
    installed_cost_per_kw::Float64
    om_cost_per_kw::Float64
    prod_factor_series::Union{Missing, Array{Real,1}}
    size_class::String
    hub_height::T where T <: Real
    wind_meters_per_sec::AbstractArray{Float64,1}
    wind_direction_degrees::AbstractArray{Float64,1}
    temperature_celsius::AbstractArray{Float64,1}
    pressure_atmospheres::AbstractArray{Float64,1}
    macrs_option_years::Int
    macrs_bonus_pct::Float64
    macrs_itc_reduction::Float64
    federal_itc_pct::Float64
    federal_rebate_per_kw::Float64
    state_ibi_pct::Float64
    state_ibi_max::Float64
    state_rebate_per_kw::Float64
    state_rebate_max::Float64
    utility_ibi_pct::Float64
    utility_ibi_max::Float64
    utility_rebate_per_kw::Float64
    utility_rebate_max::Float64
    production_incentive_per_kwh::Float64
    production_incentive_max_benefit::Float64
    production_incentive_years::Int
    production_incentive_max_kw::Float64
    can_net_meter::Bool
    can_wholesale::Bool
    can_export_beyond_nem_limit::Bool
    can_curtail::Bool

    function Wind(;
        min_kw = 0.0,
        max_kw = 1.0e9,
        installed_cost_per_kw = 0.0,
        om_cost_per_kw = 40.0,
        prod_factor_series = missing,
        size_class = "",
        wind_meters_per_sec = [],
        wind_direction_degrees = [],
        temperature_celsius = [],
        pressure_atmospheres = [],
        macrs_option_years = 5,
        macrs_bonus_pct = 0.0,
        macrs_itc_reduction = 0.5,
        federal_itc_pct = 0.26,
        federal_rebate_per_kw = 0.0,
        state_ibi_pct = 0.0,
        state_ibi_max = 1.0e10,
        state_rebate_per_kw = 0.0,
        state_rebate_max = 1.0e10,
        utility_ibi_pct = 0.0,
        utility_ibi_max = 1.0e10,
        utility_rebate_per_kw = 0.0,
        utility_rebate_max = 1.0e10,
        production_incentive_per_kwh = 0.0,
        production_incentive_max_benefit = 1.0e9,
        production_incentive_years = 1,
        production_incentive_max_kw = 1.0e9,
        can_net_meter = true,
        can_wholesale = true,
        can_export_beyond_nem_limit = true,
        can_curtail= true,
        average_elec_load = 0.0
        )
        size_class_to_hub_height = Dict(
            "residential"=> 20,
            "commercial"=> 40,
            "medium"=> 60,  # Owen Roberts provided 50m for medium size_class, but Wind Toolkit has increments of 20m
            "large"=> 80
        )
        size_class_to_installed_cost = Dict(
            "residential"=> 11950.0,
            "commercial"=> 7390.0,
            "medium"=> 4440.0,
            "large"=> 3450.0
        )

        size_class_to_itc_incentives = Dict(
            "residential"=> 0.3,
            "commercial"=> 0.3,
            "medium"=> 0.12,
            "large"=> 0.12
        )

        if size_class == ""
            if average_elec_load <= 12.5
                size_class = "residential"
            elseif average_elec_load <= 100
                size_class = "commercial"
            elseif average_elec_load <= 1000
                size_class = "medium"
            else
                size_class = "large"
            end
        elseif !(size_class in keys(size_class_to_hub_height))
            @error "Wind.size_class must be one of $(keys(size_class_to_hub_height))"
        end

        if installed_cost_per_kw == 0.0
            installed_cost_per_kw = size_class_to_installed_cost[size_class]
        end

        if federal_itc_pct == 0.3
            federal_itc_pct = size_class_to_itc_incentives[size_class]
        end

        hub_height = size_class_to_hub_height[size_class]

        new(
            min_kw,
            max_kw,
            installed_cost_per_kw,
            om_cost_per_kw,
            prod_factor_series,
            size_class,
            hub_height,
            wind_meters_per_sec,
            wind_direction_degrees,
            temperature_celsius,
            pressure_atmospheres,
            macrs_option_years,
            macrs_bonus_pct,
            macrs_itc_reduction,
            federal_itc_pct,
            federal_rebate_per_kw,
            state_ibi_pct,
            state_ibi_max,
            state_rebate_per_kw,
            state_rebate_max,
            utility_ibi_pct,
            utility_ibi_max,
            utility_rebate_per_kw,
            utility_rebate_max,
            production_incentive_per_kwh,
            production_incentive_max_benefit,
            production_incentive_years,
            production_incentive_max_kw,
            can_net_meter,
            can_wholesale,
            can_export_beyond_nem_limit,
            can_curtail
        )
    end
end


text:
module Rethinking

using PrettyTables
using StatsBase
using Optim
using Turing
using LinearAlgebra
using DataFrames

include("quap.jl")
include("precise.jl")

export quap
export precise
end


text:
using Metatheory

@testset "Reduction Basics" begin
	t = @theory begin
	    ~a + ~a --> 2*(~a)
	    ~x / ~x --> 1
	    ~x * 1 --> ~x
	end

    # basic theory to check that everything works
    @test rewrite(:(a + a), t) == :(2a)
    @test rewrite(:(a + (x * 1)), t) == :(a + x)
	@test rewrite(:(a + (a * 1)), t; order=:inner) == :(2a)
end

import Base.(+)
@testset "Extending Algebra Operators" begin
    t = @theory begin
        ~a + ~a --> 2(~a)
    end

    # Let's extend an operator from base, for sake of example
    function +(x::Symbol, y)
        rewrite(:($x + $y), t)
    end

    @test (:x + :x) == :(2x)
end

## Free Monoid

@testset "Free Monoid - Overriding identity" begin
    # support symbol literals
	symbol_monoid = @theory begin
		~a ⋅ :ε --> ~a
		:ε ⋅ ~a --> ~a
		~a::Symbol --> ~a
		~a::Symbol ⋅ ~b::Symbol => Symbol(String(a) * String(b))
		# i |> error("unsupported ", i)
	end;

    @test rewrite(:(ε ⋅ a ⋅ ε ⋅ b ⋅ c ⋅ (ε ⋅ ε ⋅ d) ⋅ e), symbol_monoid; order=:inner) == :abcde
end

## Interpolation should be possible at runtime

@testset "Calculator" begin
	calculator = @theory begin
		~x::Number ⊕ ~y::Number => ~x + ~y
		~x::Number ⊗ ~y::Number => ~x * ~y
		~x::Number ⊖ ~y::Number => ~x ÷ ~y
		~x::Symbol --> ~x
		~x::Number --> ~x
	end;
	a = 10

	@test rewrite(:(3 ⊕ 1 ⊕ $a), calculator; order=:inner) == 14
end

## Direct rules
@testset "Direct Rules" begin
    t = @theory begin
        # maps
        ~a * ~b => ((~a isa Number && ~b isa Number) ? ~a * ~b : _lhs_expr)
    end
    @test rewrite(:(3 * 1), t) == 3

    t = @theory begin
        # maps
        ~a::Number * ~b::Number => ~a * ~b
    end
    @test rewrite(:(3 * 1), t) == 3
end

## Take advantage of subtyping.
# Subtyping in Julia has been formalized in this paper
# [Julia Subtyping: A Rational Reconstruction](https://benchung.github.io/papers/jlsub.pdf)

abstract type Vehicle end
abstract type GroundVehicle <: Vehicle end
abstract type AirVehicle <: Vehicle end
struct Airplane <: AirVehicle end
struct Car <: GroundVehicle end

airpl = Airplane()
car = Car()

t = @theory begin
	~a::AirVehicle * ~b => "flies"
	~a::GroundVehicle * ~b => "doesnt_fly"
end

@testset "Subtyping" begin

	sf = rewrite(:($airpl * c), t)
	df = rewrite(:($car * c), t)

    @test sf == "flies"
    @test df == "doesnt_fly"
end

@testset "Interpolation" begin
	airpl = Airplane()
	car = Car()
	t = @theory begin
		airpl * ~b => "flies"
		car * ~b => "doesnt_fly"
	end

	sf = rewrite(:($airpl * c), t)
	df = rewrite(:($car * c), t)

    @test sf == "flies"
    @test df == "doesnt_fly"
end

@testset "Segment Variables" begin
	t = @theory begin
		f(~x, ~~y) => Expr(:call, :ok, (~~y)...)
	end

	sf = rewrite(:(f(1,2,3,4)), t)

    @test sf == :(ok(2,3,4))

    t = @theory x y begin
		f(x, y...) => Expr(:call, :ok, y...)
	end

	sf = rewrite(:(f(1,2,3,4)), t)

    @test sf == :(ok(2,3,4))
end

module NonCall
using Metatheory
t = [@rule a b (a, b) --> ok(a,b)]

test() = rewrite(:(x,y), t)
end

@testset "Non-Call expressions" begin
	@test NonCall.test() == :(ok(x,y))
end

@testset "Pattern matcher can match on both function object references and name symbols" begin
	ex = :($(+)($(sin)(x)^2, $(cos)(x)^2))
	r = @rule(sin(~x)^2 + cos(~x)^2 --> 1)

	@test r(ex) == 1
end

@testset "Pattern variable as pattern term head" begin
	foo(x) = x+2
	ex = :(($foo)(bar, 2, pazz))
	r = @rule ((~f)(~x, 2, ~y) => (~f)(2))

	@test r(ex) == 4
end

using TermInterface

using Metatheory.Syntax: @capture
@testset "Capture form" begin
    ex = :(a^a)

    #note that @test inserts a soft local scope (try-catch) that would gobble
    #the matches from assignment statements in @capture macro, so we call it
    #outside the test macro
    ret = @capture ex (~x)^(~x)
    @test ret
    @test @isdefined x
    @test x === :a

    ex = :(b^a)
    ret = @capture ex (~y)^(~y)
    @test !ret
    @test !(@isdefined y)

    ret = @capture :(a + b) (+)(~~z)
    @test ret
    @test @isdefined z
    @test all(z .=== arguments(:(a + b)))

    #a more typical way to use the @capture macro

    f(x) = if @capture x (~w)^(~w)
        w
    end

    @test f(:(b^b)) == :b
    @test isnothing(f(:(b+b)))
end


text:
module QuantumESPRESSOBase

include("Inputs/Inputs.jl")

end # module


text:
# # Basic demonstration of Julia package composability

# Notebook from [HelloJulia.jl](https://github.com/ablaom/HelloJulia.jl)

# Instantiate package environment:

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# The Unitiful package allows you to **bind physical units** to numerical data:

using Unitful

#-

a = 5.0u"m^2/s^2"

#-

# The using Measurements package allows you to **propogate uncertainties**
# in numerical computations:

using Measurements

#

b = 5.0 ± 1.2 # or measurement(5.0, 1.2)

#-

sqrt(b)

#-

# The Unitful and Measurements packages are blissfully ignorant of one
# another. That, is neither package is a dependency of the other. And yet the
# following "just works", as if by magic:
c = 5.0u"m^2/s^2" ± 0.1u"m^2/s^2"

#-

sqrt(c)


text:
struct DPM_MCMC_Chain
    gamma::Vector{Vector{Bool}}
    f::Vector{Vector{Float64}}
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DedicatedHostProperties <: SwaggerModel
    platformFaultDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: platformFaultDomain
    autoReplaceOnFailure::Any # spec type: Union{ Nothing, Bool } # spec name: autoReplaceOnFailure
    hostId::Any # spec type: Union{ Nothing, String } # spec name: hostId
    virtualMachines::Any # spec type: Union{ Nothing, Vector{SubResourceReadOnly} } # spec name: virtualMachines
    licenseType::Any # spec type: Union{ Nothing, DedicatedHostLicenseType } # spec name: licenseType
    provisioningTime::Any # spec type: Union{ Nothing, DateTime } # spec name: provisioningTime
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    instanceView::Any # spec type: Union{ Nothing, DedicatedHostInstanceView } # spec name: instanceView

    function DedicatedHostProperties(;platformFaultDomain=nothing, autoReplaceOnFailure=nothing, hostId=nothing, virtualMachines=nothing, licenseType=nothing, provisioningTime=nothing, provisioningState=nothing, instanceView=nothing)
        o = new()
        validate_property(DedicatedHostProperties, Symbol("platformFaultDomain"), platformFaultDomain)
        setfield!(o, Symbol("platformFaultDomain"), platformFaultDomain)
        validate_property(DedicatedHostProperties, Symbol("autoReplaceOnFailure"), autoReplaceOnFailure)
        setfield!(o, Symbol("autoReplaceOnFailure"), autoReplaceOnFailure)
        validate_property(DedicatedHostProperties, Symbol("hostId"), hostId)
        setfield!(o, Symbol("hostId"), hostId)
        validate_property(DedicatedHostProperties, Symbol("virtualMachines"), virtualMachines)
        setfield!(o, Symbol("virtualMachines"), virtualMachines)
        validate_property(DedicatedHostProperties, Symbol("licenseType"), licenseType)
        setfield!(o, Symbol("licenseType"), licenseType)
        validate_property(DedicatedHostProperties, Symbol("provisioningTime"), provisioningTime)
        setfield!(o, Symbol("provisioningTime"), provisioningTime)
        validate_property(DedicatedHostProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(DedicatedHostProperties, Symbol("instanceView"), instanceView)
        setfield!(o, Symbol("instanceView"), instanceView)
        o
    end
end # type DedicatedHostProperties

const _property_map_DedicatedHostProperties = Dict{Symbol,Symbol}(Symbol("platformFaultDomain")=>Symbol("platformFaultDomain"), Symbol("autoReplaceOnFailure")=>Symbol("autoReplaceOnFailure"), Symbol("hostId")=>Symbol("hostId"), Symbol("virtualMachines")=>Symbol("virtualMachines"), Symbol("licenseType")=>Symbol("licenseType"), Symbol("provisioningTime")=>Symbol("provisioningTime"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("instanceView")=>Symbol("instanceView"))
const _property_types_DedicatedHostProperties = Dict{Symbol,String}(Symbol("platformFaultDomain")=>"Int32", Symbol("autoReplaceOnFailure")=>"Bool", Symbol("hostId")=>"String", Symbol("virtualMachines")=>"Vector{SubResourceReadOnly}", Symbol("licenseType")=>"DedicatedHostLicenseType", Symbol("provisioningTime")=>"DateTime", Symbol("provisioningState")=>"String", Symbol("instanceView")=>"DedicatedHostInstanceView")
Base.propertynames(::Type{ DedicatedHostProperties }) = collect(keys(_property_map_DedicatedHostProperties))
Swagger.property_type(::Type{ DedicatedHostProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DedicatedHostProperties[name]))}
Swagger.field_name(::Type{ DedicatedHostProperties }, property_name::Symbol) =  _property_map_DedicatedHostProperties[property_name]

function check_required(o::DedicatedHostProperties)
    true
end

function validate_property(::Type{ DedicatedHostProperties }, name::Symbol, val)
    if name === Symbol("platformFaultDomain")
        Swagger.validate_param(name, "DedicatedHostProperties", :maximum, val, 2, false)
        Swagger.validate_param(name, "DedicatedHostProperties", :minimum, val, 0, false)
    end
end


text:
"""
    apply_axis!(g, a)

Internal function to apply an `Axis` object `a` in a GLE context.
"""
function apply_axis!(g::GLE, a::Axis, parent_font::String)
    parent_font = ifelse(isdef(a.textstyle.font), a.textstyle.font, parent_font)
    apply_ticks!(g, a.ticks, a.prefix, parent_font)
    if isdef(a.title)
        apply_title!(g, a.title, a.prefix, parent_font)
    end
    # XXX subticks disabled for now
    "\n\t$(a.prefix)subticks off" |> g
    #
    "\n\t$(a.prefix)axis" |> g
    a.off && ("off" |> g; return nothing)
    a.log && "log"  |> g
    isdef(a.base)   && "base $(a.base)"     |> g
    isdef(a.lwidth) && "lwidth $(a.lwidth)" |> g
    isdef(a.min)    && "min $(a.min)"       |> g
    isdef(a.max)    && "max $(a.max)"       |> g
    a.ticks.grid    && "grid"               |> g
    apply_textstyle!(g, a.textstyle, parent_font)
    return nothing
end

"""
    apply_axes!(g, a, figid)

Internal function to apply an `Axes2D` object `a` in a GLE context.
The `figid` is useful to keep track of the figure the axes belong to
which is required in the `apply_drawings` subroutine that is called.
"""
function apply_axes!(g::GLE, a::Axes2D, figid::String, axidx::Int)
    a.off && return nothing

    isdef(a.origin) && "\namove $(a.origin[1]) $(a.origin[2])" |> g
    if a.scale != ""
        scale = ifelse(isdef(a.origin), "fullsize", "scale $(a.scale)")
    end

    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    "\nbegin graph\n\t$scale"   |> g

    # graph >> math mode (crossing axis)
    a.math && "\n\tmath" |> g
    # -- size of the axes, see also layout
    isdef(a.size) && "\n\tsize $(a.size[1]) $(a.size[2])" |> g

    # graph >> apply axis (ticks, ...), passing the figure font as parent font (see issue #76)
    parent_font = Figure(figid; _noreset=true).textstyle.font
    for axis in (a.xaxis, a.x2axis, a.yaxis, a.y2axis)
        apply_axis!(g, axis, parent_font)
    end

    # graph >> apply axes title, passing the figure font as parent font
    isdef(a.title) && apply_title!(g, a.title, "", parent_font)

    # graph >> apply drawings
    apply_drawings!(g, a.drawings, figid, axidx)

    "\nend graph" |> g
    # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    # apply  legend and other floating objects
    isdef(a.legend)    && apply_legend!(g, a.legend, parent_font, figid)
    isempty(a.objects) || apply_objects!(g, a.objects, figid)
    return nothing
end

function apply_axes!(g::GLE, a::Axes3D, figid::String, axidx::Int)
# TODO
# -- title
# -- perspective

    a.off && return nothing
    axid = "a3d_$(hash(a))"

    #
    # begin object ax3d_hash
    #   begin surface
    #       size x y
    #       cube xlen 10 ylen 10 zlen 10 lstyle 9 color blue
    #       xaxis ...
    #       xtitle ...
    #       yaxis ...
    #       ytitle ...
    #       zaxis ...
    #       ztitle ...
    #       (surface)
    #   end surface
    #   (objects)
    # end object
    #
    # amove (appropriate location)
    # draw ax3d_hash.cc
    #

    "\nbegin object $axid" |> g
    "\n\tbegin surface"    |> g
    "\n\t\tsize $(a.size[1]) $(a.size[2])" |> g
    # ------------------------------------------
    # CUBE
    "\n\t\tcube"      |> g
    a.nocube && "off" |> g
    "xlen $(a.cubedims[1]) ylen $(a.cubedims[2]) zlen $(a.cubedims[3])"   |> g
    a.nocube || apply_linestyle!(g, a.linestyle)

# TODO should become apply_axis3 or something
    δx = round((a.xaxis.max - a.xaxis.min)/5, digits=1)
    "\n\t\txaxis min $(a.xaxis.min) max $(a.xaxis.max) dticks $δx" |> g
    δy = round((a.yaxis.max - a.yaxis.min)/5, digits=1)
    "\n\t\tyaxis min $(a.yaxis.min) max $(a.yaxis.max) dticks $δy" |> g
    δz = round((a.zaxis.max - a.zaxis.min)/5, digits=1)
    "\n\t\tzaxis min $(a.zaxis.min) max $(a.zaxis.max) dticks $δz" |> g

    # ROTATION
    if isdef(a.rotate)
        "\n\t\trotate $(a.rotation[1]) $(a.rotation[2]) 0" |> g
    else
        "\n\t\trotate 65 20 0" |> g
    end

    # XXX AXIS
    # parent_font = Figure(figid; _noreset=true).textstyle.font
    # for axis in (a.xaxis, a.yaxis, a.zaxis)
    #     apply_axis!(g, axis, parent_font)
    # end

    # SURFACE
    if isempty(a.drawings) || all(d->!isa(d, Surface), a.drawings)
        # NOTE if there is no surface, we MUST add dummy data otherwise ghostscript crashes.
        fd = joinpath(GP_ENV["TMP_PATH"], "$(figid)_dummy.z")
        write(fd, "! nx 2 ny 2 xmin 1 xmax 2 ymin 1 ymax 2\n1 2\n2 2\n")
        "\n\t\tdata \"$fd\""   |> g
        "\n\t\ttop off"        |> g
        "\n\t\tunderneath off" |> g
    end
    surfs = [i for i ∈ 1:length(a.drawings) if a.drawings[i] isa Surface]
    apply_drawings!(g, a.drawings[surfs], figid, axidx)
    # -----------------------------------------
    "\n\tend surface"      |> g
    apply_drawings!(g, a.drawings[setdiff(1:length(a.drawings), surfs)], figid, axidx)
    # OBJECTS
    apply_objects!(g, a.objects, figid)

    "\nend object"         |> g

    if isdef(a.origin)
        # move to center of container
        cx = a.origin[1] + a.size[1]/2
        cy = a.origin[2] + a.size[2]/2
        "\namove $cx $cy" |> g
    else
        # move to center of page
        "\namove pagewidth()/2 pageheight()/2" |> g
    end
    # draw the overall container centered
    "\ndraw $axid.cc" |> g
end


text:
include("functions-for-periodic-vectors.jl")
using LinearAlgebra, Statistics, SpecialFunctions

const default_tol_force_equilibrium = 1e-12

"Define an abstract type for all the type of particles that will be used."
abstract type AbstractParticle{d, T} end # To define an abstract type for all type of particles to be used

#########################################################################################################
#########################################################################################################
### 'MonoParticle' type and related constructors and functions. Type meant to be used in monodisperse packings
#########################################################################################################
#########################################################################################################

#=
NB:We are including the docstring of all the fields of 'MonoParticle', following the recommendations in https://docs.julialang.org/en/v1/manual/documentation/ . However, such usage produces an error in parametric struct's, just as 'MonoParticle' is (see, e.g. https://github.com/JuliaLang/julia/issues/39410). So we are leaving such docstring, but also providing the description of each field in the docstring of the struct itself.
=#

"""
    MonoParticle{d, T}

Struct used to construct a particle in a d-dimensional space.

It's meant to be used as components of *monodisperse* packings. Hence, while many of its
relevant properties are stored, its radius is instead stored in a field in 'MonoPacking' type.

The quantities this struct has are:
1. X: The particle's center (as a StaticVector{d, PeriodicNumber{T}}})
2. contact_vecs: The set of contact vectors with the particle's neighbours (as a Vector containing SVector{d,T})
3. forces: The magnitude of contact forces (Vector{T})
4. neighbours: The indices of the particle's neighbours (as Vector{Int64})

Naturally, 'contact_vecs', 'forces', and 'neighbours' have the same length.

See also: [`MonoPacking`](@ref).
"""
mutable struct MonoParticle{d, T}  <: AbstractParticle{d, T} # a struct defining a particle in a d-dimensional MONODISPERSE packing
    "Position of particle's center (SVector{d, PeriodicNumber{T}}"
    X::SVector{d, PeriodicNumber{T}}   # particle's centre
    "Vector of contact vectors, pointing towards X. Array of {SVector{d,T}}"
    contact_vecs::Vector{SVector{d,T}} # array of contact vectors (should point towards X)
    "Vector with the magnitudes of the contact forces"
    forces::Vector{T}  # list of forces magnitudes
    "List of particles' indices with which  the particle is in contact"
    neighbours::Vector{Int64}  # list of the indices, corresponding to particles in contact
end

#### Format how a `MonoParticle` is shown in IO
function Base.show(io::IO, particle::MonoParticle{d,T}) where {d, T<:Real}
    x = value.(particle.X)
    neighs = particle.neighbours
    forces = particle.forces

    println(io, d, "d MonoParticle\t centred at:  ", x)
    println(io, "neighbours: ", neighs)
    println(io, "contact forces: ", forces)
end

### Now we define some convenient constructors (and related functions) for 'MonoParticle'

#This function is suggested in the 'StaticArrays' documentation (https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#Guide-1)
"Convert a matrix to a vector with SVector elements."
function svectors(M::Matrix{T}, ::Val{d}) where {T,d}
    size(M,1) == d || error("sizes mismatch")
    isbitstype(T) || error("use for bitstypes only")
    copy(reinterpret(SVector{d,T}, vec(M)))
end
svectors(reshape(collect(1:6), (2,3)), Val(2))

function MonoParticle(X::SVector{d, PeriodicNumber{T}}, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {d, T<:Real}

    dim, z = size(contact_vecs)
    (dim!=d) && error("There is a mismatch in the dimensionality!! Dimensions of particle's position: ", d, "\t dimensions of contact vectors: ", dim)

    if z==length(fs) && z==length(neighbours)
        return MonoParticle(X, svectors(contact_vecs, Val(d)), fs, neighbours)
    else
        error("Number of contacts mismatch! There are: ", z, " contact vectors \t", length(fs), " forces and \t ", length(neighbours), " contacts indices")
    end
end
MonoParticle(SVector{3}(PeriodicNumber.(rand(3), 0.5)), rand(3,10), rand(10), rand(1:100,10))

# The following two methods are useful if one's being forgetful (or a bit lazy) about defining 'X' as a SVector of 'PeriodicNumber' elements

function MonoParticle(X::Vector{PeriodicNumber{T}}, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {T<:Real}
    d = length(X)
    MonoParticle(SVector{d}(X), contact_vecs, fs, neighbours)
end
MonoParticle(PeriodicNumber.(rand(3), 0.5), rand(3,10), rand(10), rand(1:100,10))

MonoParticle(X::Vector{T}, L::T, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {T<:Real} = MonoParticle(PeriodicNumber.(X,L), contact_vecs, fs, neighbours)

#########################################################################################################
#########################################################################################################
### 'Particle' type and related constructors and functions. Type meant to be used in polydisperse packings
#########################################################################################################
#########################################################################################################

"""
    Particle{d, T}
    This struct is used to construct a particle in a d-dimensional space. It's meant to be used as components of *bidisperse* or *polydisperse* packings. Hence, in this case also its radius is stored, along with other of its relevant properties are stored.
    The quantities this struct has are:
        1. X: The particle's center (as a StaticVector{d, PeriodicNumber{T}}})
        2. R: The particle's radius
        2. contact_vecs: The set of contact vectors with the particle's neighbours (as a Vector containing SVector{d,T})
        3. forces: The magnitude of contact forces (Vector{T})
        4. neighbours: The indices of the particle's neighbours (as Vector{Int64})
    Naturally, 'contact_vecs', 'forces', and 'neighbours' have the same length.
"""
mutable struct Particle{d, T}  <: AbstractParticle{d, T} # a struct defining a particle in a d-dimensional MONODISPERSE packing
    X::SVector{d, PeriodicNumber{T}}   # particle's centre
    R::T # particle's radius
    contact_vecs::Vector{SVector{d,T}} # array of contact vectors (should point towards X)
    forces::Vector{T}  # list of forces magnitudes
    neighbours::Vector{Int64}  # list of the indices, corresponding to particles in contact
end

#### Format how a `Particle` is shown in IO
function Base.show(io::IO, particle::Particle{d,T}) where {d, T<:Real}
    x = value.(particle.X)
    neighs = particle.neighbours
    forces = particle.forces

    println(io, d, "d Particle\t centred at:  ", x, " \t and radius R = ", particle.R)
    println(io, "neighbours: ", neighs)
    println(io, "contact forces: ", forces)
end

function Particle(X::SVector{d, PeriodicNumber{T}}, R::T, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {d, T<:Real}

    z = size(contact_vecs, 2)

    if z==length(fs) && z==length(neighbours)
        Scvs = Vector{SVector{d,T}}(undef,z)
        for k in 1:z
            Scvs[k] = SVector{d}(contact_vecs[:, k])
        end

        return Particle(X, R, Scvs, fs, neighbours)
    else
        error("Number of contacts mismatch! There are: ", z, " contact vectors \t", length(fs), " forces and \t ", length(neighbours), " contacts indices")
    end
end

# The following two methods are useful if one's being forgetful (or a bit lazy) about defining 'X' as a SVector of 'PeriodicNumber' elements
function Particle(X::Vector{PeriodicNumber{T}}, R::T, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {T<:Real}

    d, z = size(contact_vecs)

    if z==length(fs) && z==length(neighbours)
        SX = SVector{d}(X)

        Scvs = Vector{SVector{d,T}}(undef,z)
        for k in 1:z
            Scvs[k] = SVector{d}(contact_vecs[:, k])
        end

        return Particle(SX, R, Scvs, fs, neighbours)
    else
        error("Number of contacts mismatch! There are: ", z, " contact vectors \t", length(fs), " forces and \t ", length(neighbours), " contacts indices")
    end
end
Particle(PeriodicNumber.(rand(3), 0.5), rand(), rand(3,10), rand(10), rand(1:100,10))

Particle(X::Vector{T}, R::T, L::T, contact_vecs::Matrix{T}, fs::Vector{T}, neighbours::Vector{Int64}) where {T<:Real} = Particle(PeriodicNumber.(X,L), R, contact_vecs, fs, neighbours)

#########################################################################################################
#########################################################################################################
### Other functions; e.g. to compute volume
#########################################################################################################
#########################################################################################################

"""
    volume_d_ball(d::Int64, R::Real)

Compute the volume of a d-dimensional sphere of radius R (that defaults to 1)
"""
function volume_d_ball(d::Int64, R::Real=1)::Float64
    ((sqrt(π)*R)^d)/gamma(1.0 + 0.5d)
end

"""
    volume(P::Particle)

Compute the volume of 'P' (i.e. a hypersphere of d dimensions and radius R=P.R.
"""
volume(P::Particle) = volume_d_ball(length(P.X), P.R)

function check_for_overlaps(P1::Particle{d,T}, P2::Particle{d,T}, tolerance::T) where {d, T<:AbstractFloat}
    σij = P1.R + P2.R
    dist = norm(P1.X - P2.X)
    overlap=false; gap = 0.0
    if dist-σij < -tolerance
        gap = σji - dist
        overlap = true
    end
    return overlap, gap
end

# function check_for_overlaps(Xs::Vector{SVector{d, PeriodicNumber{T}}}, R::T, tolerance::T) where {d, T<:AbstractFloat}
#     σ = 2*R
#     N = length(Xs)
#     distances = distances_between_centers(Xs)
#     overlap=false; message = "No overlap is present"; particles = (0, 0)
#     for i in 1:N, j in i+1:N
#         if distances[i,j]-σ < -tolerance
#             gap = σ - distances[i,j]
#             Xi = value.(Xs[i])
#             Xj = value.(Xs[j])
#             message = "Overlap between particles $i and $j;\n Centers at: $Xi; \t $Xj. \n\t\tradius:  $R; \t Overlap = $gap ."

#             overlap = true
#             particles = (i, j)
#             break
#         end
#     end
#     return overlap, message, particles
# end

#########################################################################################################
#########################################################################################################
### Functions for generic (i.e. abstract) type of particle
#########################################################################################################
#########################################################################################################

"Return the coordination number (z) of 'P'."
get_coordination_number(P::AbstractParticle) = length(P.neighbours)
# get_coordination_number(P1); get_coordination_number(P3);

"Return the norm of each of the contact vectors of 'P'."
norm_contact_vectors(P::AbstractParticle) = norm.(P.contact_vecs)
# norm_contact_vectors(P1); norm_contact_vectors(P3);

"""
Return whether 'P' is a rattler or not; that is, if it has less than d+1 contacts.

See also: [`is_stable_particle`](@ref).
"""
function is_a_rattler(P::AbstractParticle{d, T}) where {d, T<:Real}
    z = get_coordination_number(P)
    if z<=d
        return true
    else
        return false
    end
end

"""
Return whether 'P' is a stable particle or not; that is, if it has more than d contacts.

See also: [`is_a_rattler`](@ref).
"""
is_stable_particle(P::AbstractParticle{d, T}) where {d, T<:Real} = !is_a_rattler(P)
is_stable_particle(Particle(rand(3), rand(), 0.4, rand(3,3), rand(3), collect(1:3)));

"Compute the sum of forces acting on a given particle. The output is a `StaticVector`"
function total_force(P::AbstractParticle{d, T}) where {d, T<:Real}
    fs = P.forces
    ctc_vecs = P.contact_vecs
    # tot_f = @SVector zeros(size(P.X, 1))
    tot_f = @SVector zeros(d)
    for (k, f) in enumerate(fs)
        tot_f += f*ctc_vecs[k]
    end
    return tot_f
end

"""
Test whether the sum of forces in a set of particles is smaller than a given tolerance value.

By default such tolerance is `default_tol_force_equilibrium` (a `const`) and equals 1e-12.
"""
function force_equilibrium(particles::Vector{<:AbstractParticle} ; tol_mechanical_equilibrium::Float64=default_tol_force_equilibrium)
    all(x-> x <= tol_mechanical_equilibrium, norm.(total_force.(particles)))
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
#include("../datasets/dataset_reader.jl")
#include("../datasets/Demo/datasetDemo.jl")
#include("../datasets/Demo-overlap/datasetDemoOverlap.jl")
#include("../datasets/Drug-Classification/dataset1.jl")
#include("../datasets/Income-Classification/dataset2.jl")
#include("../datasets/Wine-Classification/dataset3.jl")
#include("../datasets/Drug-Consumism-Classification/dataset4.jl")
#include("../datasets/Heart-Disease-Classification/dataset5.jl")

const MAX_DATASET_ID = 7
const DEFAULT_DATASET_PATH = "../datasets"

function f1!(ch::Chromosome, dt_setup::DatasetSetup)::Float64
    ch.fitness = 1.0
end

function f2!(ch::Chromosome, dt_setup::DatasetSetup)::Float64
    mlp = convert_to_FluxNet(ch, dt_setup.labels)
    error = error_rate(zip(dt_setup.input_train, dt_setup.output_train), mlp)
    ch.fitness = (2.0/(error+1.0)) - 1.0
end

fitness_functions = Function[f1!, f2!]
dataset_functions = Function[get_dataset0_setup, get_dataset0_overlap_setup, get_dataset1_setup,
                             get_dataset2_setup, get_dataset3_setup, get_dataset4_setup,
                             get_dataset5_setup]
dataset_paths = String["Demo", "Demo-overlap", "Drug-Classification", "Income-Classification",
                       "Wine-Classification", "Drug-Consumism-Classification",
                       "Heart-Disease-Classification"]

function get_dataset(dataset_id::Int=0; normalise::Bool=false, balance::Bool=false)
    @assert length(dataset_functions) == MAX_DATASET_ID
    @assert length(dataset_paths) == MAX_DATASET_ID
    if dataset_id > MAX_DATASET_ID || dataset_id <= 0
        get_setup = dataset_functions[1]
        path = joinpath(DEFAULT_DATASET_PATH, dataset_paths[1])
    else
        get_setup = dataset_functions[dataset_id+1]
        path = joinpath(DEFAULT_DATASET_PATH, dataset_paths[dataset_id+1])
    end
    setup = get_setup(path, normalise=normalise, balance=balance, separate=true)
    return setup
end

struct Evaluation <: Function
    dataset_setup::DatasetSetup
    fitness_function::Function
    function Evaluation(dataset_id::Int=0, function_id::Int=1; normalise::Bool=false,
                        balance::Bool=false)
        dataset_setup = get_dataset(dataset_id, normalise=normalise, balance=balance)
        if function_id < 1 || function_id > length(fitness_functions)
            fitness = fitness_functions[1]
        else
            fitness = fitness_functions[function_id]
        end
        new(dataset_setup, fitness)
    end
end

(e::Evaluation)(ch::Chromosome) = e.fitness_function(ch, e.dataset_setup)

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT