text:
const APMorVM = Union{AbstractPeriodicArray, AbstractVecOrMat}
"""
    ps([PMT::Type,] A::PM1, B::PM2, C::PM3, D::PM4) -> psys::PeriodicStateSpace

Construct a `PeriodicStateSpace` object from a quadruple of periodic matrix objects.

The periodic matrix objects `A`, `B`, `C`, `D` specifies the periodic matrices
`A(t)`, `B(t)`, `C(t)` and `D(t)` of a linear periodic time-varying state space model
in the continuous-time form

     dx(t)/dt = A(t)x(t) + B(t)u(t) ,
     y(t)     = C(t)x(t) + D(t)u(t) ,

or in the discrete-time form

     x(t+1)  = A(t)x(t) + B(t)u(t) ,
     y(t)    = C(t)x(t) + D(t)u(t) ,

where `x(t)`, `u(t)` and `y(t)` are the system state vector,
system input vector and system output vector, respectively,
and `t` is the continuous or discrete time variable.
The system matrices satisfy `A(t) = A(t+T₁)`, `B(t) = B(t+T₂)`, `C(t) = C(t+T₃)`, `D(t) = D(t+T₄)`,
i.e., are periodic with periods `T₁`, `T₂`, `T₃` and `T₄`, respectively.
The different periods must be commensurate (i.e., their ratios must be rational numbers with
numerators and denominators up to at most 4 decimal digits).
The periodic matrix objects `A`, `B`, `C`, `D` can have different types `PM1`, `PM2`, `PM3`, `PM4`, respectively,
where for a contiuous-time system
`PM1`, `PM2`, `PM3`, `PM4` must be one of the supported continuous-time periodic matrix types, i.e.,
[`PeriodicFunctionMatrix`](@ref), [`PeriodicSymbolicMatrix`](@ref),
[`HarmonicArray`](@ref), [`FourierFunctionMatrix`](@ref) or [`PeriodicTimeSeriesMatrix`](@ref),
while for a discrete-time system  `PM1`, `PM2`, `PM3`, `PM4` must be one of the supported discrete-time periodic matrix types, i.e.,
[`PeriodicMatrix`](@ref) or [`PeriodicArray`](@ref).
Any of the objects `A`, `B`, `C`, `D` can be also specified as a real matrix or vector of appropriate size.

If `PMT` is not specified, the resulting `psys` has periodic matrices of the same type `PT`, such that `PT` is either the common type
of all matrix objects or `PT = PeriodicFunctionMatrix` for a continuous-time system or  `PT = PeriodicMatrix` for a discrete-time system.
If `PMT` is specified, the resulting `psys` has periodic matrices of the same type `PMT`.

Other convenience constructors are implemented as follows:

    ps([PMT::Type,] A, B, C) -> psys::PeriodicStateSpace

to construct a `PeriodicStateSpace` object for a quadruple of the form `(A,B,C,0)`;

    ps(D) -> psys::PeriodicStateSpace

to construct a `PeriodicStateSpace` object for a quadruple of the form `([],[],[],D)`;

    ps([PMT::Type,] A, B, C, D, period) -> psys::PeriodicStateSpace

to construct a `PeriodicStateSpace` object with a desired `period` for a quadruple `(A,B,C,D)`;
all objects `A`, `B`, `C`, `D` can be also specified as real matrices or vectors of appropriate sizes.

    ps([PMT::Type,] A, B, C, period) -> psys::PeriodicStateSpace

to construct a `PeriodicStateSpace` object with a desired `period` for a quadruple `(A,B,C,0)`;
all objects `A`, `B`, `C` can be also specified as real matrices or vectors of appropriate sizes.

    ps(sys::DescriptorStateSpace, period) -> psys::PeriodicStateSpace

to construct a `PeriodicStateSpace` object with a desired `period` for a quadruple `(sys.A,sys.B,sys.C,sys.D)`
(`sys.E = I` is assumed).
"""
function ps(A::APMorVM, B::APMorVM, C::APMorVM, D::APMorVM)
   PMT = (typeof(A),typeof(B),typeof(C),typeof(D))
   PMTM = trues(4)
   cont = nothing
   for i = 1:4
       PMT[i] <: AbstractVecOrMat && continue
       PMTM[i] = false
       t = PMT[i].parameters[1] == :c
       if isnothing(cont)
          cont = t
       else
          cont == t ||  error("all matrix objects must be either continuous or discrete")
       end
   end
   all(PMTM) && error("period and sample time must be specified in the case of only matrix inputs")
   period = promote_period(A,B,C,D)
   PMTS = cont ? promote_cpmtype(A,B,C,D) : promote_dpmtype(A,B,C,D)
   return PeriodicStateSpace(PMTM[1] ? PMTS(A,period) : (PMT[1] <: PMTS ? A : convert(PMTS,A)),
                             PMTM[2] ? PMTS(B,period) : (PMT[2] <: PMTS ? B : convert(PMTS,B)),
                             PMTM[3] ? PMTS(C,period) : (PMT[3] <: PMTS ? C : convert(PMTS,C)),
                             PMTM[4] ? PMTS(D,period) : (PMT[4] <: PMTS ? D : convert(PMTS,D)))

end
function ps(PMT::Type, A::APMorVM, B::APMorVM, C::APMorVM, D::APMorVM)
    period = promote_period(A,B,C,D)
    ps(typeof(A) <: AbstractVecOrMat ? PMT(A,period) : convert(PMT,A),
       typeof(B) <: AbstractVecOrMat ? PMT(B,period) : convert(PMT,B),
       typeof(C) <: AbstractVecOrMat ? PMT(C,period) : convert(PMT,C),
       typeof(D) <: AbstractVecOrMat ? PMT(D,period) : convert(PMT,D))
end
ps(A::APMorVM, B::APMorVM, C::APMorVM, D::APMorVM, period::Real) =
    ps(set_period(A,period), set_period(B,period), set_period(C,period), set_period(D,period))
ps(PMT::Type, A::APMorVM, B::APMorVM, C::APMorVM, D::APMorVM, period::Real) =
    ps(typeof(A) <: AbstractVecOrMat ? A : convert(PMT,set_period(A,period)),
       typeof(B) <: AbstractVecOrMat ? B : convert(PMT,set_period(B,period)),
       typeof(C) <: AbstractVecOrMat ? C : convert(PMT,set_period(C,period)),
       typeof(D) <: AbstractVecOrMat ? D : convert(PMT,set_period(D,period)))

ps(A::APMorVM, B::APMorVM, C::APMorVM) = ps(A,B,C,zeros(Float64,size(C,1)[1],size(B,2)[1]))
function ps(PMT::Type, A::APMorVM, B::APMorVM, C::APMorVM)
   period = promote_period(A,B,C)
   ps(typeof(A) <: AbstractVecOrMat ? PMT(A,period) : convert(PMT,A),
      typeof(B) <: AbstractVecOrMat ? PMT(B,period) : convert(PMT,B),
      typeof(C) <: AbstractVecOrMat ? PMT(C,period) : convert(PMT,C),
      PMT(zeros(Float64,size(C,1)[1],size(B,2)[1]),period))
end
ps(A::APMorVM, B::APMorVM, C::APMorVM, period::Real) =
    ps(set_period(A,period), set_period(B,period), set_period(C,period), zeros(size(C,1),size(B,2)))
ps(PMT::Type, A::APMorVM, B::APMorVM, C::APMorVM, period::Real) =
    ps(PMT, set_period(A,period), set_period(B,period), set_period(C,period))

function ps(D::PM) where {PM <: AbstractPeriodicArray}
    p, m = size(D,1), size(D,2)
    PMT = typeof(D).name.wrapper
    ps(convert(PMT,PeriodicFunctionMatrix(zeros(0,0),D.period)),
       convert(PMT,PeriodicFunctionMatrix(zeros(0,m),D.period)),
       convert(PMT,PeriodicFunctionMatrix(zeros(p,0),D.period)), D)
end
function ps(sys::DST, period::Real) where {DST <: DescriptorStateSpace}
    sys.E == I || error("only standard state-spece models supported")
    Ts = sys.Ts
    if Ts == 0
       ps(PeriodicFunctionMatrix(sys.A,period), PeriodicFunctionMatrix(sys.B,period),
          PeriodicFunctionMatrix(sys.C,period),PeriodicFunctionMatrix(sys.D,period))
    else
      r = rationalize(period/abs(Ts))
      denominator(r) == 1 || error("incommensurate period and sample time")
      d = Ts == -1 ? 1 : numerator(r)
      ps(PeriodicMatrix(sys.A, period; nperiod = d), PeriodicMatrix(sys.B,period; nperiod = d),
         PeriodicMatrix(sys.C,period; nperiod = d), PeriodicMatrix(sys.D,period; nperiod = d))
    end
end
islti(psys::PeriodicStateSpace) = isconstant(psys.A) && isconstant(psys.B) && isconstant(psys.C) && isconstant(psys.D)

function ps(A::AbstractVecOrMat, B::AbstractVecOrMat, C::AbstractVecOrMat, D::AbstractVecOrMat, period::Real; Ts::Union{Real,Missing} = missing)
    if ismissing(Ts)
       ps(PeriodicFunctionMatrix(A,period), PeriodicFunctionMatrix(B,period),
          PeriodicFunctionMatrix(C,period), PeriodicFunctionMatrix(D,period))
    else
       Ts == -1 || Ts > 0 || error("sampling time must be -1 or positive")
       r = rationalize(period/abs(Ts))
       denominator(r) == 1 || error("incommensurate period and sample time")
       d = numerator(r)
       ps(PeriodicMatrix(A, period; nperiod = d), PeriodicMatrix(B,period; nperiod = d),
          PeriodicMatrix(C,period; nperiod = d), PeriodicMatrix(D,period; nperiod = d))
    end
end
ps(A::AbstractVecOrMat, B::AbstractVecOrMat, C::AbstractVecOrMat, period::Real; Ts::Union{Real,Missing} = missing) =
   ps(A, B, C, zeros(size(C,1),size(B,2)), period; Ts)
isct(psys::PeriodicStateSpace) = iscontinuous(psys.A)
isdt(psys::PeriodicStateSpace) = !iscontinuous(psys.A)


text:
"""
A pointer to a `T` stored inside a Blob.
"""
struct Blob{T}
    base::Ptr{Nothing}
    offset::Int64
    limit::Int64

    function Blob{T}(base::Ptr{Nothing}, offset::Int64, limit::Int64) where {T}
        @assert isbitstype(T)
        new(base, offset, limit)
    end
end

function Blob{T}(blob::Blob) where T
    Blob{T}(getfield(blob, :base), getfield(blob, :offset), getfield(blob, :limit))
end

function assert_same_allocation(blob1::Blob, blob2::Blob)
    @assert getfield(blob1, :base) == getfield(blob2, :base) "These blobs do not share the same allocation: $blob1 - $blob2"
end

function Base.pointer(blob::Blob{T}) where T
    convert(Ptr{T}, getfield(blob, :base) + getfield(blob, :offset))
end

function Base.:+(blob::Blob{T}, offset::Integer) where T
    Blob{T}(getfield(blob, :base), getfield(blob, :offset) + offset, getfield(blob, :limit))
end

function Base.:-(blob1::Blob, blob2::Blob)
    assert_same_allocation(blob1, blob2)
    getfield(blob1, :offset) - getfield(blob2, :offset)
end

@inline function boundscheck(blob::Blob{T}) where T
    @boundscheck begin
        if (getfield(blob, :offset) < 0) || (getfield(blob, :offset) + self_size(T) > getfield(blob, :limit))
            throw(BoundsError(blob))
        end
    end
end

Base.@propagate_inbounds function Base.getindex(blob::Blob{T}) where T
    boundscheck(blob)
    unsafe_load(blob)
end

# TODO(jamii) do we need to align data?
"""
    self_size(::Type{T}, args...) where {T}

The number of bytes needed to allocate `T` itself.

Defaults to `sizeof(T)`.
"""
@generated function self_size(::Type{T}) where T
    @assert isconcretetype(T)
    if isempty(fieldnames(T))
        quote
            $(Expr(:meta, :inline))
            $(sizeof(T))
        end
    else
        quote
            $(Expr(:meta, :inline))
            $(+(0, @splice i in 1:length(fieldnames(T)) begin
                self_size(fieldtype(T, i))
            end))
        end
    end
end

function blob_offset(::Type{T}, i::Int) where {T}
    +(0, @splice j in 1:(i-1) begin
        self_size(fieldtype(T, j))
    end)
end

@generated function Base.getindex(blob::Blob{T}, ::Type{Val{field}}) where {T, field}
    i = findfirst(isequal(field), fieldnames(T))
    @assert i != nothing "$T has no field $field"
    quote
        $(Expr(:meta, :inline))
        Blob{$(fieldtype(T, i))}(blob + $(blob_offset(T, i)))
    end
end

Base.@propagate_inbounds function Base.setindex!(blob::Blob{T}, value::T) where T
    boundscheck(blob)
    unsafe_store!(blob, value)
end

# if the value is the wrong type, try to convert it (just like setting a field normally)
Base.@propagate_inbounds function Base.setindex!(blob::Blob{T}, value) where T
    setindex!(blob, convert(T, value))
end

@generated function Base.unsafe_load(blob::Blob{T}) where {T}
    if isempty(fieldnames(T))
        quote
            $(Expr(:meta, :inline))
            unsafe_load(pointer(blob))
        end
    else
        quote
            $(Expr(:meta, :inline))
            $(Expr(:new, T, @splice (i, field) in enumerate(fieldnames(T)) quote
                unsafe_load(getindex(blob, $(Val{field})))
            end))
        end
    end
end

@generated function Base.unsafe_store!(blob::Blob{T}, value::T) where {T}
    if isempty(fieldnames(T))
        quote
            $(Expr(:meta, :inline))
            unsafe_store!(pointer(blob), value)
            value
        end
    elseif T <: Tuple
        quote
            $(Expr(:meta, :inline))
            $(@splice (i, field) in enumerate(fieldnames(T)) quote
                unsafe_store!(getindex(blob, $(Val{field})), value[$field])
            end)
            value
        end
    else
        quote
            $(Expr(:meta, :inline))
            $(@splice (i, field) in enumerate(fieldnames(T)) quote
                unsafe_store!(getindex(blob, $(Val{field})), value.$field)
            end)
            value
        end
    end
end

# if the value is the wrong type, try to convert it (just like setting a field normally)
function Base.unsafe_store!(blob::Blob{T}, value) where {T}
    unsafe_store!(blob, convert(T, value))
end

# syntax sugar

function Base.propertynames(::Blob{T}, private=false) where T
    fieldnames(T)
end

function Base.getproperty(blob::Blob{T}, field::Symbol) where T
    getindex(blob, Val{field})
end

function Base.setproperty!(blob::Blob{T}, field::Symbol, value) where T
    setindex!(blob, Val{field}, value)
end

function rewrite_address(expr)
    if !(expr isa Expr)
        esc(expr)
    elseif expr.head == :.
        (object, field) = expr.args
        if field isa QuoteNode
            fieldname = field.value
        elseif field isa Expr && field.head == :quote
            fieldname = field.args[1]
        else
            error("Impossible?")
        end
        :(getindex($(rewrite_address(object)), $(Val{fieldname})))
    elseif expr.head == :ref
        object = expr.args[1]
        :(getindex($(rewrite_address(object)), $(map(esc, expr.args[2:end])...)))
    elseif expr.head == :macrocall
        rewrite_address(macroexpand(expr))
    else
        error("Don't know how to compute address for $expr")
    end
end

"""
    @a blob.x[2].y

Get a `Blob` pointing at the *address* of `blob.x[2].y`.
"""
macro a(expr)
    rewrite_address(expr)
end

function rewrite_value(expr)
    if (expr isa Expr) && (expr.head == :(=))
        if length(expr.args) == 2
            :(unsafe_store!($(rewrite_address(expr.args[1])), $(esc(expr.args[2]))))
        else
            error("Don't know how to compute assignment $expr")
        end
    else
        :(unsafe_load($(rewrite_address(expr))))
    end
end

"""
    @v blob.x[2].y

Get the *value* at `blob.x[2].y`.

    @v blob.x[2].y = 42

Set the *value* at `blob.x[2].y`.

NOTE macros bind tightly, so:

    # invalid syntax
    @v blob.x[2].y < 42

    # valid syntax
    (@v blob.x[2].y) < 42
"""
macro v(expr)
    rewrite_value(expr)
end

# patch pointers on the fly during load/store!

function self_size(::Type{Blob{T}}) where T
    sizeof(Int64)
end

@inline function Base.unsafe_load(blob::Blob{Blob{T}}) where {T}
    offset = unsafe_load(Blob{Int64}(blob))
    Blob{T}(getfield(blob, :base), getfield(blob, :offset) + offset, getfield(blob, :limit))
end

@inline function Base.unsafe_store!(blob::Blob{Blob{T}}, value::Blob{T}) where {T}
    assert_same_allocation(blob, value)
    offset = getfield(value, :offset) - getfield(blob, :offset)
    unsafe_store!(Blob{Int64}(blob), offset)
end

ext: jl


text:
"""
`module  JAC.FormFactor`
    ... a submodel of JAC that contains all methods for computing atomic form factors for some level(s).
"""
module FormFactor

    using Printf, ..Basics, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..RadialIntegrals, ..TableStrings

    """
    `struct  FormFactor.Settings  <:  AbstractPropertySettings`  ... defines a type for the details and parameters of computing alpha-variation parameters.

        + qList                    ::Array{Float64,1} ... List of q-values in [a.u.]
        + printBefore              ::Bool             ... True if a list of selected levels is printed before the
                                                          actual computations start.
        + levelSelection           ::LevelSelection   ... Specifies the selected levels, if any.
    """
    struct Settings  <:  AbstractPropertySettings
        qList                      ::Array{Float64,1}
        printBefore                ::Bool
        levelSelection             ::LevelSelection
    end

    """
    `FormFactor.Settings()`  ... constructor for an `empty` instance of FormFactor.Settings for the computation of atomic form factors.
    """
    function Settings()
        Settings(Float64[0., 0.1, 1.0, 10.], false, LevelSelection() )
    end

    # `Base.show(io::IO, settings::FormFactor.Settings)`  ... prepares a proper printout of the variable settings::FormFactor.Settings.
    function Base.show(io::IO, settings::FormFactor.Settings)
        println(io, "qList:                    $(settings.qList)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `struct  FormFactor.Outcome`  ... defines a type to keep the outcome of a form-factor computation, such as the standard form factor
                                      as well other results.

        + level                     ::Level              ... Atomic level to which the outcome refers to.
        + qValues                   ::Array{Float64,1}   ... momentum transfer |q|.
        + standardFs                ::Array{Float64,1}   ... standard atomic form factor F(q) for an (assumed) spherical charge distr.
        + modifiedFs                ::Array{Float64,1}   ... modified atomic form factor F(q) for such a charge distribution.
    """
    struct Outcome
        level                       ::Level
        qValues                     ::Array{Float64,1}
        standardFs                  ::Array{Float64,1}
        modifiedFs                  ::Array{Float64,1}
    end

    """
    `FormFactor.Outcome()`
        ... constructor for an `empty` instance of FormFactor.Outcome for the computation of atomic form factors.
    """
    function Outcome()
        Outcome(Level(), Float64[], Float64[], Float64[])
    end

    # `Base.show(io::IO, outcome::FormFactor.Outcome)`  ... prepares a proper printout of the variable outcome::FormFactor.Outcome.
    function Base.show(io::IO, outcome::FormFactor.Outcome)
        println(io, "level:                   $(outcome.level)  ")
        println(io, "qValues:                 $(outcome.qValues)  ")
        println(io, "standardFs:              $(outcome.standardFs)  ")
        println(io, "modifiedFs:              $(outcome.modifiedFs)  ")
    end

    """
    `FormFactor.amplitude(q::Float64, finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)`
        ... to compute the momentum transfer amplitude  <(alpha_f J_f, kappa) J_i || T^(1) || alpha_i J_i>  for the given
            final and initial level. A value::ComplexF64 is returned.
    """
    function amplitude(q::Float64, finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)
        nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
        printstyled("Compute (single-electron) momentum transfer matrix of dimension $nf x $ni in the final- and initial-state bases " *
                    "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
        matrix = zeros(ComplexF64, nf, ni)
        #
        printstyled("done. \n", color=:light_green)
        printstyled("\nWarning -- FormFactor.amplitude():: Not yet implemented.", color=:cyan)
        #
        amplitude = 4.0 + 5.0*im
        #
        if  display
            println("   (Single-electron) Momentum transfer amplitude:  "                                                 *
                    "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || T^(1) ($q a.u.) ||"   *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = $amplitude  ")
        end

        return( amplitude )
    end

    """
    `FormFactor.standardF(q::Float64, level::Level, grid::Radial.Grid)`
         ... to compute the (real) Fourier transform of the (spherically-symmetric) charge distribution. A value::Float64 is returned.
    """
    function standardF(q::Float64, level::Level, grid::Radial.Grid)
        wa = zeros( grid.NoPoints )
        # Compute the total density
        for  sh in level.basis.subshells
            occ  = Basics.computeMeanSubshellOccupation(sh, [level])
            orb  = level.basis.orbitals[sh]
            nrho = length(orb.P)
            for    i = 1:nrho   wa[i] = wa[i] + occ * (orb.P[i]^2 + orb.Q[i]^2)    end
        end
        # Compute the full integrant; the factor 4pi * r^2 is already in the density
        if     q == 0.
        else             for    i = 2:grid.NoPoints   wa[i] = wa[i] / grid.r[i] * sin(q*grid.r[i]) / q    end
        end
        sF = RadialIntegrals.V0(wa, grid.NoPoints, grid::Radial.Grid)

        return( sF )
    end

    """
    `FormFactor.modifiedF(q::Float64, level::Level, grid::Radial.Grid)`
        ... to compute the (real) Fourier transform of the (spherically-symmetric) charge distribution but including a correction
            due to the local (DFS) potential, relative to the rest mass of the electron. A value::Float64 is returned.
    """
    function modifiedF(q::Float64, level::Level, grid::Radial.Grid)
        wa  = zeros( grid.NoPoints )
        wc  = Defaults.getDefaults("speed of light: c")
        pot = Basics.compute("radial potential: Dirac-Fock-Slater", grid, level)
        # Compute the total density
        for  sh in level.basis.subshells
            occ  = Basics.computeMeanSubshellOccupation(sh, [level])
            orb  = level.basis.orbitals[sh]
            ei   = abs(orb.energy)
            nrho = length(orb.P)
            #
            for    i = 2:nrho
                wx = wc^2 / (wc^2 + ei + pot.Zr[i]/grid.r[i])
                if     q == 0.   wa[i] = wa[i] + occ * (orb.P[i]^2 + orb.Q[i]^2) * wx
                else             wa[i] = wa[i] + occ * (orb.P[i]^2 + orb.Q[i]^2) / grid.r[i] * sin(q*grid.r[i]) / q * wx
                end
            end
        end
        mF = RadialIntegrals.V0(wa, grid.NoPoints, grid::Radial.Grid)

        return( mF )
    end

    """
    `FormFactor.computeAmplitudesProperties(outcome::FormFactor.Outcome, nm::Nuclear.Model, grid::Radial.Grid, settings::FormFactor.Settings)`
        ... to compute all form factors, etc. for the given outcome; a newOutcome::FormFactor.Outcome is returned.
    """
    function computeAmplitudesProperties(outcome::FormFactor.Outcome, nm::Nuclear.Model, grid::Radial.Grid, settings::FormFactor.Settings)
        standardFs = Float64[];    modifiedFs = Float64[]
        for  i = 1:length(outcome.qValues)
            sF = FormFactor.standardF(outcome.qValues[i]::Float64, outcome.level::Level, grid::Radial.Grid)
            mF = FormFactor.modifiedF(outcome.qValues[i]::Float64, outcome.level::Level, grid::Radial.Grid)
            push!( standardFs, sF);    push!( modifiedFs, mF)
        end
        newOutcome = FormFactor.Outcome(outcome.level, outcome.qValues, standardFs, modifiedFs)
        return( newOutcome )
    end

    """
    `FormFactor.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::FormFactor.Settings; output=true)`
        ... to compute (as selected) the alpha-variation parameters for the levels of the given multiplet and as specified
            by the given settings. The results are printed in neat tables to screen but nothing is returned otherwise.
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::FormFactor.Settings; output=true)
        println("")
        printstyled("FormFactor.computeOutcomes(): The computation of the atomic form factors starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------- \n", color=:light_green)
        #
        outcomes = FormFactor.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    FormFactor.displayOutcomes(outcomes)    end
        # Calculate all amplitudes and requested properties
        newOutcomes = FormFactor.Outcome[]
        for  outcome in outcomes
            newOutcome = FormFactor.computeAmplitudesProperties(outcome, nm, grid, settings)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        FormFactor.displayResults(stdout, newOutcomes)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    FormFactor.displayResults(iostream, newOutcomes)   end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `FormFactor.determineOutcomes(multiplet::Multiplet, settings::FormFactor.Settings)`
        ... to determine a list of Outcomes's for the computation of the alpha-variation parameters for the given
            multiplet. It takes into account the particular selections and settings. An Array{FormFactor.Outcome,1}
            is returned. Apart from the level specification, all physical properties are set to zero during the
            initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::FormFactor.Settings)
        outcomes = FormFactor.Outcome[]
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, FormFactor.Outcome(level, settings.qList, Float64[], Float64[] ) )
            end
        end
        return( outcomes )
    end

    """
    `FormFactor.displayOutcomes(outcomes::Array{FormFactor.Outcome,1})`
        ... to display a list of levels that have been selected for the computations. A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise.
    """
    function  displayOutcomes(outcomes::Array{FormFactor.Outcome,1})
        nx = 86
        println(" ")
        println("  Selected FormFactor levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * "  q-values  ...";                                                      sb = sb * "   [a.u.]"
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            if  length(outcome.qValues) >= 1   sa = sa * @sprintf("%.2e", outcome.qValues[1] ) * ",  "   end
            if  length(outcome.qValues) >= 2   sa = sa * @sprintf("%.2e", outcome.qValues[2] ) * ",  "   end
            if  length(outcome.qValues) >= 3   sa = sa * @sprintf("%.2e", outcome.qValues[3] ) * ",  "   end
            if  length(outcome.qValues) >  3   sa = sa * "...    "   else   sa = sa * "    "             end
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `FormFactor.displayResults(stream::IO, outcomes::Array{FormFactor.Outcome,1})`
        ... to display the energies, M_ms and F-parameters, etc. for the selected levels. A neat table is printed but
            nothing is returned otherwise.
    """
    function  displayResults(stream::IO, outcomes::Array{FormFactor.Outcome,1})
        nx = 88
        println(stream, " ")
        println(stream, "  Standard and modified atomic form factors [note F^(standard) (0) = N_e]:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=3)
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(14, "q [a.u.]";       na=1)
        sa = sa * TableStrings.center(14, "standard F";     na=1)
        sa = sa * TableStrings.center(14, "modified F";     na=4)
        sb = sb * TableStrings.center(14, "    " ; na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            sa = sa * @sprintf("%.5e", outcome.qValues[1])                                       * "    "
            sa = sa * @sprintf("%.5e", outcome.standardFs[1])                                    * "    "
            sa = sa * @sprintf("%.5e", outcome.modifiedFs[1])                                    * "    "
            println(stream, sa )
            for  i = 2:length(outcome.qValues)
                sa = repeat(" ", 47)
                sa = sa * @sprintf("%.5e", outcome.qValues[i])                                   * "    "
                sa = sa * @sprintf("%.5e", outcome.standardFs[i])                                * "    "
                sa = sa * @sprintf("%.5e", outcome.modifiedFs[i])                                * "    "
                println(stream, sa )
            end
        end
        println(stream, "  ", TableStrings.hLine(nx), "\n\n")
        #
        return( nothing )
    end

end # module


text:
tab(n::Int) = " "^n
ignore_type(type) = type in [Any, String] || type <: AbstractString

function has_args(cmd::LeafCommand)
    !isempty(cmd.args) || !isnothing(cmd.vararg)
end

section(io) = print(io, "\n\n")
function section(io::IO, title)
    section(io)
    printstyled(io, title; bold=true)
    section(io)
end

Base.@kwdef mutable struct Indent
    text::Int = 0
    desc::Int = 0
end

Base.@kwdef mutable struct Color
    name::Symbol = :light_blue
    args::Symbol = :light_magenta
    dash::Symbol = :light_cyan
end

Base.@kwdef mutable struct Terminal
    width::Int = displaysize(stdout)[2]
    left::Int = floor(Int, 0.4 * width) # left column max width
    right::Int = floor(Int, 0.4 * width) # right column max width

    color::Color = Color()
    indent::Indent = Indent()

    brief::Bool = true
end

Base.show(io::IO, ::MIME"text/plain", cmd::ComoniconExpr) = print_cmd(io, cmd)

function Base.show(io::IO, ::MIME"text/plain", cmd::Description)
    printstyled(io, "brief:\n"; color=:light_black)
    println(io, cmd.brief)
    printstyled(io, "content:\n"; color=:light_black)
    print(io, cmd.content)
end

print_cmd(cmd) = print_cmd(stdout, cmd)
print_cmd(io::IO, cmd) = print_cmd(io, cmd, Terminal())
print_cmd(cmd, t::Terminal) = print_cmd(stdout, cmd, t)

function print_cmd(io::IO, arg::Argument, t::Terminal)
    color = t.color.args
    arg.require || printstyled(io, "["; color)
    arg.require && printstyled(io, "<"; color)
    printstyled(io, arg.name; color)
    if !ignore_type(arg.type)
        printstyled(io, "::", arg.type; color)
    end
    arg.vararg && printstyled(io, "..."; color)
    arg.require && printstyled(io, ">"; color)
    arg.require || printstyled(io, "]"; color)
    return
end

print_cmd(io::IO, cmd::Flag, t::Terminal) = _print_dash(io, cmd, t)

function print_cmd(io::IO, cmd::Option, t::Terminal)
    _print_dash(io, cmd, t)
    isnothing(cmd.hint) ||
        printstyled(io, tab(1), "<", cmd.hint, ">"; color=t.color.args)
    return
end

function _print_dash(io::IO, cmd::Union{Option, Flag}, t::Terminal)
    color = t.color.dash
    if cmd.short
        printstyled(io, "-", first(cmd.name), ", "; color)
    end
    printstyled(io, "--", cmd.name; color)
end

function print_content(io::IO, desc::Description, t::Terminal)
    isnothing(desc.content) || print_within(io, desc.content, t.width, 0)
    return
end

function print_cmd(io::IO, cmd::Entry, t::Terminal)
    section(io)
    printstyled(io, tab(2), cmd.root.name; color=t.color.name, bold=true)
    isnothing(cmd.version) || print(io, " v", cmd.version)
    section(io)
    t.brief = false
    # print description ahead
    print_content(io, cmd.root.description, t)
    section(io, "Usage")
    print_head(io, cmd, t)
    print_body(io, cmd, t)
end

function print_cmd(io::IO, cmd::NodeCommand, t::Terminal)
    section(io)
    print_head(io, cmd, t)
    section(io)
    print_content(io, cmd.description, t)
    print_body(io, cmd, t)
end

function print_cmd(io::IO, cmd::LeafCommand, t::Terminal)
    section(io)
    print_head(io, cmd, t)
    section(io)
    print_content(io, cmd.description, t)
    print_body(io, cmd, t)
end

print_head(io::IO, cmd::Entry, t::Terminal) = print_head(io, cmd.root, t)

function print_head(io::IO, cmd::NodeCommand, t::Terminal)
    print(io, tab(2))
    print_signature(io, cmd, t)
end

function print_head(io::IO, cmd::LeafCommand, t::Terminal)
    print(io, tab(2))
    print_name(io, cmd, t)
    if has_args(cmd)
        printstyled(io, tab(1), "<args>"; color=t.color.args)
    end

    isempty(cmd.args) || printstyled(io, tab(1), "[options]"; color = t.color.dash)
    isempty(cmd.flags) || printstyled(io, tab(1), "[flags]"; color = t.color.dash)
    return
end

function print_name(io::IO, cmd, t::Terminal)
    printstyled(io, cmd.name; color=t.color.name, bold=true)
end

function print_signature(io::IO, cmd, t::Terminal)
    print_cmd(io, cmd, t)
end

function print_signature(io::IO, cmd::NodeCommand, t::Terminal)
    print_name(io, cmd, t)
    printstyled(io, tab(1), "<command>"; color=t.color.name)
end

function print_signature(io::IO, cmd::LeafCommand, t::Terminal)
    print_name(io, cmd, t)
    for each in cmd.args
        print(io, tab(1))
        print_cmd(io, each, t)
    end

    if !isnothing(cmd.vararg)
        print(io, tab(1))
        print_cmd(io, cmd.vararg, t)
    end
end

function print_body(io::IO, cmd::Entry, t::Terminal)
    print_body(io, cmd.root, t)
    version_flag = "-V, --version"
    printstyled(io, tab(2), version_flag; color=t.color.dash)
    print_indent_content(io, "print version information", t, length(version_flag)+2)
    println(io)
end

function print_body(io::IO, cmd::NodeCommand, t::Terminal)
    section(io, "Commands")
    for each in values(cmd.subcmds)
        print_sig_brief(io, each, t)
        println(io)
    end
    section(io, "Flags")
    print_help(io, t)
    return
end

function print_body(io::IO, cmd::LeafCommand, t::Terminal)
    if has_args(cmd)
        section(io, "Args")
    end

    for each in cmd.args
        print_sig_brief(io, each, t)
        println(io)
    end

    if !isnothing(cmd.vararg)
        print_sig_brief(io, cmd.vararg, t)
        println(io)
    end

    if !isempty(cmd.options)
        section(io, "Options")
        for each in unique(values(cmd.options))
            print_sig_brief(io, each, t)
            println(io)
        end
    end

    section(io, "Flags")

    for each in unique(values(cmd.flags))
        print_sig_brief(io, each, t)
        section(io)
    end
    print_help(io, t)
end

function print_help(io::IO, t::Terminal)
    help_flag = "-h, --help"
    printstyled(io, tab(2), help_flag; color=t.color.dash)
    print_indent_content(io, "print this help message", t, length(help_flag)+2)
    println(io)
end

function print_sig_brief(io::IO, cmd, t::Terminal)
    buf = IOBuffer()
    print_signature(buf, cmd, t)
    s = String(take!(buf))

    print(io, tab(2))
    print_signature(io, cmd, t)
    isnothing(cmd.description.brief) && return

    print_indent_content(io, cmd.description.brief, t, length(s)+2)
    return
end

function print_indent_content(io::IO, text::String, t::Terminal, firstline::Int)
    middle = t.width - t.left - t.right
    lines = splitlines(text, t.width)
    isempty(lines) && return
    print(io, tab(t.left - firstline + middle), lines[1])
    for i in 2:length(lines)
        print(io, tab(t.width - t.right), lines[i])
        if i !== lastindex(lines)
            println(io)
        end
    end
    return
end

function print_within(io::IO, text::String, width::Int, indent::Int)
    lines = splitlines(text, width - indent)
    for i in eachindex(lines)
        print(io, tab(indent), lines[i])
        if i !== lastindex(lines)
            println(io)
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, DataFrames
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 2.2"

# ╔═╡ 6ffea8a0-1ea1-11eb-1991-61c054bcac44
begin
	A = Set(['a','e','i','o','u'])
	B = Set(['x','y','z'])
	omega = 'a':'z'

	N = 10^6
end;

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	df = DataFrame()
	local mcEst1, mcEst2
	for _ in 1:5
    	mcEst1 = sum([in(sample(omega),A) || in(sample(omega),B) for _ in 1:N])/N
    	mcEst2 = sum([in(sample(omega),union(A,B)) for _ in 1:N])/N
		append!(df, DataFrame(:mcEst1 => mcEst1, :mcEst2 => mcEst2))
	end
	df
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 2.2"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═6ffea8a0-1ea1-11eb-1991-61c054bcac44
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
import GridWorlds
import GridWorlds: GW
import ReinforcementLearningBase
import ReinforcementLearningBase: RLBase
import BenchmarkTools
const BT = BenchmarkTools
import Random
import Dates

const STEPS_PER_RESET = 100
const NUM_RESETS = 100
const SEED = 0

const rng = Random.MersenneTwister(SEED)
const information = Dict()

ENVS = [GW.EmptyRoomDirected,
        GW.EmptyRoomUndirected,
        GW.GridRoomsDirected,
        GW.GridRoomsUndirected,
        GW.SequentialRoomsDirected,
        GW.SequentialRoomsUndirected,
        GW.MazeDirected,
        GW.MazeUndirected,
        GW.GoToTargetDirected,
        GW.GoToTargetUndirected,
        GW.DoorKeyDirected,
        GW.DoorKeyUndirected,
        GW.CollectGemsDirected,
        GW.CollectGemsUndirected,
        GW.DynamicObstaclesDirected,
        GW.DynamicObstaclesUndirected,
        GW.SokobanDirected,
        GW.SokobanUndirected,
        GW.Snake,
        GW.Catcher,
        GW.TransportDirected,
        GW.TransportUndirected,
       ]

function run_experiment(rng, Env, num_resets, steps_per_reset)
    env = Env(rng = rng)

    for _ in 1:num_resets
        RLBase.reset!(env)
        total_reward = 0
        for _ in 1:steps_per_reset
            state = RLBase.state(env)
            action = rand(RLBase.action_space(env))
            env(action)
            is_terminated = RLBase.is_terminated(env)
            reward = RLBase.reward(env)
            total_reward += reward
        end
    end

    return nothing
end

function format_benchmark(str::String)
    l = split(str, "\n")
    deleteat!(l, (1, 4, 9))
    return strip.(l)
end

function write_benchmarks(information; file = "benchmark.md")
    io = open(file, "w")

    write(io, "Date: " * Dates.format(Dates.now(), "yyyy_mm_dd_HH_MM_SS") * "\n")
    write(io, "# List of Environments\n")

    for Env in ENVS
        name = Env.body.body.name.name
        write(io, "  1. [$(String(name))](#$(lowercase(String(name))))\n")
    end

    write(io, "\n")
    write(io, "# Benchmarks\n\n")

    for Env in ENVS
        name = Env.body.body.name.name
        env_benchmark = information[name]

        write(io, "# $(String(name))\n\n")

        write(io, "#### Run uniformly random policy, NUM_RESETS = $(NUM_RESETS), STEPS_PER_RESET = $(STEPS_PER_RESET), TOTAL_STEPS = $(NUM_RESETS * STEPS_PER_RESET)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:run_experiment]))
            write(io, line * "\n\n")
        end

        write(io, "#### $(String(Symbol(Env)))()\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:instantiation]))
            write(io, line * "\n\n")
        end

        write(io, "#### RLBase.reset!(env)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:reset!]))
            write(io, line * "\n\n")
        end

        write(io, "#### RLBase.state(env)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:state]))
            write(io, line * "\n\n")
        end

        write(io, "#### RLBase.action_space(env)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:action_space]))
            write(io, line * "\n\n")
        end

        write(io, "#### RLBase.is_terminated(env)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:is_terminated]))
            write(io, line * "\n\n")
        end

        write(io, "#### RLBase.reward(env)\n\n")
        for line in format_benchmark(repr("text/plain", env_benchmark[:reward]))
            write(io, line * "\n\n")
        end

        for action in keys(env_benchmark[:action_info])
            write(io, "#### env($action)\n\n")
            for line in format_benchmark(repr("text/plain", env_benchmark[:action_info][action]))
                write(io, line * "\n\n")
            end
        end

    end

    close(io)
end

# compile everything once
for Env in ENVS
    run_experiment(rng, Env, NUM_RESETS, STEPS_PER_RESET)
end

@info "First run (for compilation) is complete"

for Env in ENVS

    env_benchmark = Dict()

    env_benchmark[:run_experiment] = BT.@benchmark run_experiment($(Ref(rng))[], $(Ref(Env))[], $(Ref(NUM_RESETS))[], $(Ref(STEPS_PER_RESET))[])

    env_benchmark[:instantiation] = BT.@benchmark $(Ref(Env))[](rng = $(Ref(rng))[])

    env = Env(rng = rng)

    env_benchmark[:reset!] = BT.@benchmark RLBase.reset!($(Ref(env))[])
    env_benchmark[:state] = BT.@benchmark RLBase.state($(Ref(env))[])
    env_benchmark[:action_space] = BT.@benchmark RLBase.action_space($(Ref(env))[])
    env_benchmark[:is_terminated] = BT.@benchmark RLBase.is_terminated($(Ref(env))[])
    env_benchmark[:reward] = BT.@benchmark RLBase.reward($(Ref(env))[])

    action_info = Dict()
    for action in RLBase.action_space(env)
        action_info[Symbol(action)] = BT.@benchmark $(Ref(env))[]($(Ref(action))[])
    end
    env_benchmark[:action_info] = action_info

    name = Env.body.body.name.name
    information[name] = env_benchmark

    @info "$(name) benchmark complete"
end

write_benchmarks(information)


text:
export stylesheet

function stylesheet(href::String; args...) :: String
  Genie.Renderer.Html.link(href=href, rel="stylesheet", args...)
end


text:
# Standard isentropic vortex test case.  For a more complete description of
# the setup see for Example 3 of:
#
# @article{ZHOU2003159,
#   author = {Y.C. Zhou and G.W. Wei},
#   title = {High resolution conjugate filters for the simulation of flows},
#   journal = {Journal of Computational Physics},
#   volume = {189},
#   number = {1},
#   pages = {159--179},
#   year = {2003},
#   doi = {10.1016/S0021-9991(03)00206-7},
#   url = {https://doi.org/10.1016/S0021-9991(03)00206-7},
# }
#
# This version runs the isentropic vortex as a stand alone test (no dependence
# on CLIMA moist thermodynamics)

using MPI
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E)
const statenames = ("ρ", "U", "V", "W", "E")
const γ_exact = 7 // 5
if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
  using Random
end

# preflux computation
@inline function preflux(Q, _...)
  γ::eltype(Q) = γ_exact
  @inbounds ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
  ρinv = 1 / ρ
  u, v, w = ρinv * U, ρinv * V, ρinv * W
  ((γ-1)*(E - ρinv * (U^2 + V^2 + W^2) / 2), u, v, w, ρinv)
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
  γ::eltype(Q) = γ_exact
  @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# physical flux function
eulerflux!(F, Q, QV, aux, t) =
eulerflux!(F, Q, QV, aux, t, preflux(Q)...)

@inline function eulerflux!(F, Q, QV, aux, t, P, u, v, w, ρinv)
  @inbounds begin
    ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]

    F[1, _ρ], F[2, _ρ], F[3, _ρ] = U          , V          , W
    F[1, _U], F[2, _U], F[3, _U] = u * U  + P , v * U      , w * U
    F[1, _V], F[2, _V], F[3, _V] = u * V      , v * V + P  , w * V
    F[1, _W], F[2, _W], F[3, _W] = u * W      , v * W      , w * W + P
    F[1, _E], F[2, _E], F[3, _E] = u * (E + P), v * (E + P), w * (E + P)
  end
end

# initial condition
const halfperiod = 5
function isentropicvortex!(Q, t, x, y, z)
  DFloat = eltype(Q)

  γ::DFloat    = γ_exact
  uinf::DFloat = 2
  vinf::DFloat = 1
  Tinf::DFloat = 1
  λ::DFloat    = 5

  xs = x - uinf*t
  ys = y - vinf*t

  # make the function periodic
  xtn = floor((xs+halfperiod)/(2halfperiod))
  ytn = floor((ys+halfperiod)/(2halfperiod))
  xp = xs - xtn*2*halfperiod
  yp = ys - ytn*2*halfperiod

  rsq = xp^2 + yp^2

  u = uinf - λ*(1//2)*exp(1-rsq)*yp/π
  v = vinf + λ*(1//2)*exp(1-rsq)*xp/π
  w = zero(DFloat)

  ρ = (Tinf - ((γ-1)*λ^2*exp(2*(1-rsq))/(γ*16*π*π)))^(1/(γ-1))
  p = ρ^γ
  U = ρ*u
  V = ρ*v
  W = ρ*w
  E = p/(γ-1) + (1//2)*ρ*(u^2 + v^2 + w^2)

  if integration_testing
    @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E] = ρ, U, V, W, E
  else
    @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E] =
    10+rand(), rand(), rand(), rand(), 10+rand()
  end

end

function main(mpicomm, DFloat, topl::AbstractTopology{dim}, N, timeend,
              ArrayType, dt) where {dim}

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  # spacedisc = data needed for evaluating the right-hand side function
  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = _nstate,
                           flux! = eulerflux!,
                           numerical_flux! = (x...) ->
                           NumericalFluxes.rusanov!(x..., eulerflux!,
                                                    wavespeed,
                                                    preflux))

  # This is a actual state/function that lives on the grid
  initialcondition(Q, x...) = isentropicvortex!(Q, DFloat(0), x...)
  Q = MPIStateArray(spacedisc, initialcondition)

  lsrk = LowStorageRungeKutta(spacedisc, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf """Update
  simtime = %.16e
  runtime = %s
  norm(Q) = %.16e""" ODESolvers.gettime(lsrk) Dates.format(convert(Dates.DateTime, Dates.now()-starttime[]), Dates.dateformat"HH:MM:SS") energy
    end
  end

  #= Paraview calculators:
  P = (0.4) * (E  - (U^2 + V^2 + W^2) / (2*ρ) - 9.81 * ρ * coordsZ)
  theta = (100000/287.0024093890231) * (P / 100000)^(1/1.4) / ρ
  =#
  step = [0]
  mkpath("vtk")
  cbvtk = GenericCallbacks.EveryXSimulationSteps(100) do (init=false)
    outprefix = @sprintf("vtk/isentropicvortex_%dD_mpirank%04d_step%04d",
                         dim, MPI.Comm_rank(mpicomm), step[1])
    @debug "doing VTK output" outprefix
    DGBalanceLawDiscretizations.writevtk(outprefix, Q, spacedisc, statenames)
    step[1] += 1
    nothing
  end

  # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  if integration_testing
    Qe = MPIStateArray(spacedisc,
                       (Q, x...) -> isentropicvortex!(Q, DFloat(timeend), x...))
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf/eng0 engf-eng0 errf errf / engfe
  else
    @info @sprintf """Finished
    norm(Q)            = %.16e
    norm(Q) / norm(Q₀) = %.16e
    norm(Q) - norm(Q₀) = %.16e""" engf engf/eng0 engf-eng0
  end
  integration_testing ? errf : (engf / eng0)
end

function run(mpicomm, dim, Ne, N, timeend, DFloat, dt)
  ArrayType = Array

  brickrange = ntuple(j->range(DFloat(-halfperiod); length=Ne[j]+1,
                               stop=halfperiod), dim)
  topl = BrickTopology(mpicomm, brickrange, periodicity=ntuple(j->true, dim))
  main(mpicomm, DFloat, topl, N, timeend, ArrayType, dt)
end

using Test
let
  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())
  mpicomm = MPI.COMM_WORLD
  if MPI.Comm_rank(mpicomm) == 0
    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
    global_logger(ConsoleLogger(stderr, loglevel))
  else
    global_logger(NullLogger())
  end

  if integration_testing
    timeend = 1
    numelem = (5, 5, 1)

    polynomialorder = 4

    expected_error = Array{Float64}(undef, 2, 3) # dim-1, lvl
    expected_error[1,1] = 5.7115689019456495e-01
    expected_error[1,2] = 6.9418982796523573e-02
    expected_error[1,3] = 3.2927550219067014e-03
    expected_error[2,1] = 1.8061566743070110e+00
    expected_error[2,2] = 2.1952209848920567e-01
    expected_error[2,3] = 1.0412605646145325e-02
    lvls = size(expected_error, 2)

    for DFloat in (Float64,) #Float32)
      for dim = 2:3
        err = zeros(DFloat, lvls)
        for l = 1:lvls
          Ne = ntuple(j->2^(l-1) * numelem[j], dim)
          dt = 1e-2 / Ne[1]
          nsteps = ceil(Int64, timeend / dt)
          dt = timeend / nsteps
          err[l] = run(mpicomm, dim, Ne, polynomialorder, timeend, DFloat, dt)
          @test err[l] ≈ DFloat(expected_error[dim-1, l])
        end
        @info begin
          msg = ""
          for l = 1:lvls-1
            rate = log2(err[l]) - log2(err[l+1])
            msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
          end
          msg
        end
      end
    end
  else
    numelem = (3, 4, 5)
    dt = 1e-3
    timeend = 2dt

    polynomialorder = 4

    mpicomm = MPI.COMM_WORLD

    check_engf_eng0 = Dict{Tuple{Int64, Int64, DataType}, AbstractFloat}()
    check_engf_eng0[2, 1, Float64] = 9.9999795068862996e-01
    check_engf_eng0[3, 1, Float64] = 9.9999641494886327e-01
    check_engf_eng0[2, 3, Float64] = 9.9999876109562658e-01
    check_engf_eng0[3, 3, Float64] = 9.9999654059181553e-01

    for DFloat in (Float64,) #Float32)
      for dim = 2:3
        Random.seed!(0)
        engf_eng0 = run(mpicomm, dim, numelem[1:dim], polynomialorder, timeend,
                        DFloat, dt)
        @test check_engf_eng0[dim, MPI.Comm_size(mpicomm), DFloat] ≈ engf_eng0
      end
    end
  end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
module CategoricalAlgebra

include("ShapeDiagrams.jl")
include("Matrices.jl")
include("FinSets.jl")
include("FinRelations.jl")
include("Permutations.jl")
include("CSets.jl")
include("Graphs.jl")

end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export ChemfilesError

struct ChemfilesError <: Exception
    message::AbstractString
end

Base.show(io::IO, e::ChemfilesError) = show(io, "Chemfiles error: $(e.message)")

"""
Get the last error message from the chemfiles runtime.
"""
function last_error()
    unsafe_string(lib.chfl_last_error())
end

"""
Clear any error messages stored by the chemfiles runtime.
"""
function clear_errors()
    __check(lib.chfl_clear_errors())
end

"""
Default warning callback for Chemfiles
"""
function __default_warning_callback(message::String)
    @warn "$message"
end

# Store the warning callback in a global to prevent garbage collection
WARNING_CALLBACK = nothing

function _warning_callback_adaptator(message)
    try
        WARNING_CALLBACK(unsafe_string(message))
    catch error
        @warn "caught $error in warning callback"
    end
end

"""
Set the global warning ``callback`` to be used for each warning event.

The ``callback`` function must take a ``String`` and return nothing.
"""
function set_warning_callback(callback::Function)
    global WARNING_CALLBACK
    WARNING_CALLBACK = callback
    ptr = @cfunction(_warning_callback_adaptator, Cvoid, (Ptr{UInt8},))
    __check(lib.chfl_set_warning_callback(ptr))
end

"""
Read configuration data from the file at `path`.

By default, chemfiles reads configuration from any file name `.chemfilesrc`
in the current directory or any parent directory. This function can be used
to add data from another configuration file.

This function will fail if there is no file at `path`, or if the file is
incorectly formatted. Data from the new configuration file will overwrite
any existing data.
"""
function add_configuration(path::String)
    __check(lib.chfl_add_configuration(pointer(path)))
end


text:
# save to disk/read from disk
savfile1 = "test.evt"

printstyled("  read/write of EventTraceData with compression\n", color=:light_green)
SeisIO.KW.comp = 0x02
S = convert(EventTraceData, randSeisData())
wseis(savfile1, S)
R = rseis(savfile1)[1]
@test R == S

SeisIO.KW.comp = 0x01
S = convert(EventTraceData, randSeisData())
C = convert(EventChannel, SeisChannel())
nx = SeisIO.KW.n_zip*2
C.t = [1 0; nx 0]
C.x = randn(nx)
push!(S, C)
wseis(savfile1, S)
R = rseis(savfile1)[1]
@test R == S

rm(savfile1)


text:
for F in (:isfinite, :isinf, :isnan, :issubnormal)
    @eval begin
        $F(x::T) where {T<:MarkableFloat64} = $F(Float64(x))
        $F(x::T) where {T<:MarkableFloat32} = $F(Float32(x))
        $F(x::T) where {T<:MarkableFloat16} = $F(Float16(x))
    end
end


text:
using Mimi
using MimiDICE2016R2

DICE = MimiDICE2016R2.get_model()
run(DICE)

explore(DICE)


text:
"""
    deleteat!(ds::Dataset, inds)

Delete rows specified by `inds` from a `Dataset` `ds` in place and return it.

`inds` must be:
a vector of sorted and unique integers, a boolean vector, an integer, or `Not`.

# Examples
```jldoctest
julia> ds = Dataset(a=1:3, b=4:6)
3×2 Dataset
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     1      4
   2 │     2      5
   3 │     3      6

julia> deleteat!(ds, 2)
2×2 Dataset
 Row │ a      b
     │ Int64  Int64
─────┼──────────────
   1 │     1      4
   2 │     3      6
```

"""
function Base.deleteat!(ds::Dataset, inds)

# Modify Dataset
    if !isempty(inds) && size(ds, 2) == 0
        throw(BoundsError(ds, (inds, :)))
    end

    # we require ind to be stored and unique like in Base
    # otherwise an error will be thrown and the data set will get corrupted
    return _delete!_helper(ds, inds)
end

# Modify Dataset
function Base.deleteat!(ds::Dataset, inds::AbstractVector{Bool})
    if length(inds) != size(ds, 1)
        throw(BoundsError(ds, (inds, :)))
    end
    # drop = _findall(inds)
    return _delete!_helper(ds, inds)
end

# Modify Dataset
Base.deleteat!(ds::Dataset, inds::Not) = deleteat!(ds, axes(ds, 1)[inds])

# Modify Dataset
function _delete!_helper(ds::Dataset, drop)
    cols = _columns(ds)
    isempty(cols) && return ds

    n = nrow(ds)
    col1 = cols[1]
    deleteat!(col1, drop)
    newn = length(col1)

    for i in 2:length(cols)
        col = cols[i]
        if length(col) == n
            deleteat!(col, drop)
        end
    end

    for i in 1:length(cols)
        # this should never happen, but we add it for safety
        @assert length(cols[i]) == newn corrupt_msg(ds, i)
    end
    _reset_grouping_info!(ds)
    _modified(_attributes(ds))
    return ds
end

"""
    empty!(ds::Dataset)

Remove all rows from `ds`, making each of its columns empty.
"""
function Base.empty!(ds::Dataset)

# Modify Dataset
    foreach(empty!, _columns(ds))
    _reset_grouping_info!(ds)
    _modified(_attributes(ds))
    return ds
end

"""
    append!(ds::Dataset, ds2::AbstractDataset; cols::Symbol=:setequal,
            promote::Bool=(cols in [:union, :subset]))
    append!(ds::Dataset, table; cols::Symbol=:setequal,
            promote::Bool=(cols in [:union, :subset]))

Add the rows of `ds2` to the end of `ds`. If the second argument `table` is not an
`AbstractDataset` then it is converted using `Dataset(table, copycols=false)`
before being appended.

The exact behavior of `append!` depends on the `cols` argument:
* If `cols == :setequal` (this is the default)
  then `ds2` must contain exactly the same columns as `ds` (but possibly in a
  different order).
* If `cols == :orderequal` then `ds2` must contain the same columns in the same
  order (for `AbstractDict` this option requires that `keys(row)` matches
  `propertynames(ds)` to allow for support of ordered dicts; however, if `ds2`
  is a `Dict` an error is thrown as it is an unordered collection).
* If `cols == :intersect` then `ds2` may contain more columns than `ds`, but all
  column names that are present in `ds` must be present in `ds2` and only these
  are used.
* If `cols == :subset` then `append!` behaves like for `:intersect` but if some
  column is missing in `ds2` then a `missing` value is pushed to `ds`.
* If `cols == :union` then `append!` adds columns missing in `ds` that are present
  in `ds2`, for columns present in `ds` but missing in `ds2` a `missing` value
  is pushed.

If `promote=true` and element type of a column present in `ds` does not allow
the type of a pushed argument then a new column with a promoted element type
allowing it is freshly allocated and stored in `ds`. If `promote=false` an error
is thrown.

The above rule has the following exceptions:
* If `ds` has no columns then copies of columns from `ds2` are added to it.
* If `ds2` has no columns then calling `append!` leaves `ds` unchanged.

Please note that `append!` must not be used on a `Dataset` that contains
columns that are aliases (equal when compared with `===`).

# See also

Use [`push!`](@ref) to add individual rows to a data set and [`vcat`](@ref)
to vertically concatenate data sets.

# Examples
```jldoctest
julia> ds1 = Dataset(A=1:3, B=1:3)
3×2 Dataset
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3

julia> ds2 = Dataset(A=4.0:6.0, B=4:6)
3×2 Dataset
 Row │ A        B
     │ Float64  Int64
─────┼────────────────
   1 │     4.0      4
   2 │     5.0      5
   3 │     6.0      6

julia> append!(ds1, ds2);

julia> ds1
6×2 Dataset
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     4      4
   5 │     5      5
   6 │     6      6
```
"""
function Base.append!(ds1::Dataset, ds2::AbstractDataset; cols::Symbol=:setequal,
                      promote::Bool=(cols in [:union, :subset]))
# appending ds2 to ds1 should keep the formats of ds1, and metadata of ds1, if ds1 is
# empty then formats of ds2 transfer to ds1, but not metadata
# Modify Dataset
    if !(cols in (:orderequal, :setequal, :intersect, :subset, :union))
        throw(ArgumentError("`cols` keyword argument must be " *
                            ":orderequal, :setequal, :intersect, :subset or :union)"))
    end

    if ncol(ds1) == 0
        for (n, v) in pairs(eachcol(ds2))
            f_col = getformat(ds2, n)
            ds1[!, n] = copy(v) # make sure ds1 does not reuse ds2
            setformat!(ds1, n => f_col)
        end
        _copy_grouping_info!(ds1, ds2)
        setinfo!(ds1, _attributes(ds2).meta.info[])
        return ds1
    end
    ncol(ds2) == 0 && return ds1

    if cols == :orderequal && _names(ds1) != _names(ds2)
        wrongnames = symdiff(_names(ds1), _names(ds2))
        if isempty(wrongnames)
            mismatches = findall(_names(ds1) .!= _names(ds2))
            @assert !isempty(mismatches)
            throw(ArgumentError("Columns number " *
                                join(mismatches, ", ", " and ") *
                                " do not have the same names in both passed " *
                                "data sets and `cols == :orderequal`"))
        else
            mismatchmsg = " Column names :" *
            throw(ArgumentError("Column names :" *
                                join(wrongnames, ", :", " and :") *
                                " were found in only one of the passed data sets " *
                                "and `cols == :orderequal`"))
        end
    elseif cols == :setequal
        wrongnames = symdiff(_names(ds1), _names(ds2))
        if !isempty(wrongnames)
            throw(ArgumentError("Column names :" *
                                join(wrongnames, ", :", " and :") *
                                " were found in only one of the passed data sets " *
                                "and `cols == :setequal`"))
        end
    elseif cols == :intersect
        wrongnames = setdiff(_names(ds1), _names(ds2))
        if !isempty(wrongnames)
            throw(ArgumentError("Column names :" *
                                join(wrongnames, ", :", " and :") *
                                " were found in only in destination data set " *
                                "and `cols == :intersect`"))
        end
    end

    nrows, ncols = size(ds1)
    targetrows = nrows + nrow(ds2)
    current_col = 0
    # in the code below we use a direct access to _columns because
    # we resize the columns so temporarily the `Dataset` is internally
    # inconsistent and normal data set indexing would error.

    # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    current_modified_time = _attributes(ds1).meta.modified[]
    try
        for (j, n) in enumerate(_names(ds1))
            current_col += 1
            format_of_cur_col = getformat(ds1, n)
            if hasproperty(ds2, n)
                ds2_c = _columns(ds2)[index(ds2)[n]]
                S = eltype(ds2_c)
                ds1_c = _columns(ds1)[j]
                T = eltype(ds1_c)
                if S <: T || !promote || promote_type(S, T) <: T
                    # if S <: T || promote_type(S, T) <: T this should never throw an exception
                    append!(ds1_c, ds2_c)
                else
                    newcol = similar(ds1_c, promote_type(S, T), targetrows)
                    copyto!(newcol, 1, ds1_c, 1, nrows)
                    copyto!(newcol, nrows+1, ds2_c, 1, targetrows - nrows)
                    firstindex(newcol) != 1 && _onebased_check_error()
                    _columns(ds1)[j] = newcol
                end
            else
                if Missing <: eltype(ds1[!, j])
                    resize!(_columns(ds1)[j], targetrows)
                    ds1[nrows+1:targetrows, j] .= missing
                elseif promote
                    newcol = similar(_columns(ds1)[j], Union{Missing, eltype(ds1[!, j])},
                                     targetrows)
                    copyto!(newcol, 1, _columns(ds1)[j], 1, nrows)
                    newcol[nrows+1:targetrows] .= missing
                    firstindex(newcol) != 1 && _onebased_check_error()
                    _columns(ds1)[j] = newcol
                else
                    throw(ArgumentError("promote=false and source data set does " *
                                        "not contain column :$n, while destination " *
                                        "column does not allow for missing values"))
                end
            end
            setformat!(ds1, n => format_of_cur_col)
            _reset_grouping_info!(ds1)
            _modified(_attributes(ds1))
        end
        current_col = 0
        for col in _columns(ds1)
            current_col += 1
            @assert length(col) == targetrows
        end
        if cols == :union
            for n in setdiff(_names(ds2), _names(ds1))
                newcol = similar(_columns(ds2)[index(ds2)[n]], Union{Missing, eltype(ds2[!, n])},
                                 targetrows)
                @inbounds newcol[1:nrows] .= missing
                copyto!(newcol, nrows+1, _columns(ds2)[index(ds2)[n]], 1, targetrows - nrows)
                ds1[!, n] = newcol
                _modified(_attributes(ds1))
                _reset_grouping_info!(ds1)
            end
        end
    catch err
        # Undo changes in case of error
        for col in _columns(ds1)
            resize!(col, nrows)
        end
        # go back to original modified time
        _attributes(ds1).meta.modified[] = current_modified_time
        @error "Error adding value to column :$(_names(ds1)[current_col])."
        rethrow(err)
    end
    return ds1
end

# TODO needs more works on how to handle formats??
# push! always reset the grouping information (??)
# Modify Dataset
function Base.push!(ds::Dataset, row::Union{AbstractDict, NamedTuple};
                    cols::Symbol=:setequal,
                    promote::Bool=(cols in [:union, :subset]))
    # push keep formats
    possible_cols = (:orderequal, :setequal, :intersect, :subset, :union)
    if !(cols in possible_cols)
        throw(ArgumentError("`cols` keyword argument must be any of :" *
                            join(possible_cols, ", :")))
    end

    nrows, ncols = size(ds)
    targetrows = nrows + 1
    # here the formats should be kept, setproperty! modifies time
    if ncols == 0 && row isa NamedTuple
        for (n, v) in pairs(row)
            format_of_cur_col = getformat(ds, n)
            setproperty!(ds, n, fill!(Tables.allocatecolumn(typeof(v), 1), v))
            setformat!(ds, n => format_of_cur_col)
        end
        _reset_grouping_info!(ds)
        return ds
    end

    old_row_type = typeof(row)
    if row isa AbstractDict && keytype(row) !== Symbol &&
        (keytype(row) <: AbstractString || all(x -> x isa AbstractString, keys(row)))
        row = (;(Symbol.(keys(row)) .=> values(row))...)
    end

    # in the code below we use a direct access to _columns because
    # we resize the columns so temporarily the `Dataset` is internally
    # inconsistent and normal data set indexing would error.
    if cols == :union
        current_modified = _attributes(ds).meta.modified[]
        if row isa AbstractDict && keytype(row) !== Symbol && !all(x -> x isa Symbol, keys(row))
            throw(ArgumentError("when `cols == :union` all keys of row must be Symbol"))
        end
        for (i, colname) in enumerate(_names(ds))
            format_of_cur_col = getformat(ds, colname)
            col = _columns(ds)[i]
            if haskey(row, colname)
                val = row[colname]
            else
                val = missing
            end
            S = typeof(val)
            T = eltype(col)
            if S <: T || promote_type(S, T) <: T
                push!(col, val)
            elseif !promote
                try
                    push!(col, val)
                catch err
                    setformat!(ds, colname => format_of_cur_col)
                    for col in _columns(ds)
                        resize!(col, nrows)
                    end
                    _attributes(ds).meta.modified[] = current_modified
                    @error "Error adding value to column :$colname."
                    rethrow(err)
                end
            else
                newcol = similar(col, promote_type(S, T), targetrows)
                copyto!(newcol, 1, col, 1, nrows)
                newcol[end] = val
                firstindex(newcol) != 1 && _onebased_check_error()
                _columns(ds)[i] = newcol
                setformat!(ds, colname => format_of_cur_col)
            end
        end
        for (colname, col) in zip(_names(ds), _columns(ds))
            if length(col) != targetrows
                for col2 in _columns(ds)
                    resize!(col2, nrows)
                end
                _attributes(ds).meta.modified[] = current_modified
                throw(AssertionError("Error adding value to column :$colname"))
            end
        end
        for colname in setdiff(keys(row), _names(ds))
            val = row[colname]
            S = typeof(val)
            if nrows == 0
                newcol = [val]
            else
                newcol = Tables.allocatecolumn(Union{Missing, S}, targetrows)
                fill!(newcol, missing)
                newcol[end] = val
            end
            ds[!, colname] = newcol
        end
        _modified(_attributes(ds))
        _reset_grouping_info!(ds)
        return ds
    end

    if cols == :orderequal
        if old_row_type <: Dict
            throw(ArgumentError("passing `Dict` as `row` when `cols == :orderequal` " *
                                "is not allowed as it is unordered"))
        elseif length(row) != ncol(ds) || any(x -> x[1] != x[2], zip(keys(row), _names(ds)))
            throw(ArgumentError("when `cols == :orderequal` pushed row must " *
                                "have the same column names and in the " *
                                "same order as the target data set"))
        end
    elseif cols === :setequal
        # Only check for equal lengths if :setequal is selected,
        # as an error will be thrown below if some names don't match
        if length(row) != ncols
            # an explicit error is thrown as this was allowed in the past
            throw(ArgumentError("`push!` with `cols` equal to `:setequal` " *
                                "requires `row` to have the same number of elements " *
                                "as the number of columns in `ds`."))
        end
    end
    current_col = 0
    current_modified = _attributes(ds).meta.modified[]
    try
        for (col, nm) in zip(_columns(ds), _names(ds))
            format_of_cur_col = getformat(ds, nm)
            current_col += 1
            if cols === :subset
                val = get(row, nm, missing)
            else
                val = row[nm]
            end
            S = typeof(val)
            T = eltype(col)
            if S <: T || !promote || promote_type(S, T) <: T
                push!(col, val)
            else
                newcol = similar(col, promote_type(S, T), targetrows)
                copyto!(newcol, 1, col, 1, nrows)
                newcol[end] = val
                firstindex(newcol) != 1 && _onebased_check_error()
                _columns(ds)[columnindex(ds, nm)] = newcol
                setformat!(ds, nm => format_of_cur_col)
            end
        end
        current_col = 0
        for col in _columns(ds)
            current_col += 1
            @assert length(col) == targetrows
        end
    catch err
        for col in _columns(ds)
            resize!(col, nrows)
        end
        _attributes(ds).meta.modified[] = current_modified
        @error "Error adding value to column :$(_names(ds)[current_col])."
        rethrow(err)
    end
    _modified(_attributes(ds))
    _reset_grouping_info!(ds)
    return ds
end

"""
    push!(ds::Dataset, row::Union{Tuple, AbstractArray}; promote::Bool=false)
    push!(ds::Dataset, row::Union{DatasetRow, NamedTuple, AbstractDict};
          cols::Symbol=:setequal, promote::Bool=(cols in [:union, :subset]))

Add in-place one row at the end of `ds` taking the values from `row`.

Column types of `ds` are preserved, and new values are converted if necessary.
An error is thrown if conversion fails.

If `row` is neither a `DatasetRow`, `NamedTuple` nor `AbstractDict` then
it must be a `Tuple` or an `AbstractArray`
and columns are matched by order of appearance. In this case `row` must contain
the same number of elements as the number of columns in `ds`.

If `row` is a `DatasetRow`, `NamedTuple` or `AbstractDict` then
values in `row` are matched to columns in `ds` based on names. The exact behavior
depends on the `cols` argument value in the following way:
* If `cols == :setequal` (this is the default)
  then `row` must contain exactly the same columns as `ds` (but possibly in a
  different order).
* If `cols == :orderequal` then `row` must contain the same columns in the same
  order (for `AbstractDict` this option requires that `keys(row)` matches
  `propertynames(ds)` to allow for support of ordered dicts; however, if `row`
  is a `Dict` an error is thrown as it is an unordered collection).
* If `cols == :intersect` then `row` may contain more columns than `ds`,
  but all column names that are present in `ds` must be present in `row` and only
  they are used to populate a new row in `ds`.
* If `cols == :subset` then `push!` behaves like for `:intersect` but if some
  column is missing in `row` then a `missing` value is pushed to `ds`.
* If `cols == :union` then columns missing in `ds` that are present in `row` are
  added to `ds` (using `missing` for existing rows) and a `missing` value is
  pushed to columns missing in `row` that are present in `ds`.

If `promote=true` and element type of a column present in `ds` does not allow
the type of a pushed argument then a new column with a promoted element type
allowing it is freshly allocated and stored in `ds`. If `promote=false` an error
is thrown.

As a special case, if `ds` has no columns and `row` is a `NamedTuple` or
`DatasetRow`, columns are created for all values in `row`, using their names
and order.

Please note that `push!` must not be used on a `Dataset` that contains columns
that are aliases (equal when compared with `===`).

# Examples
```jldoctest
julia> ds = Dataset(A=1:3, B=1:3);

julia> push!(ds, (true, false))
4×2 Dataset
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     1      0

julia> push!(ds, ds[1, :])
5×2 Dataset
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     1      0
   5 │     1      1

julia> push!(ds, (C="something", A=true, B=false), cols=:intersect)
6×2 Dataset
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      1
   2 │     2      2
   3 │     3      3
   4 │     1      0
   5 │     1      1
   6 │     1      0

julia> push!(ds, Dict(:A=>1.0, :C=>1.0), cols=:union)
7×3 Dataset
 Row │ A        B        C
     │ Float64  Int64?   Float64?
─────┼─────────────────────────────
   1 │     1.0        1  missing
   2 │     2.0        2  missing
   3 │     3.0        3  missing
   4 │     1.0        0  missing
   5 │     1.0        1  missing
   6 │     1.0        0  missing
   7 │     1.0  missing        1.0

julia> push!(ds, NamedTuple(), cols=:subset)
8×3 Dataset
 Row │ A          B        C
     │ Float64?   Int64?   Float64?
─────┼───────────────────────────────
   1 │       1.0        1  missing
   2 │       2.0        2  missing
   3 │       3.0        3  missing
   4 │       1.0        0  missing
   5 │       1.0        1  missing
   6 │       1.0        0  missing
   7 │       1.0  missing        1.0
   8 │ missing    missing  missing
```
"""
function Base.push!(ds::Dataset, row::Any; promote::Bool=false)

# Modify Dataset
    if !(row isa Union{Tuple, AbstractArray})
        # an explicit error is thrown as this was allowed in the past
        throw(ArgumentError("`push!` does not allow passing collections of type " *
                            "$(typeof(row)) to be pushed into a Dataset. Only " *
                            "`Tuple`, `AbstractArray`, `AbstractDict`, `DatasetRow` " *
                            "and `NamedTuple` are allowed."))
    end
    nrows, ncols = size(ds)
    targetrows = nrows + 1
    if length(row) != ncols
        msg = "Length of `row` does not match `Dataset` column count."
        throw(DimensionMismatch(msg))
    end
    current_col = 0
    current_modified = _attributes(ds).meta.modified[]
    try
        for (i, (col, val)) in enumerate(zip(_columns(ds), row))
            current_col += 1
            format_of_cur_col = getformat(ds, current_col)
            S = typeof(val)
            T = eltype(col)
            if S <: T || !promote || promote_type(S, T) <: T
                push!(col, val)
            else
                newcol = Tables.allocatecolumn(promote_type(S, T), targetrows)
                copyto!(newcol, 1, col, 1, nrows)
                newcol[end] = val
                firstindex(newcol) != 1 && _onebased_check_error()
                _columns(ds)[i] = newcol
                setformat!(ds, i => format_of_cur_col)
            end
        end
        current_col = 0
        for col in _columns(ds)
            current_col += 1
            @assert length(col) == targetrows
        end
    catch err
        #clean up partial row
        for col in _columns(ds)
            resize!(col, nrows)
        end
        _attributes(ds).meta.modified[] = current_modified
        @error "Error adding value to column :$(_names(ds)[current_col])."
        rethrow(err)
    end
    _modified(_attributes(ds))
    _reset_grouping_info!(ds)
    ds
end


text:
module PoissonBox

using TiledIteration
import LinearAlgebra
import IterativeSolvers

function CartesianIndices_inner(x::AbstractArray)
    ranges = map(axes(x)) do ax
        (first(ax)+1):(last(ax)-1)
    end
    CartesianIndices(ranges)
end

function forward!(out,u)
    if axes(out) != axes(u)
        msg = """
        Axes do not match:
        axes(out) == $(axes(out))
        axes(u) == $(axes(u))
        """
        throw(DimensionMismatch(msg))
    end
    all_inds = CartesianIndices(u)
    inner_inds = CartesianIndices_inner(u)
    begin
        for I in inner_inds
            out[I] = laplace(u, I)
        end
        for I in EdgeIterator(all_inds, inner_inds)
            out[I] = u[I]
        end
    end
    return out
end
function forward(u)
    out = similar(u, float(eltype(u)))
    forward!(out, u)
end

Base.@propagate_inbounds function laplace(u::AbstractArray{<:Any,1},I::CartesianIndex)
    (i,) = Tuple(I)
    (u[i+1] + u[i-1])/2 - u[i]
end

Base.@propagate_inbounds function laplace(u::AbstractArray{<:Any,2},I::CartesianIndex)
    (i,j) = Tuple(I)
    (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i, j-1])/4 - u[i,j]
end
Base.@propagate_inbounds function laplace(u::AbstractArray{<:Any,3},I::CartesianIndex)
    (i,j,k) = Tuple(I)
    (u[i+1,j  ,k  ] + u[i-1,j   ,k  ] +
     u[i  ,j+1,k  ] + u[i  , j-1,k  ] +
     u[i  ,j  ,k+1] + u[i  , j  ,k-1]
    )/6 - u[i,j,k]
end

################################################################################
#### LaplaceDirichlet
################################################################################
struct LaplaceDirichlet{T,N}
    size::NTuple{N,Int}
end

make_eltype_float(x::AbstractArray{<:AbstractFloat}) = x
function make_eltype_float(x::AbstractArray)
    T = float(eltype(x))
    out = similar(x, T)
    copy!(out, x)
end

function poisson_problem(arr::AbstractArray{T,N}) where {T,N}
    op = LaplaceDirichlet{float(T),N}(size(arr))
    v = reshape(make_eltype_float(arr), length(arr))
    op, v
end
function reshape_solution(op, sol)
    reshape(sol, op.size)
end

function LinearAlgebra.mul!(out, op::LaplaceDirichlet, v::AbstractVector)
    forward!(reshape(out, op.size), reshape(v, op.size))
    out
end
function Base.:(*)(op::LaplaceDirichlet, v::AbstractVector)
    T = typeof(one(eltype(v))/one(eltype(op)))
    out = similar(v,T)
    mul!(out, op, v)
end
function Base.eltype(::Type{LaplaceDirichlet{T,N}}) where {T,N}
    T
end
function Base.size(op::LaplaceDirichlet, d)
    (1,2)[d] # for appropriate error
    prod(op.size)
end

function Base.collect(::Type{T}, op::LaplaceDirichlet) where {T}
    # useful for debugging
    n = size(op,1)
    out = Matrix{T}(undef, n,n)
    v = zeros(Int,n)
    for i in 1:n
        v[i] = 1
        if i > 1
            v[i-1] = 0
        end
        mul!(view(out,:,i), op, v)
    end
    out
end
function Base.collect(op::LaplaceDirichlet)
    Base.collect(eltype(op), op)
end

function solve(solver, y::AbstractArray, solver_args...; solver_kw...)
    op, v = poisson_problem(y)
    sol = solver(op, v, solver_args...; solver_kw...)
    reshape_solution(op, sol)
end
function solve(y::AbstractArray)
    solve(IterativeSolvers.bicgstabl, y, 2)
end

end


text:
"""
Submodule of `SirenSeq`.
Used for rendering `SirenSeq` expression to a `.pdf` file.
Must have `musescore` 2.0 installed.
"""
module Render

export renderMidi

using SirenSeq
using SirenSeq.defaultBpm

"""
	renderMidi(path=\"temp\", name=path)

Renders \"`path`.mid\" to \"`name`.pdf\" using musescore. \\
Make sure musescore 2.0 or higher is installed.
"""
function renderMidi(; path::AbstractString="temp", name::AbstractString=path)
	print("rendering: $(name).pdf ... \t")
	try
		run(`mscore -o $(name).pdf $(path).mid`)
		print("done\n")
	catch err
		print("failed!\n")
		throw(err)
	end
	return
end

"""
	renderMidi(x::Expi; path=\"temp\", name=path, tpq=96, bpm=defaultBpm)

Creates a midi file caled \"`path`.mid\" (do not include the \".mid\" part) using audio expression `x`. \\
`tpq` is the midi ticks per quarter note.  `bpm` is the intial beats per minute (may be changed within `x` at any time). \\
Renders \"`path`.mid\" to \"`name`.pdf\" using musescore. \\
Make sure musescore 2.0 or higher is installed.
"""
function renderMidi(x; path::AbstractString="temp", name=path, tpq::Int=96, bpm::Int=defaultBpm)
	makeMidi(x;path=path,tpq=tpq,bpm=bpm)
	sleep(0.2)
	renderMidi(path=path,name=name)
	return
end

end


text:
"""
    Rimu
Random Integrator for Many-Body Quantum Systems
"""
module Rimu

using Arrow
using DataFrames
using DataStructures
using LinearAlgebra
using OrderedCollections # for LittleDict
using Parameters
using Reexport
using Setfield
using SplittablesBase
using StaticArrays
using StatsBase
using ThreadsX

@reexport using Distributed

include("helpers.jl") # non MPI-dependent helper functions

include("ConsistentRNG.jl")
@reexport using .ConsistentRNG
include("Interfaces/Interfaces.jl")
@reexport using .Interfaces
include("StochasticStyles/StochasticStyles.jl")
@reexport using .StochasticStyles
include("DictVectors/DictVectors.jl")
@reexport using .DictVectors
include("BitStringAddresses/BitStringAddresses.jl")
@reexport using .BitStringAddresses
include("Hamiltonians/Hamiltonians.jl")
@reexport using .Hamiltonians
include("RimuIO.jl")
@reexport using .RimuIO
include("StatsTools/StatsTools.jl")
@reexport using .StatsTools

export lomc!
export FciqmcRunStrategy, RunTillLastStep
export MemoryStrategy, NoMemory, DeltaMemory, ShiftMemory
export ShiftStrategy, LogUpdate, LogUpdateAfterTargetWalkers
export DontUpdate, DoubleLogUpdate, DoubleLogUpdateAfterTargetWalkers
export ReportingStrategy, ReportDFAndInfo, ReportToFile
export ReplicaStrategy, NoStats, AllOverlaps
export PostStepStrategy, Projector, ProjectedEnergy, SignCoherence, WalkerLoneliness, Timer
export TimeStepStrategy, ConstantTimeStep, OvershootControl
export threadedWorkingMemory, localpart, walkernumber

include("strategies_and_params/fciqmcrunstrategy.jl")
include("strategies_and_params/memorystrategy.jl")
include("strategies_and_params/poststepstrategy.jl")
include("strategies_and_params/replicastrategy.jl")
include("strategies_and_params/reportingstrategy.jl")
include("strategies_and_params/shiftstrategy.jl")
include("strategies_and_params/timestepstrategy.jl")
include("strategies_and_params/threadingstrategy.jl")
include("strategies_and_params/deprecated.jl")

include("apply_memory_noise.jl")
include("lomc.jl")                  # top level

include("RMPI/RMPI.jl")

# Modules for parallel computing not exported by default for now
include("EmbarrassinglyDistributed.jl")

# analysis tool not reexported (to be deprecated)
include("Blocking.jl")

end # module


text:
using BinDeps

@BinDeps.setup

fluidsynth = library_dependency("fluidsynth", aliases=["libfluidsynth"])

provides(AptGet, "libfluidsynth-dev", fluidsynth)

@static if Sys.isapple()
    using Homebrew
    provides(Homebrew.HB, "fluid-synth", fluidsynth, os = :Darwin)
end

@BinDeps.install Dict(:fluidsynth => :libfluidsynth)


text:
using Compose
using Test
using EzXML
using Colors
using Measures

@testset "Issue 267" begin
    global c = compose(context(), fill(["red", "blue"]),
                [context(), fill("green"), circle([0.25, 0.75], [0.5], [0.25])])
    img = SVG(8cm, 6cm, false)
    draw(img, c)
    svgxml = root(parsexml(String(img.out.data)))
    # get all grouped values that have the fill attribute
    fillcolors = nodecontent.(findall("//ns:g[@fill]/@fill", svgxml, ["ns"=>namespace(svgxml)]))
    # there should only be a single color (because green should clobber the red
    # and blue colors)
    @test length(fillcolors) == 1
    # make sure it's green
    @test fillcolors[1] == "#008000"
end


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets, ScientificTypes, PrettyPrinting

@load LinearRegressor pkg=MLJLinearModels
@load RidgeRegressor pkg=MLJLinearModels
@load LassoRegressor pkg=MLJLinearModels

hitters = dataset("ISLR", "Hitters")
@show size(hitters)
names(hitters) |> pprint

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
# This file is a part of RunStatistics.jl, licensed under the MIT License (MIT).

using RunStatistics
using Test

@testset "squares_approx" begin

    T_obs = 20
    N = 30
    n = 100
    epsp = 10^(-2)

    @test squares_pvalue_approx(T_obs, N, n, epsp) ≈ 0.24164150728705625
    @test squares_cdf_approx(T_obs, N, n, epsp) ≈ 0.7583584927129438
end

# I measured the timing for N=96 with the exact and approximate solution. See if julia is really faster. It seems some experiment with threading was done but without threads, I'd be surprised if it's faster than the C++ code built in release mode


text:
module RefHull
export verify_hull

using ConicHulls.Dets.incomplete_det

function eliminate_duplicates{T<:Integer}(generators::Vector{Vector{T}})
    found = Set{NTuple{length(generators[1]), T}}()
    gs = Vector{T}[]
    for g in generators
        gd = div(g, gcd(g...))
        key = tuple(gd...)
        if !(key in found)
            push!(gs, g)
            push!(found, key)
        end
    end
    gs
end

type SortedFacet{T}
    generators::Vector{Vector{T}}
    positive::Bool
    plane::Vector{T}
end

function enumerate_planes{T}(NC::Int, generators::Vector{Vector{T}})
    NF = NC - 1
    extreme = ObjectIdDict()
    facets = SortedFacet{T}[]
    for ks in combinations(1:length(generators), NF)
        gs = generators[ks]
        plane = incomplete_det(gs...)

        if all(plane .== 0); continue; end

        positive = negative = npos = nneg = 0
        for (k,g) in enumerate(generators)
            if k in ks; continue; end

            d = dot(plane, g)
            if d > 0;     positive = k; npos += 1;
            elseif d < 0; negative = k; nneg += 1;
            end
            if (npos > 1) && (nneg > 1); break; end
        end
        if npos == 1; extreme[generators[positive]] = true; end
        if nneg == 1; extreme[generators[negative]] = true; end
        if npos == nneg == 0
            error("Generators do not span the the space")
        end
        if npos == 0 || nneg == 0
            positive = (npos == 0)
            push!(facets, SortedFacet(gs, positive, positive ? plane : -plane))
        end
    end
    # Filter out non-extreme generators
    gs = Vector{T}[]
    for g in generators
        if haskey(extreme, g); push!(gs, g); end
    end
    # Filter out facets with non-extreme generators
    fs = SortedFacet{T}[]
    for facet in facets
        if all([haskey(extreme, g) for g in facet.generators])
            push!(fs, facet)
        end
    end
    gs, fs
end

function filter_coplanar_facets{T}(generators::Vector{Vector{T}}, facets::Vector{SortedFacet{T}})
    NC = length(generators[1])
    order = ObjectIdDict()
    for (k,g) in enumerate(generators); order[g] = k; end

    fs = SortedFacet{T}[]
    for facet in facets
        j = first(find(facet.plane .!= 0))
        dominated = false
        for g in generators
            if g in facet.generators; continue; end

            d = dot(facet.plane,g)
            @assert d <= 0
            if d < 0; continue; end

            # Coplanar
            gs = Vector{T}[facet.generators..., g]

            ks = Int[order[g] for g in facet.generators]
            push!(ks, order[g])

            for i in sortperm(ks, rev=true)
                p = incomplete_det(gs[1:(i-1)]...,gs[(i+1):end]...)
                if all(p .== 0); continue; end

                sameway = (p[j] > 0) == (facet.plane[j] > 0)
                if (sameway == iseven(NC-i)) == facet.positive
                    dominated = true;
                else
                    @assert i < NC # i == NC ==> should always be dominated
                end
                break
            end
            if dominated; break; end
        end
        if !dominated; push!(fs, facet); end
    end
    fs
end

function conichull{T<:Integer}(NC::Int, generators::Vector{Vector{T}})
    # Verify input data
    for g in generators; @assert length(g) == NC; end

    generators         = eliminate_duplicates(generators)
    generators, facets = enumerate_planes(NC, generators)
    facets             = filter_coplanar_facets(generators, facets)

    generators, facets
end

using ConicHulls: Common, Hulls

function verify_hull(hull::ConicHull, generators::Vector)
    orig_generators = generators
    validate(hull)

    gmap = ObjectIdDict()
    for g in generators; gmap[g.x] = g; end
    generators = [g.x for g in generators]
    generators, facets = RefHull.conichull(nconic(hull), generators)

    hfacets = {}
    for facet in hull.facets
        gs, even = get_canonical_winding(facet)
        push!(hfacets, (gs, even))
    end
    sort!(hfacets,
          lt=(w1,w2)->(
              for (g1,g2) in zip(w1[1], w2[1]);
                  if g1.id < g2.id; return true; end;
                  if g2.id > g2.id; return false; end;
              end;
              return false;))

    try
        @assert length(hull.facets) == length(facets)
        facetset = Set()
        for facet in facets
            push!(facetset, tuple([gmap[g] for g in facet.generators]..., facet.positive))
        end
        for facet in hull.facets
            gs, even = get_canonical_winding(facet)
            key = tuple(gs..., even)
            @assert key in facetset
        end
    catch e
        @show length(hull.facets), length(facets)

        println("hull.facets:")
#        for facet in hull.facets; println("  ", facet.generators); end
        for (gens, pos) in hfacets; println("  ", gens, ": ", pos); end
        println("facets:")
        for facet in facets; println("  ", [gmap[g] for g in facet.generators], ": ", facet.positive); end
        println("generators:")
        for g in orig_generators; println("  ", g, ": ", g.x); end
        println("filtered generators:")
        for g in generators; println("  ", gmap[g]); end

        rethrow(e)
    end
end

end # module


text:
# compact end efficient sampler for 2D distribution

struct MatrixAliasSampler{IndexType<:Integer, ProbType<:Real}
  acceptances_mat::Matrix{ProbType}
  aliases_mat::Matrix{IndexType}
  # each column of both matrices represent an alias sampler
end

function MatrixAliasSampler(probs::AbstractMatrix{T}, IdxType::Type=Int) where T<:Real
  n = size(probs, 1)
  @assert size(probs, 2) == n
  aliases_mat = zeros(IdxType, n, n)
  acceptances_mat = zeros(T, n, n)

  smalls = Vector{IdxType}(undef, n)
  larges = Vector{IdxType}(undef, n)

  for i in 1:n
    weights = view(probs, :, i)
    aliases = view(aliases_mat, :, i)
    acceptances = view(acceptances_mat, :, i)

    setup_alias_sampler!(weights, acceptances, aliases, smalls, larges, sum(weights))
  end
  MatrixAliasSampler{IdxType, T}(acceptances_mat, aliases_mat)
  end

function sample(rng::AbstractRNG, m::MatrixAliasSampler, source::Integer)
  N, M = m.aliases_mat |> size

  @assert N==M
  @assert N == size(m.acceptances_mat,1) == size(m.acceptances_mat, 2)
  acceptances = @view m.acceptances_mat[:, source]
  aliases = @view m.aliases_mat[:, source]

  Idx = eltype(aliases)
  idx = rand(rng, Idx(1):Idx(N))
  if acceptances[idx] >= 1.0
    return idx
  end
  rand(rng) < acceptances[idx] ? idx : aliases[idx]
end


text:
function ABCDQR_simula(nt::Int,seed::Int)
	#
	# simulated model to test the functions
	#
	# nt: number of time instants
	# seed: seed for the random numbers

	# random seed
	# ------------------------------
	Random.seed!(seed)

	# simulated system
	m = simula01()

	# vibration parameters
	# ----------------------------
	# time interval (s)
	dt = m["dt"]

	# state-space matrices
	# --------------------------------------
	A = m["A"]
	B = m["B"]
	C = m["C"]
	D = m["D"]

	# input definition
	# ---------------------------------------
	# system input
	u = 100*randn(2,nt) # if 100 is not used, Q and R are too little

	# state equation, no noise
	# ---------------------------------
	x = zeros(4,nt)
	for t=1:nt-1
		x[:,t+1] = A*x[:,t] + B*u[:,t]
	end
	P = cov(x')
	Q = 0.01*P # std approx 0.1
	ch = cholesky(Q)
	w = ch.L*randn(4,nt)

	# state equation with noise
	# ------------------------------------------
	x = zeros(4,nt)
	for t = 1:nt-1
		x[:,t+1] = A*x[:,t] + B*u[:,t] + w[:,t]
	end

	# observation equation
	# -------------------------------------
	y = C*x + D*u
	L = cov(y')
	R = 0.01*L # std approx 0.1
	ch = cholesky(R)
	v = ch.L*randn(2,nt)
	y = y + v

	#
	S = cov(w',v')

	# starting values
	m1 = C\y[:,1:1]
	P1 = 0*A

	m = Dict("dt"=>dt,
	"y"=>y,"u"=>u,
	"A"=>A,"B"=>B,"C"=>C,"D"=>D,
	"Q"=>Q,"R"=>R,"S"=>S,
	"m1"=>m1,"P1"=>P1)

	return m

end


text:
module AtomGraphs

using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using Serialization
using Xtals
#rc[:paths][:crystals] = @__DIR__ # so that Xtals.jl knows where things are
using MolecularGraph

import ChemistryFeaturization: elements

include("atomgraph.jl")
export AtomGraph, elements, visualize

include("graph_building.jl")
export inverse_square, exp_decay

end


text:
# module operators

export ⊕, ⊗

###############################
# operators
###############################

#=
[1] (exponentiation)
^ ↑ ↓ ⇵ ⟰ ⟱ ⤈ ⤉ ⤊ ⤋ ⤒ ⤓ ⥉ ⥌ ⥍ ⥏ ⥑ ⥔ ⥕ ⥘ ⥙ ⥜ ⥝ ⥠ ⥡ ⥣ ⥥ ⥮ ⥯ ￪ ￬

[2] (unary)
+ - ! ~ ¬ √ ∛ ∜ ⋆ ± ∓

[3] (bitshift)
<< >> >>>

[4] (fraction)
//

[5] (multiplication)
* / ÷ % & ⋅ ∘ × |\\| ∩ ∧ ⊗ ⊘ ⊙ ⊚ ⊛ ⊠ ⊡ ⊓ ∗ ∙ ∤ ⅋ ≀ ⊼ ⋄ ⋆ ⋇ ⋉ ⋊ ⋋ ⋌ ⋏ ⋒ ⟑ ⦸ ⦼ ⦾ ⦿ ⧶ ⧷ ⨇ ⨰ ⨱ ⨲ ⨳ ⨴ ⨵ ⨶ ⨷ ⨸ ⨻ ⨼ ⨽ ⩀ ⩃ ⩄ ⩋ ⩍ ⩎ ⩑ ⩓ ⩕ ⩘ ⩚ ⩜ ⩞ ⩟ ⩠ ⫛ ⊍ ▷ ⨝ ⟕ ⟖ ⟗ ⨟

[6] (addition)
+ - |\|| ⊕ ⊖ ⊞ ⊟ |++| ∪ ∨ ⊔ ± ∓ ∔ ∸ ≏ ⊎ ⊻ ⊽ ⋎ ⋓ ⧺ ⧻ ⨈ ⨢ ⨣ ⨤ ⨥ ⨦ ⨧ ⨨ ⨩ ⨪ ⨫ ⨬ ⨭ ⨮ ⨹ ⨺ ⩁ ⩂ ⩅ ⩊ ⩌ ⩏ ⩐ ⩒ ⩔ ⩖ ⩗ ⩛ ⩝ ⩡ ⩢ ⩣

=#

# [1]
# f^x == f(x)
Base.:^(f, x) = x | functionalize(f)

# [2]
# ~f == functionalize(f)
# -f == a->functionalize(f)(a...)
# +f == a->functionalize(f).(a)
Base.:~(f) = functionalize(f)
Base.:-(f) = a -> functionalize(f)(a...)
Base.:+(f) = a -> functionalize(f).(a)

# [3]
# f // x == f(x)
Base.://(f, x) = x | functionalize(f)

# [4]
# f & g == f / g == x -> g(f(x))
Base.:&(f, g) = functionalize(g) ∘ functionalize(f)
Base.:/(f, g) = functionalize(g) ∘ functionalize(f)

# [5]
# x | f == f(x)
Base.:|(x, f) = (g = functionalize(f); Core.applicable(g, x) ? g(x) : g(x...))

# args(x...;ka...) | f == f(x...;ka...)
Base.:|(x::_Args, f) = x(f)

to_rng(i::Integer) = i>=1 ? (1:i) : (-1:-1:i)
to_rng(x) = x

###############################
# zip
###############################

@inline ⊕(x, y) = zip(to_rng(x), to_rng(y))
@inline ⊕(x::Base.Iterators.Zip, y) = zip(x.is..., to_rng(y))
@inline ⊕(x::Base.Iterators.Zip, y::Base.Iterators.Zip) = zip(x, y)
@inline ⊕(x...) = zip(x...)

###############################
# product
###############################

@inline ⊗(x, y) = Iterators.product(to_rng(x), to_rng(y))
@inline ⊗(x::Base.Iterators.ProductIterator, y) =
    Iterators.product(x.iterators..., to_rng(y))
@inline ⊗(x::Base.Iterators.ProductIterator, y::Base.Iterators.ProductIterator) =
    Iterators.product(x, y)
@inline ⊗(x...) = Iterators.product(x...)

#=
###############################
# concatenate
###############################

@inline ⊞(x) = ⊞(x...)
@inline ⊞(x::Tuple, y::Tuple) = tpljoin(x,y)
@inline ⊞(x::AbstractArray, y::AbstractArray) = aryjoin(x,y)

###############################
# fgen, fary, tfary, pfary
###############################

@inline ⩅(itr, fnc) = fgen(fnc, itr)
# @inline ⩂(fnc, itr) = fgen(fnc, itr)
@inline ⋓(itr, fnc) = fary(fnc, itr)
@inline ⩐(itr, fnc) = fary(fnc, itr; prgrs=true)
@inline ⩏(itr, fnc) = tfary(fnc, itr)
@inline ⊔(itr, fnc) = tfary(fnc, itr; prgrs=true)
@inline ⩖(itr, fnc) = pfary(fnc, itr)
@inline ⩒(itr, fnc) = pfary(fnc, itr; prgrs=true)
=#

# end