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

name = "Librsvg"
version = v"2.42.2"

# Collection of sources required to build librsvg
sources = [
    ArchiveSource("https://download.gnome.org/sources/librsvg/$(version.major).$(version.minor)/librsvg-$(version).tar.xz",
                  "0c550a0bffef768a436286116c03d9f6cd3f97f5021c13e7f093b550fac12562"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/librsvg-*/

# Pango's `.la` files are wrong on certain platforms, and we don't need them anyway
rm -f ${prefix}/lib/*.la

atomic_patch -p1 "${WORKSPACE}/srcdir/patches/librsvg_link_order.patch"
autoreconf -f -i

# We need this for bootstrapping purposes
apk add gdk-pixbuf

FLAGS=()
if [[ "${target}" == *-apple-* ]]; then
    # We purposefully use an old binutils, so we must disable -Bsymbolic
    FLAGS+=(--disable-Bsymbolic)
fi

# cssparser must be upgraded as it doesn't build properly anymore
sed -i.bak -e 's&cssparser = "0.23"&cssparser = "0.25"&' rust/Cargo.toml
(cd rust && cargo vendor)

LDFLAGS="-L${prefix}/lib -L${prefix}/lib64 -Wl,-rpath,${prefix}/lib -Wl,-rpath,${prefix}/lib64" ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --disable-static \
    --enable-pixbuf-loader \
    --disable-introspection \
    --disable-gtk-doc-html \
    --enable-shared \
    "${FLAGS[@]}"

if [[ ${target} == *mingw* ]]; then
    # pass static rust package to linker
    sed -i "s/^deplibs_check_method=.*/deplibs_check_method=\"pass_all\"/g" libtool
    # add missing crt libs (ws2_32 and userenv) to LIBRSVG_LIBS
    sed -i "s/^LIBRSVG_LIBS = .*/& -lws2_32 -luserenv/g" Makefile
fi

# Don't try to unwind on i686-w64-mingw32, just panic because rust doesn't know how to SLJL
# https://github.com/rust-lang/rust/issues/12859#issuecomment-185081071
if [[ ${target} == i686-w64-mingw32 ]]; then
    export RUSTFLAGS="-C panic=abort"
fi

# Manually build rust library because rust target doesn't match C target
(cd rust && PKG_CONFIG_ALLOW_CROSS=1 cargo build --release)

RUST_LIB="$(pwd)/rust/target/${rust_target}/release/librsvg_internals.a"
make RUST_LIB="${RUST_LIB}" -j${nproc}
make RUST_LIB="${RUST_LIB}" install
"""

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

# The products that we will ensure are always built
products = [
    # The main event
    LibraryProduct("librsvg-2", :librsvg),

    # This is named `.so` even on darwin, so do it as a FileProduct.....sigh
    FileProduct(["lib/gdk-pixbuf-2.0/2.10.0/loaders/libpixbufloader-svg.so",
                 "lib/gdk-pixbuf-2.0/2.10.0/loaders/libpixbufloader-svg.dll"], :libpixbufloader_svg),
    #LibraryProduct("libpixbufloader-svg", :libpixbufloader_svg, ["lib/gdk-pixbuf-2.0/2.10.0/loaders"]),

    # And to round it out, let's get an executable as well!
    ExecutableProduct("rsvg-convert", :rsvg_convert),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("gdk_pixbuf_jll"),
    Dependency("Pango_jll", v"1.42.4"; compat="1.42.4"),
    Dependency("Libcroco_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :rust])

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
# SET MODEL PARAMETERS

#*========== MAIN PARAMETERS ===========
# beyond what fraction of line limit do we penalize?
const delta_line = 0.95

# how much costlier are the regular generator reserves compared to normal costs?
const reservesCostFactor = 1.5

# penalty parameter for exceeding reserves in the recourse model
const gamma_res = 10.0

# how do the wind generator reserve costs scale compared to the cheapest regular generator?
const windReservesCostFactor = 0.1

#*========== OTHER PARAMETERS ===========
# parameter for smoothing saturation function (relative to bounds)
const tau_sat = 1E-04

# parameter for smooth positive function approximation
const tau_pos = 1E-04

# tolerance for solving recourse system of nonlinear equations
const slackTolerance = 1E-06

# display precision
const displayPrecision = 4

# factor by which to scale wind forecast standard deviation to compute first-stage wind generator upper bounds (pmax_FS)
const windBoundFactor = 5.0

ext: jl

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
# Normal evaluation and normalization

#################################################
#
#   Vector norms
#
#################################################

function vnorm(x::NumericArray, p::Real)
    p > 0 || throw(ArgumentError("p must be positive."))
    p == 1 ? sumabs(x) :
    p == 2 ? sqrt(sumsq(x)) :
    isinf(p) ? maxabs(x) :
    sum(FixAbsPow(p), x) .^ inv(p)
end

vnorm(x::NumericArray) = vnorm(x, 2)

function vnorm!(dst::NumericArray, x::NumericArray, p::Real, dims::DimSpec)
    fill!(dst, zero(eltype(dst)))
    p > 0 || throw(ArgumentError("p must be positive."))
    p == 1 ? sumabs!(dst, x, dims) :
    p == 2 ? map1!(SqrtFun(), sumsq!(dst, x, dims)) :
    isinf(p) ? maxabs!(dst, x, dims) :
    map1!(FixAbsPow(inv(p)), sum!(dst, FixAbsPow(p), x, dims))
end

function vnorm(x::NumericArray, p::Real, dims::DimSpec)
    tt = promote_type(fptype(eltype(x)), fptype(typeof(p)))
    r = Array(tt, Base.reduced_dims(size(x), dims))
    vnorm!(r, x, p, dims)
end

# vnormdiff

function vnormdiff(x::DenseArrOrNum, y::DenseArrOrNum, p::Real)
    p > 0 || throw(ArgumentError("p must be positive."))
    p == 1 ? sumabsdiff(x, y) :
    p == 2 ? sqrt(sumsqdiff(x, y)) :
    isinf(p) ? maxabsdiff(x, y) :
    sumfdiff(FixAbsPow(p), x, y) .^ inv(p)
end

vnormdiff(x::DenseArrOrNum, y::DenseArrOrNum) = vnormdiff(x, y, 2)

function vnormdiff!(dst::NumericArray, x::DenseArrOrNum, y::DenseArrOrNum, p::Real, dims::DimSpec)
    fill!(dst, zero(eltype(dst)))
    p > 0 || throw(ArgumentError("p must be positive."))
    p == 1 ? sumabsdiff!(dst, x, y, dims) :
    p == 2 ? map1!(SqrtFun(), sumsqdiff!(dst, x, y, dims)) :
    isinf(p) ? maxabsdiff!(dst, x, y, dims) :
    map1!(FixAbsPow(inv(p)), sumfdiff!(dst, FixAbsPow(p), x, y, dims))
end

function vnormdiff(x::DenseArrOrNum, y::DenseArrOrNum, p::Real, dims::DimSpec)
    tt = promote_type(fptype(promote_type(eltype(x), eltype(y))), fptype(typeof(p)))
    r = Array(tt, Base.reduced_dims(size(x), dims))
    vnormdiff!(r, x, y, p, dims)
end

#################################################
#
#   Normalization
#
#################################################

# normalize an array as a whole

function normalize!(dst::ContiguousRealArray, x::ContiguousRealArray, p::Real)
    return map!(Multiply(), dst, x, inv(vnorm(x, p)))
end
normalize!(x::ContiguousRealArray, p::Real) = normalize!(x, x, p)
normalize(x::ContiguousRealArray, p::Real) = x * inv(vnorm(x, p))

# normalize along specific dimension

function normalize!(dst::NumericArray, x::NumericArray, p::Real, dims::DimSpec)
    p > 0 || throw(ArgumentError("p must be positive."))
    length(dst) == length(x) || throw(ArgumentError("Inconsistent argument dimensions!"))

    if length(dims) == 1 && dims[1] == 1
        n = succ_length(size(x), 1)
        if p == 1
            for j = 1:n
                xj = view(x,:,j)
                map!(Multiply(), view(dst,:,j), xj, inv(sumabs(xj)))
            end
        elseif p == 2
            for j = 1:n
                xj = view(x,:,j)
                map!(Multiply(), view(dst,:,j), xj, inv(sqrt(sumsq(xj))))
            end
        elseif isinf(p)
            for j = 1:n
                xj = view(x,:,j)
                map!(Multiply(), view(dst,:,j), xj, inv(maxabs(xj)))
            end
        else
            for j = 1:n
                xj = view(x,:,j)
                u = sum(FixAbsPow(p), xj) .^ inv(p)
                map!(Multiply(), view(dst,:,j), xj, inv(u))
            end
        end

    elseif length(dims) == 2 && dims[1] == 1 && dims[2] == 2
        n = succ_length(size(x), 2)
        if p == 1
            for j = 1:n
                xj = view(x,:,:,j)
                map!(Multiply(), view(dst,:,:,j), xj, inv(sumabs(xj)))
            end
        elseif p == 2
            for j = 1:n
                xj = view(x,:,:,j)
                map!(Multiply(), view(dst,:,:,j), xj, inv(sqrt(sumsq(xj))))
            end
        elseif isinf(p)
            for j = 1:n
                xj = view(x,:,:,j)
                map!(Multiply(), view(dst,:,:,j), xj, inv(maxabs(xj)))
            end
        else
            for j = 1:n
                xj = view(x,:,:,j)
                u = sum(FixAbsPow(p), xj) .^ inv(p)
                map!(Multiply(), view(dst,:,:,j), xj, inv(u))
            end
        end

    else
        broadcast!(.*, dst, x, rcp!(vnorm(x, p, dims)))
    end
    dst
end

normalize!(x::ContiguousRealArray, p::Real, dims::DimSpec) = normalize!(x, x, p, dims)
normalize(x::ContiguousRealArray, p::Real, dims::DimSpec) = normalize!(similar(x), x, p, dims)


text:
using FractionalTransforms
using Test

@testset "Test FRCT" begin
    @test isapprox(real(frct([1,2,3], 0.5, 2)), [1.707106781186547, 1.5606601717798205, -0.3535533905932727], atol=1e-5)
    @test isapprox(imag(frct([1,2,3], 0.5, 2)), [0.9874368670764581, -1.3964466094067267, -0.6982233047033652], atol=1e-5)
end


text:
# Return codes adapted from papi.h
@enum(RetCode,
      OK = 0,
      EINVAL = -1,      # Invalid argument
      ENOMEM = -2,      # Insufficient memory
      ESYS   = -3,      # A System C library call failed
      ECMP   = -4,      # Not supported by component
      ECLOST = -5,      # Access to the counters was lost or interrupted
      EBUG   = -6,      # Internal error, please send mail to the developers
      ENOEVNT = -7,     # Event does not exist
      ECNFLCT = -8,     # Event exists, but cannot be counted due to counter resource limitations
      ENOTRUN = -9,     # EventSet is currently not running
      EISRUN  = -10,    # EventSet is currently counting
      ENOEVST = -11,    # No such EventSet Available
      ENOTPRESET = -12, # Event in argument is not a valid preset
      ENOCNTR = -13,    # Hardware does not support performance counters
      EMISC = -14,      # Unknown error code
      EPERM = -15,      # Permission level does not permit operation
      ENOINIT = -16,    # PAPI hasn't been initialized yet
      ENOCMP = -17,     # Component Index isn't set
      ENOSUPP = -18,    # Not supported
      ENOIMPL = -19,    # Not implemented
      EBUF = -20,       # Buffer size exceeded
      EINVAL_DOM = -21, # EventSet domain is not supported for the operation
      EATTR = -22,      # Invalid or missing event attributes
      ECOUNT = -23,     # Too many events or attributes */attributes
      ECOMBO =  -24,    # Bad combination of  features
)

const errmsg = let msgs = Dict{RetCode,ASCIIString}(
        OK => "OK",
        EINVAL => "Invalid argument",
        ENOMEM => "Insufficient memory",
        ESYS => "A System C library call failed",
        ECMP => "Not supported by component",
        ECLOST => "Access to the counters was lost or interrupted",
        EBUG => "Internal error, please send mail to the developers",
        ENOEVNT => "Event does not exist",
        ECNFLCT => "Event exists, but cannot be counted due to counter resource limitations",
        ENOTRUN => "EventSet is currently not running",
        EISRUN => "EventSet is currently counting",
        ENOEVST => "No such EventSet Available",
        ENOTPRESET => "Event in argument is not a valid preset",
        ENOCNTR => "Hardware does not support performance counters",
        EMISC => "Unknown error code",
        EPERM => "Permission level does not permit operation",
        ENOINIT => "PAPI hasn't been initialized yet",
        ENOCMP => "Component Index isn't set",
        ENOSUPP => "Not supported",
        ENOIMPL => "Not implemented",
        EBUF => "Buffer size exceeded",
        EINVAL_DOM => "EventSet domain is not supported for the operation",
        EATTR => "Invalid or missing event attributes",
        ECOUNT => "Too many events or attributes",
        ECOMBO => "Bad combination of  features"
    )
    errmsg(r::RetCode) = msgs[r]
end


text:
function augment_with_forward_solution(
    estimates;
    solver,
    control_system,
    player_cost_models_gt,
    x0,
    T,
    verbose = false,
    kwargs...,
)
    @showprogress pmap(enumerate(estimates)) do (ii, estimate)
        # overwrite the weights of the ground truth model with the weights of the estimate.
        player_cost_models_est =
            map(player_cost_models_gt, estimate.player_weights) do cost_model_gt, weights
                merge(cost_model_gt, (; weights))
            end

        # solve the forward game at this point
        converged, forward_solution, forward_opt_model =
            solve_game(solver, control_system, player_cost_models_est, x0, T; kwargs...)

        converged ||
            verbose && @warn "Forward solution augmentation did not converge on observation $ii."

        merge(estimate, (; forward_solution.x, forward_solution.u, converged))
    end
end


text:
function sparse(X::NTuple{N,InvOrAdjOrVecc}) where {N}
    # only works if there are sparse methods for each Vecchia Factor in X
	foldr(*, map(sparse, X))
end

function Matrix(X::NTuple{N,InvOrAdjOrVecc}) where {N}
	foldr(*, map(Matrix, X))
end

# sparse, Matrix and show for Ridiagonal
# =================================================

# ## sparse
function sparse(R::Ridiagonal{T}) where {T}
    n           = size(R,1)
    block_sizes = block_size(R,1)
    Rmat = PseudoBlockArray(
        spdiagm(fill(one(T),n)),
        block_sizes,
        block_sizes
    )
    for i=1:length(R.data)
        Rmat[Block(i+1,i)] .= R.data[i]
    end
    sparse(Rmat)
end

function sparse(adjR::Adjoint{<:Any,<:Ridiagonal})
	adjoint(sparse(adjR.parent))
end

# No sparse representations for invR or inv_adj_R
# sparse(invR::Inv{<:Any,<:Ridiagonal})
# sparse(inv_adj_R::Inv{<:Any,<:Adjoint{<:Any,<:Ridiagonal}})

# ## Matrix
Matrix(R::Ridiagonal) = Matrix(sparse(R))

function Matrix(adjR::Adjoint{<:Any,<:Ridiagonal})
	Matrix(adjoint(sparse(adjR.parent)))
end

function Matrix(invR::Inv{<:Any,<:Ridiagonal})
	inv(Matrix(invR.parent))
end

function Matrix(inv_adj_R::Inv{<:Any,<:Adjoint{<:Any,<:Ridiagonal}})
	inv(adjoint(Matrix(inv_adj_R.parent.parent)))
end

# ## Show
function Base.show(io::IO, m::MIME"text/plain", R::Ridiagonal)
    # show(io, m, sparse(R))
    println(io, "Ridiagonal:")
    X        = sparse(R)
    io       = IOContext(io, :typeinfo => eltype(X))
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, X)
end

function Base.show(io::IO, m::MIME"text/plain", adjR::Adjoint{<:Any,<:Ridiagonal})
    # show(io, m, sparse(R))
    println(io, "Adjoint Ridiagonal:")
    X        = sparse(adjR)
    io       = IOContext(io, :typeinfo => eltype(X))
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, X)
end

# sparse, matrix and show for Midiagonal
# =================================================

# ## Sparse
sparse(M::Midiagonal) = sparse(mortar(Diagonal(M.data)))

function sparse(adjM::Adjoint{<:Any,<:Midiagonal})
	adjoint(sparse(adjM.parent))
end

function sparse(invM::Inv{<:Any,<:Midiagonal})
	sparse(Midiagonal(map(inv,invM.parent)))
end

function sparse(inv_adj_R::Inv{<:Any,<:Adjoint{<:Any,<:Midiagonal}})
	sparse(Midiagonal(map(x->inv(x'),invM.parent)))
end

# ## Matrix
Matrix(M::Midiagonal) = Matrix(sparse(M))

Matrix(M::Adjoint{<:Any,<:Midiagonal}) = Matrix(sparse(M))

Matrix(M::Inv{<:Any,<:Adjoint{<:Any,<:Midiagonal}}) = Matrix(sparse(M))

# ## show

function Base.show(io::IO, m::MIME"text/plain", M::Midiagonal)
    println(io, "Midiagonal:")
    X        = sparse(M)
    io       = IOContext(io, :typeinfo => eltype(X))
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, X)
end

# sparse, matrix and show for Piv
# =================================================

function sparse(p::Piv)
    n = length(p.perm)
    sparse(1:n, p.perm, fill(true,n))
end

Matrix(p::Piv) = Matrix(sparse(p))

# Note: we don't need any more methods for adjoint or inv of Piv
# since they always just return another Piv for the return values.

function Base.show(io::IO, m::MIME"text/plain", p::Piv)
    println(io, "Piv permutation type:")
    X        = sparse(p)
    io       = IOContext(io, :typeinfo => eltype(X))
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, X)
end


text:
__precompile__()

module FredApi

export
    set_api_key,
    get_symbols,
    search_symbol,
    get_release,
    get_metadata,
    get_category,
    get_sources,
    get_source

using DataFrames
using TimeSeries
using Dates
using HTTP
using JSON

include("set_api_key.jl")
include("get_symbols.jl")
include("search_symbol.jl")
include("get_release.jl")
include("get_metadata.jl")
include("get_category.jl")
include("get_sources.jl")
include("get_source.jl")

end


text:
using FermiCG
using Documenter

pages = [
    "Home" => "index.md",
    "Installation Instructions" => "installation_instructions.md",
    # "Code Basics" => "basics.md",
    # "Grids" => "grids.md",
    # "Problem" => "problem.md",
    # "GPU" => "gpu.md",
    "Examples" => ["cmf.md","fci.md"],
#    "Library" => [
#        "Contents" => "library/outline.md",
#        "Public" => "library/public.md",
#        "Private" => "library/internals.md",
#        "Function index" => "library/function_index.md",
#        ],
"Functions" => Any[
                   #"FermiCG" => "library/FermiCG.md",
                   "CMF" => "library/CMFs.md",
                   "TPSCI" => "library/TPSCI.md",
                   "BST" => "library/BST.md",
                   "Internals" => "library/Internals.md",
                   "ClusteredTerms" => "library/ClusteredTerms.md",
                   "States" => "library/States.md",
                   "PyscfFunctions" => "library/PyscfFunctions.md",
                   "StringCI" => "library/StringCI.md",
                   "Utils" => "library/Utils.md",
                  ],
]

#####
##### Generate examples
#####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

examples = [
    "test_cmf.jl",
    "test_fci.jl",
]

# for example in examples
#   example_filepath = joinpath(EXAMPLES_DIR, example)
#   withenv("GITHUB_REPOSITORY" => "FourierFlows/FourierFlowsDocumentation") do
#     example_filepath = joinpath(EXAMPLES_DIR, example)
#     Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true)
#     Literate.notebook(example_filepath, OUTPUT_DIR, documenter=true)
#     Literate.script(example_filepath, OUTPUT_DIR, documenter=true)
#   end
# end

makedocs(
    modules=[FermiCG],
    authors="Nick Mayhall <nmayhall@vt.edu> and contributors",
    #repo="https://github.com/nmayhall-vt/FermiCG/blob/{commit}{path}#L{line}",
    sitename="FermiCG",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nmayhall-vt.github.io/FermiCG/stable",
        assets=String[],
    ),
    #html_prettyurls = !("local" in ARGS),
    pages=pages,
)

deploydocs(
    repo="github.com/nmayhall-vt/FermiCG.git",
    branch = "gh-pages",
    devbranch = "main",
    #push_preview = true,
    target= "build",
    deps = Deps.pip("pygments", "mkdocs", "python-markdown-math")
)

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
"""
```
impulse_responses(m, horizon, shock_ind, shock_size, z₀; deviations = true)
impulse_responses(m, horizon, shock_ind, shock_size; deviations = true)
```

calculates impulse responses according to a risk-adjusted linearization,
given a size for the shock, the index of the shock, and the initial state.
The second method assumes that the initial state is the stochastic steady state.

### Inputs
- `m::RiskAdjustedLinearization`: a solved risk-adjusted linearization of a dynamic economy
- `z₀::AbstractVector`: the initial state from which the economy begins
- `horizon::Int`: number of periods to be simulated
- `shock_ind::Int`: index of the shock that should be nonzero (other shocks are zero)
- `shock_size::Number`: size of the shock (can be positive or negative)

### Keywords
- `deviations::Bool`: if true, the impulse responses are returned in deviations from steady state.

### Outputs
- `states`: a matrix of the simulated path of states `z`, with type specified by the array type of `z₀`
- `jumps`: a matrix of the simulated path of jump variables `y`, with type specified by the array type of `z₀`
"""
function impulse_responses(m::RiskAdjustedLinearization, horizon::Int, shock_ind::Int,
                           shock_size::Number, z₀::AbstractVector;
                           deviations::Bool = true)

    # Create shock vector and output matrices
    shock            = zeros(eltype(z₀), m.Nε, 1)
    shock[shock_ind] = shock_size
    states           = similar(z₀, m.Nz, horizon)
    jumps            = similar(z₀, m.Ny, horizon)

    # Calculate state after impact
    states[:, 1], jumps[:, 1] = simulate(m, 1, shock, z₀)

    # Simulate with no other shocks drawn
    states[:, 2:end], jumps[:, 2:end] = simulate(m, horizon - 1, (@view states[:, 1]))

    if deviations
        return states .- m.z, jumps .- m.y
    else
        return states, jumps
    end
end

function impulse_responses(m::RiskAdjustedLinearization, horizon::Int, shock_ind::Int, shock_size::Number;
                           deviations::Bool = true)

    # Calculate starting at stochastic steady state
    return impulse_responses(m, horizon, shock_ind, shock_size, m.z; deviations = deviations)
end


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPMa27

import ..MadNLPHSL:
    @kwdef, Logger, @debug, @warn, @error, libhsl,
    AbstractOptions, AbstractLinearSolver, set_options!, SparseMatrixCSC, SubVector, StrideOneVector,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, findIJ, nnz

const ma27_default_icntl = Int32[
    6,6,0,2139062143,1,32639,32639,32639,32639,14,9,8,8,9,10,32639,32639,32639,32689,24,11,9,8,9,10,0,0,0,0,0]
const ma27_default_cntl  = [.1,1.0,0.,0.,0.]
const INPUT_MATRIX_TYPE = :csc

@kwdef mutable struct Options <: AbstractOptions
    ma27_pivtol::Float64 = 1e-8
    ma27_pivtolmax::Float64 = 1e-4
    ma27_liw_init_factor::Float64 = 5.
    ma27_la_init_factor::Float64 =5.
    ma27_meminc_factor::Float64 =2.
end

mutable struct Solver <: AbstractLinearSolver
    csc::SparseMatrixCSC{Float64,Int32}
    I::Vector{Int32}
    J::Vector{Int32}

    icntl::Vector{Int32}
    cntl::Vector{Float64}

    info::Vector{Int32}

    a::Vector{Float64}
    a_view::StrideOneVector{Float64}
    la::Int32
    ikeep::Vector{Int32}

    iw::Vector{Int32}
    liw::Int32
    iw1::Vector{Int32}
    nsteps::Vector{Int32}
    w::Vector{Float64}
    maxfrt::Vector{Int32}

    opt::Options
    logger::Logger
end

ma27ad!(n::Cint,nz::Cint,I::StrideOneVector{Cint},J::StrideOneVector{Cint},
        iw::Vector{Cint},liw::Cint,ikeep::Vector{Cint},iw1::Vector{Cint},
        nsteps::Vector{Cint},iflag::Cint,icntl::Vector{Cint},cntl::Vector{Cdouble},
        info::Vector{Cint},ops::Cdouble) = ccall(
            (:ma27ad_,libhsl),
            Nothing,
            (Ref{Cint},Ref{Cint},Ptr{Cint},Ptr{Cint},
             Ptr{Cint},Ref{Cint},Ptr{Cint},Ptr{Cint},
             Ptr{Cint},Ref{Cint},Ptr{Cint},Ptr{Cdouble},
             Ptr{Cint},Ref{Cdouble}),
            n,nz,I,J,iw,liw,ikeep,iw1,nsteps,iflag,icntl,cntl,info,ops)

ma27bd!(n::Cint,nz::Cint,I::StrideOneVector{Cint},J::StrideOneVector{Cint},
        a::StrideOneVector{Cdouble},la::Cint,iw::Vector{Cint},liw::Cint,
        ikeep::Vector{Cint},nsteps::Vector{Cint},maxfrt::Vector{Cint},iw1::Vector{Cint},
        icntl::Vector{Cint},cntl::Vector{Cdouble},info::Vector{Cint}) = ccall(
            (:ma27bd_,libhsl),
            Nothing,
            (Ref{Cint},Ref{Cint},Ptr{Cint},Ptr{Cint},
             Ptr{Cdouble},Ref{Cint},Ptr{Cint},Ref{Cint},
             Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},
             Ptr{Cint},Ptr{Cdouble},Ptr{Cint}),
            n,nz,I,J,a,la,iw,liw,ikeep,nsteps,maxfrt,iw1,icntl,cntl,info)

ma27cd!(n::Cint,a::Vector{Cdouble},la::Cint,iw::Vector{Cint},
        liw::Cint,w::Vector{Cdouble},maxfrt::Vector{Cint},rhs::Vector{Cdouble},
        iw1::Vector{Cint},nsteps::Vector{Cint},icntl::Vector{Cint},
        info::Vector{Cint}) = ccall(
            (:ma27cd_,libhsl),
            Nothing,
            (Ref{Cint},Ptr{Cdouble},Ref{Cint},Ptr{Cint},
             Ref{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble},
             Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
            n,a,la,iw,liw,w,maxfrt,rhs,iw1,nsteps,icntl,info)

function Solver(csc::SparseMatrixCSC;
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),kwargs...)

    set_options!(opt,option_dict,kwargs)

    I,J = findIJ(csc)
    nz=Int32(nnz(csc))

    liw= Int32(2*(2*nz+3*csc.n+1))
    iw = Vector{Int32}(undef,liw)
    ikeep= Vector{Int32}(undef,3*csc.n)
    iw1  = Vector{Int32}(undef,2*csc.n)
    nsteps=Int32[1]
    iflag =Int32(0)

    icntl= copy(ma27_default_icntl)
    icntl[1:2] .= 0
    cntl = copy(ma27_default_cntl)
    cntl[1] = opt.ma27_pivtol

    info = Vector{Int32}(undef,20)
    ma27ad!(Int32(csc.n),nz,I,J,iw,liw,ikeep,iw1,nsteps,Int32(0),icntl,cntl,info,0.)
    info[1]<0 && throw(SymbolicException())

    la = ceil(Int32,max(nz,opt.ma27_la_init_factor*info[5]))
    a = Vector{Float64}(undef,la)
    a_view = view(a,1:nnz(csc))
    liw= ceil(Int32,opt.ma27_liw_init_factor*info[6])
    resize!(iw,liw)
    maxfrt=Int32[1]

    return Solver(csc,I,J,icntl,cntl,info,a,a_view,la,ikeep,iw,liw,
                  iw1,nsteps,Vector{Float64}(),maxfrt,opt,logger)
end

function factorize!(M::Solver)
    M.a_view.=M.csc.nzval
    while true
        ma27bd!(Int32(M.csc.n),Int32(nnz(M.csc)),M.I,M.J,M.a,M.la,
                M.iw,M.liw,M.ikeep,M.nsteps,M.maxfrt,
                M.iw1,M.icntl,M.cntl,M.info)
        if M.info[1] == -3
            M.liw = ceil(Int32,M.opt.ma27_meminc_factor*M.liw)
            resize!(M.iw, M.liw)
            @debug(M.logger,"Reallocating memory: liw ($(M.liw))")
        elseif M.info[1] == -4
            M.la = ceil(Int32,M.opt.ma27_meminc_factor*M.la)
            resize!(M.a,M.la)
            @debug(M.logger,"Reallocating memory: la ($(M.la))")
        elseif M.info[1] < 0
            throw(FactorizationException())
        else
            break
        end
    end
    return M
end

function solve!(M::Solver,rhs::StrideOneVector{Float64})
    length(M.w)<M.maxfrt[1] && resize!(M.w,M.maxfrt[1])
    length(M.iw1)<M.nsteps[1] && resize!(M.iw1,M.nsteps[1])
    ma27cd!(Int32(M.csc.n),M.a,M.la,M.iw,M.liw,M.w,M.maxfrt,rhs,
            M.iw1,M.nsteps,M.icntl,M.info)
    M.info[1]<0 && throw(SolveException())
    return rhs
end

is_inertia(::Solver) = true
function inertia(M::Solver)
    rank = M.info[1]==3 ? M.info[2] : rank = M.csc.n
    return (rank-M.info[15],M.csc.n-rank,M.info[15])
end

function improve!(M::Solver)
    if M.cntl[1] == M.opt.ma27_pivtolmax
        @debug(M.logger,"improve quality failed.")
        return false
    end
    M.cntl[1] = min(M.opt.ma27_pivtolmax,M.cntl[1]^.75)
    @debug(M.logger,"improved quality: pivtol = $(M.cntl[1])")
    return true
end

introduce(::Solver)="ma27"

end # module


text:
module LaplaceBeltramiCubedSphereTests
  # This program solves the Laplace-Beltrami PDE on a unit sphere.
  # In strong form, the problem can be stated as:
  # find u s.t.
  #   -Δ_S(u) = f on \Omega
  # Where Δ_S is the spherical Laplacian defined as Δ_S = ∇_S⋅(∇_S(u)),
  # with ∇_S() and  ∇_S⋅ being the spherical gradient and spherical
  # divergence, resp.

  # Given a known analytical solution u(x), the program manufactures f(x)
  # by applying -Δ_S to u, and then performs an error convergence analysis.
  # For the discretization of the sphere, it uses the so-called cubed sphere
  # mesh
  using Test
  using Gridap
  import Gridap.Fields: ∇
  using GridapGeosciences
  using Plots

  include("ConvergenceAnalysisTools.jl")

  function u(x)
    sin(π*x[1])*cos(π*x[2])*exp(x[3])
  end
  function uθϕ(θϕ)
    u(θϕ2xyz(θϕ))
  end
  function f(x)
    -(laplacian_unit_sphere(uθϕ)(xyz2θϕ(x)))
  end
  function ∇(::typeof(u))
    gradient_unit_sphere(uθϕ)
  end

  function solve_laplace_beltrami(model,order,degree)
     V = FESpace(model,ReferenceFE(lagrangian,Float64,order); conformity=:H1)
     U = TrialFESpace(V)

     Ω = Triangulation(model)
     dΩ = Measure(Ω,degree)

     a(u,v) = ∫(∇(v)⋅∇(u))dΩ
     b(v)   = ∫(v*f)dΩ

     op = AffineFEOperator(a,b,U,V)
     uh = solve(op)

     e=u-uh
     ∇e = ∇(uh)-∇(u)∘xyz2θϕ

     # H1-norm
     # sqrt(sum(∫( e*e + (∇e)⋅(∇e) )dΩ))
     sqrt(sum(∫( (∇e)⋅(∇e) )dΩ))#,uh
  end

  #model=CubedSphereDiscreteModel(50)
  #e,uh=solve_laplace_beltrami(model,1,8)
  #writevtk(Triangulation(model),"u",nsubcells=4,cellfields=["u"=>u,"uh"=>uh])

  @time ahs1,ak1errors,as1=convergence_study(solve_laplace_beltrami,generate_n_values(2),1,8)
  @test round(as1,digits=1) ≈ 1.0

  @time ahs2,ak2errors,as2=convergence_study(solve_laplace_beltrami,generate_n_values(2),2,12)
  @test round(as2,digits=1) ≈ 2.0

  #  plot([ahs1,ahs2],[ak1errors,ak2errors],
  #  xaxis=:log, yaxis=:log,
  #  label=["k=1 analytical map" "k=2 analytical map"],
  #  shape=:auto,
  #  xlabel="h",ylabel="error H1 semi-norm", legend=:bottomright)

end


text:
using ArgParse
using CSV
using DataFrames
using DelimitedFiles
using LinearAlgebra
using Mamba
using MicrobiomeMixedModels
using Random
using YAML

function parse_commandline()
    s = ArgParseSettings("Fit a model to real-world data.")
    @add_arg_table! s begin
        "input"
            help = "Directory from which to read X.csv, Y.csv, and Z.csv."
        "output"
            help = "Directory to which results are written."
        "--iters", "-i"
            arg_type = Int
            default = 500
            help = "Number of MCMC iterations."
        "--burnin", "-b"
            arg_type = Int
            default = 0
            help = "Number of initial MCMC iterations to remove as burn-in."
        "--thin", "-t"
            arg_type = Int
            default = 1
            help = "Retain one of every 'thin' iterations."
        "--chains", "-c"
            arg_type = Int
            default = 1
            help = "Number of MCMC chains to run."
        "--seed"
            arg_type = Int
            default = 1
            help = "Reseed the random number generator."
        "--monitor"
            nargs = '*'
            help = "YAML defining the nodes to monitor through MCMC."
        "--inits"
            nargs = '*'
            help = "YAML defining initial values of stochastic nodes."
        "--hyper"
            nargs = '*'
            help = "YAML defining hyperparameters of stochastic node priors."
        "--factors"
            arg_type = Int
            default = 0
            help = "Number of factors to use in fitting the MIMIX model."
        "--loadings"
            default = "DL"
            help = "Distributions on loadings (G for Gaussian priors, DL for Dirichlet-Laplace priors)."
        "--permanova"
            help = "Run PERMANOVA with vegan::adonis() in R."
            action = :store_true
        "--post-pred-check"
            help = "Perform posterior predictive checks, requires at least one of :Y, :ϕ, or :θ are monitored."
            action = :store_true
    end
    return parse_args(s)
end

function read_data(dir; L=0)
    println("Reading X.csv")
    X = convert(Matrix{Int}, DataFrame!(CSV.File(joinpath(dir, "X.csv"), datarow=1)))
    println("Reading Y.csv")
    Y = convert(Matrix{Int}, DataFrame!(CSV.File(joinpath(dir, "Y.csv"), datarow=1)))
    println("Reading Z.csv")
    Z = convert(Matrix{Int}, DataFrame!(CSV.File(joinpath(dir, "Z.csv"), datarow=1)))
    N, K = size(Y)
    m = sum(Y, dims=2)
    q = maximum(vec(Z))
    data = Dict{Symbol, Any}(
        :X => X,
        :Y => Y,
        :Z => Z,
        :N => N,
        :K => K,
        :m => m,
        :q => q,
        :L => L
    )
    data[:p] = size(X, 2)
    data[:num_blocking_factors] = size(Z, 2)
    data[:blocking_factor] = Dict{Int, Int}()
    for level in unique(Z)
        bf = findall(vec(any(Z .== level, dims=1)))
        @assert length(bf) == 1
        data[:blocking_factor][level] = bf[1]
    end

    return data
end

args = parse_commandline()

input = abspath(args["input"])
@assert isdir(input)

output = abspath(args["output"])
Random.seed!(args["seed"])

if args["permanova"]
    data = read_data(input)

    println("Beginning PERMANOVA test with vegan::adonis in R.")
    sim = permanova(data)
    global_results = DataFrame(pval = sim[:pval])

    global_test_path = joinpath(output, "global-test.tsv")
    println("Writing global test results to $global_test_path")
    CSV.write(global_test_path, global_results, delim='\t')
else  # mimix
    @assert 0 < args["iters"]   "Iters must be positive"
    @assert 0 <= args["burnin"] "Burn-in must be non-negative"
    @assert 0 < args["thin"]    "Thin must be positive"
    @assert 0 < args["chains"]  "Chains must be positive"

    factors = args["factors"]
    loadings = args["loadings"]
    if (factors > 0) & (loadings == "DL")
        model_type = MIMIX(factors)
    elseif (factors > 0) & (loadings == "G")
        model_type = MIMIXGaussian(factors)
    elseif factors == 0
        model_type = MIMIXNoFactors()
    else
        ValueError("--factors requires a non-negative integer")
    end

    monitor_conf = load_config(AbstractString[abspath(data_path) for data_path in args["monitor"]])
    hyper_conf = load_config(AbstractString[abspath(data_path) for data_path in args["hyper"]])
    inits_conf = load_config(AbstractString[abspath(data_path) for data_path in args["inits"]])

    if args["post-pred-check"] && !monitor_conf[:θ] && (!monitor_conf[:θ_mean] || !monitor_conf[:θ_var])
        error("Posterior prediction checks require monitoring θ or both θ_mean and θ_var, please edit the monitor config.")
    end

    model = get_model(model_type, monitor_conf, hyper_conf)

    data = read_data(input; L = factors)

    inits = get_inits(model_type, inits_conf, data)
    inits = [inits for _ in 1:args["chains"]]

    println("Beginning MCMC sampling")
    mcmc_kwargs = Dict(Symbol(key) => args[key] for key in ["burnin", "thin", "chains"])
    sim = mcmc(model, data, inits, args["iters"]; mcmc_kwargs...)

    monitor_Y = pop!(monitor_conf, :Y)
    monitor_ϕ = pop!(monitor_conf, :ϕ)
    monitor_θ = pop!(monitor_conf, :θ)
    monitor_θ_mean = pop!(monitor_conf, :θ_mean)
    monitor_θ_var = monitor_conf[:θ_var]
    monitor_Λ = pop!(monitor_conf, :Λ)

    for (param, monitor) in monitor_conf
        if monitor
            post = get_post(sim, data, param)
            param = String(param)
            for (letter, name) in MicrobiomeMixedModels.latin
                if startswith(param, letter)
                    param = replace(param, letter => name)
                end
            end
            post_path = joinpath(output, "$param.tsv")
            println("Saving posterior samples of $param to $post_path")
            CSV.write(post_path, post, delim='\t')
        end
    end

    if monitor_Λ
        # Lambda is large, so rather than save every sample we save the posterior mean
        Λ = get_post(sim, data, :Λ)
        Λ_values = convert(Matrix, Λ)
        Λ_postmean = reshape(vec(mean(Λ_values, dims=1)), data[:L], data[:K])
        post_path = joinpath(output, "Lambda-postmean.tsv")
        println("Saving posterior mean of Lambda to $post_path")
        writedlm(post_path, Λ_postmean)

        # We can also store the posterior mean of Lambda' * Lambda
        # ΛtΛs = Matrix{Float64}[]
        # for i in 1:size(Λ_values, 1)
        #     Λ_mat = reshape(Λ_values[i, :], data[:L], data[:K])
        #     push!(ΛtΛs, transpose(Λ_mat) * Λ_mat)
        # end
        # ΛtΛ_postmean = sum(ΛtΛs, dims=3)
        # ΛtΛ_postmean ./= size(Λ_values, 1)
        # post_path = joinpath(output, "LambdatLambda-postmean.tsv")
        # println("Saving posterior mean of Lambda' * Lambda to $post_path")
        # writedlm(post_path, ΛtΛ_postmean)
    end

    if args["post-pred-check"]
        # Posterior predictive checks on sparsity, overdispersion, and ecological diversity (alpha and beta)

        function shannon_diversity(x::Vector{T}) where T <: Real
            p = x ./ sum(x)
            p = p[p .> 0.]
            return -sum(p .* log.(p))
        end

        function simpson_diversity(x::Vector{T}) where T <: Real
            p = x ./ sum(x)
            return sum(abs2, p)
        end

        function braycurtis_diversity(x::Vector{T}, y::Vector{T}) where T <: Real
            c = zero(eltype(x))
            s = zero(eltype(x))
            for i in eachindex(x, y)
                c += min(x[i], y[i])
                s += x[i] + y[i]
            end
            return 2c / s
        end

        function jaccard_diversity(x::Vector{T}, y::Vector{T}) where T <: Real
            p = 0
            q = 0
            r = 0
            z = zero(T)
            for i in eachindex(x, y)
                if (x[i] > z) && (y[i] > z)
                    p += 1
                elseif (x[i] > z) && (y[i] == z)
                    q += 1
                elseif (x[i] == z) && (y[i] > z)
                    r += 1
                end
            end
            return p / (p + q + r)
        end

        function calculate_mean_beta_div_by_group(Y, groups, diversity_fn)
            @assert size(Y, 1) == length(groups)
            mean_div = zeros(maximum(group_ids))
            for group_id in unique(group_ids)
                Y_group = Y[group_id .== group_ids, :]
                pairwise_div = Float64[]
                for i in 1:size(Y_group, 1)
                    for j in i+1:size(Y_group, 1)
                        push!(pairwise_div, diversity_fn(Y_group[i, :], Y_group[j, :]))
                    end
                end
                mean_div[group_id] = mean(pairwise_div)
            end
            return mean_div
        end

        matchrows(A::Matrix, B::Matrix) = mapslices(findall, hcat([all(A .== transpose(B[i, :]), dims=2) for i in 1:size(B, 1)]...), dims=2)

        full_combos = hcat(data[:X], data[:Z])
        group_ids = vec(matchrows(full_combos, unique(full_combos, dims=1)))  # each group is unique combo of treatment, site, and blocks

        obs_max_count = vec(maximum(data[:Y], dims=2) ./ sum(data[:Y], dims=2))
        obs_mean_count = vec(var(data[:Y], dims=2))
        obs_prop_eq_zero = vec(mean(data[:Y] .== 0.0, dims=2))
        obs_prop_leq_one = vec(mean(data[:Y] .<= 1.0, dims=2))
        obs_prop_leq_two = vec(mean(data[:Y] .<= 2.0, dims=2))
        obs_shannon_div = mapslices(shannon_diversity, data[:Y], dims=2)
        obs_simpson_div = mapslices(simpson_diversity, data[:Y], dims=2)
        obs_braycurtis_div = reshape(calculate_mean_beta_div_by_group(data[:Y], group_ids, braycurtis_diversity), maximum(group_ids), 1)
        obs_jaccard_div = reshape(calculate_mean_beta_div_by_group(data[:Y], group_ids, jaccard_diversity), maximum(group_ids), 1)
        writedlm(joinpath(output, "obs-max-count.tsv"), obs_max_count)
        writedlm(joinpath(output, "obs-mean-count.tsv"), obs_mean_count)
        writedlm(joinpath(output, "obs-prop-eq-zero.tsv"), obs_prop_eq_zero)
        writedlm(joinpath(output, "obs-prop-leq-one.tsv"), obs_prop_leq_one)
        writedlm(joinpath(output, "obs-prop-leq-two.tsv"), obs_prop_leq_two)
        writedlm(joinpath(output, "obs-shannon-div.tsv"), obs_shannon_div)
        writedlm(joinpath(output, "obs-simpson-div.tsv"), obs_simpson_div)
        writedlm(joinpath(output, "obs-braycurtis-div.tsv"), obs_braycurtis_div)
        writedlm(joinpath(output, "obs-jaccard-div.tsv"), obs_jaccard_div)

        post_pred_checks = Dict{String, Matrix}()

        if monitor_θ
            θ_post = convert(Matrix, get_post(sim, data, :θ))
            post_pred_checks["theta_mcmc"] = θ_post
        end

        if monitor_θ_mean && monitor_θ_var
            θ_mean_post = convert(Matrix, get_post(sim, data, :θ_mean))
            θ_var_post = convert(Matrix, get_post(sim, data, :θ_var))
            iters = size(θ_mean_post, 1)
            θ_post = zeros(iters, data[:N] * data[:K])
            for iter in 1:iters
                θ_mean = reshape(θ_mean_post[iter, :], data[:N], data[:K])
                θ_new = [rand(MvNormal(θ_mean[i, :], sqrt.(θ_var_post[iter, :]))) for i in 1:data[:N]]
                θ_post[iter, :] = vec(transpose(hcat(θ_new...)))
            end
            post_pred_checks["theta_prior"] = θ_post
        end

        for post_pred_check in post_pred_checks
            name, θ_post = post_pred_check

            n_iter = size(θ_post, 1)
            max_count = zeros(n_iter, data[:N])
            mean_count = zeros(n_iter, data[:N])
            prop_eq_zero = zeros(n_iter, data[:N])
            prop_leq_one = zeros(n_iter, data[:N])
            prop_leq_two = zeros(n_iter, data[:N])
            shannon_div = zeros(n_iter, data[:N])
            simpson_div = zeros(n_iter, data[:N])
            braycurtis_div = zeros(n_iter, maximum(group_ids))
            jaccard_div = zeros(n_iter, maximum(group_ids))
            for iter in 1:n_iter
                ϕ = clr(reshape(θ_post[iter, :], data[:N], data[:K]))
                Y_pred_iter = transpose(hcat([rand(Multinomial(data[:m][i], ϕ[i, :])) for i in 1:data[:N]]...))
                max_count[iter, :] = vec(maximum(Y_pred_iter, dims=2) ./ sum(Y_pred_iter, dims=2))
                mean_count[iter, :] = vec(var(Y_pred_iter, dims=2))
                prop_eq_zero[iter, :] = vec(mean(Y_pred_iter .== 0.0, dims=2))
                prop_leq_one[iter, :] = vec(mean(Y_pred_iter .<= 1.0, dims=2))
                prop_leq_two[iter, :] = vec(mean(Y_pred_iter .<= 2.0, dims=2))
                shannon_div[iter, :] = mapslices(shannon_diversity, Y_pred_iter, dims=2)
                simpson_div[iter, :] = mapslices(simpson_diversity, Y_pred_iter, dims=2)
                braycurtis_div[iter, :] = calculate_mean_beta_div_by_group(Y_pred_iter, group_ids, braycurtis_diversity)
                jaccard_div[iter, :] = calculate_mean_beta_div_by_group(Y_pred_iter, group_ids, jaccard_diversity)
            end

            if !isdir(joinpath(output, name))
                mkdir(joinpath(output, name))
            end

            writedlm(joinpath(output, name, "post-pred-max-count.tsv"), max_count)
            writedlm(joinpath(output, name, "post-pred-mean-count.tsv"), mean_count)
            writedlm(joinpath(output, name, "post-pred-prop-eq-zero.tsv"), prop_eq_zero)
            writedlm(joinpath(output, name, "post-pred-prop-leq-one.tsv"), prop_leq_one)
            writedlm(joinpath(output, name, "post-pred-prop-leq-two.tsv"), prop_leq_two)
            writedlm(joinpath(output, name, "post-pred-shannon-div.tsv"), shannon_div)
            writedlm(joinpath(output, name, "post-pred-simpson-div.tsv"), simpson_div)
            writedlm(joinpath(output, name, "post-pred-braycurtis-div.tsv"), braycurtis_div)
            writedlm(joinpath(output, name, "post-pred-jaccard-div.tsv"), jaccard_div)
        end
    end
end

max_stars_repo_name: nsgrantham/mimix

max_issues_repo_name: nsgrantham/mimix

max_forks_repo_name: nsgrantham/mimix


text:
# Interpolation Lagrange

#Lagrange Polynomial 1D
function lagrangePoly(order::Int, i::Int, ξ::Float64)
    n = order+1

    ξ_vec = collect(range(-1, stop=1, length=n))
    if order != 1
        ξ_vec = [ξ_vec[1];ξ_vec[end];ξ_vec[2:end-1]]
    end

    Li = 1
    ξi = ξ_vec[i]
    for j=1:n
        ξj = ξ_vec[j]
        if i!=j
            Li *= (ξ - ξj)/(ξi - ξj)
        end
    end
    Li
end

##################################
# Lagrange dim 1 RefLine order>1 #
##################################
using Ferrite
import Ferrite: getnbasefunctions, nvertexdofs, faces, reference_coordinates, value

# shape
struct RefLine <: Ferrite.AbstractRefShape end

# Lagrange {dim,shape,order} -> dim=1, shape=RefLine , order >= 1
getnbasefunctions(::Lagrange{1,RefLine,order}) where {order} = order+1  #
nvertexdofs(::Lagrange{1,RefLine,order}) where {order} = 1  #
faces(::Lagrange{1,RefLine,order}) where {order} = ((1,), (2,)) #tuple([(i,) for i in 1:order+1]...)

function reference_coordinates(::Lagrange{1,RefLine,order}) where {order}
    order > 20 && throw(ArgumentError("The order must be less than 20"))
    n=order+1
    ξ_vec = collect(range(-1, stop=1, length=n))
    if order != 1
        ξ_vec = [ξ_vec[1];ξ_vec[end];ξ_vec[2:end-1]]
    end
    return [Vec{1, Float64}((ξ_vec[i],)) for i in 1:n]
end

function value(ip::Lagrange{1,RefLine,order}, i::Int, ξ::Vec{dim}) where {dim,order}
    order > 20 && throw(ArgumentError("The order must be less than 20"))
    ξ_x = ξ[1]
    return lagrangePoly(order, i, ξ_x)
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

"""
    initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D)

A Gaussian pulse in the density with constant velocity and pressure; reduces the
compressible Euler equations to the linear advection equations.
"""
function initial_condition_density_pulse(x, t, equations::CompressibleEulerEquations3D)
  rho = 1 + exp(-(x[1]^2 + x[2]^2 + x[3]^2))/2
  v1 = 1
  v2 = 1
  v3 = 1
  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_v3 = rho * v3
  p = 1
  rho_e = p/(equations.gamma - 1) + 1/2 * rho * (v1^2 + v2^2 + v3^2)
  return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e)
end
initial_condition = initial_condition_density_pulse

volume_flux = flux_ranocha
solver = DGSEM(polydeg=3, surface_flux=flux_ranocha,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-2.0, -2.0, -2.0)
coordinates_max = ( 2.0,  2.0,  2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                n_cells_max=100_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

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

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.1)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
using JuMP, Cbc, CSV

file = "campaignData.csv"

data_import = CSV.read(file, header=false)
campData = Matrix{Int}(data_import)

totalsColumn = size(campData,2)
numOptions = size(campData,2)-1

customersData = campData[1,1:numOptions]
resourcesData = campData[2:end,1:numOptions]
totalsData = campData[2:end,totalsColumn]

adModel = Model(Cbc.Optimizer)

@variable(adModel, options[1:numOptions], binary=true)
@objective(adModel, Max, sum(options[:].*customersData[:]))

for r = 1:size(resourcesData,1)
    @constraint(adModel, sum(resourcesData[r,:].*options[:])<=totalsData[r])
end
@constraints(adModel, begin
    options[6]-(options[4]+options[5])<=0
    options[2]+options[4]<=1
end)

stats = JuMP.optimize!(adModel)
result = JuMP.value.(options)

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
"""
    ==(a::Node, b::Node)

Test Node equality. The parent, children and siblings are not tested, only their id is.
"""
function Base.:(==)(a::T, b::T) where {T<:Node}
    isequal(a.name, b.name) &&
        isequal(a.id, b.id) &&
        isequal(
            isroot(a) ? nothing : parent(a).id,
            isroot(b) ? nothing : parent(b).id
        ) &&
        isequal(
            a.children !== nothing ? keys(a.children) : nothing,
            a.children !== nothing ? keys(a.children) : nothing
        ) &&
        isequal(
            a.siblings !== nothing ? keys(a.siblings) : nothing,
            a.siblings !== nothing ? keys(a.siblings) : nothing
        ) &&
        isequal(a.MTG, b.MTG) &&
        isequal(a.attributes, b.attributes)
end

"""
    ==(a::Node, b::Node)

Test AbstractNodeMTG equality.
"""
function Base.:(==)(a::T, b::T) where {T<:AbstractNodeMTG}
    isequal(a.link, b.link) &&
        isequal(a.symbol, b.symbol) &&
        isequal(a.index, b.index) &&
        isequal(a.scale, b.scale)
end


text:
# ---
# title: 939. Minimum Area Rectangle
# id: problem939
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/minimum-area-rectangle/description/>
# hidden: true
# ---
#
# Given a set of points in the xy-plane, determine the minimum area of a
# rectangle formed from these points, with sides parallel to the x and y axes.
#
# If there isn't any rectangle, return 0.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,1],[1,3],[3,1],[3,3],[2,2]]
#     Output: 4
#
#
# **Example 2:**
#
#
#
#     Input: [[1,1],[1,3],[3,1],[3,3],[4,1],[4,3]]
#     Output: 2
#
#
#
#
# **Note** :
#
#   1. `1 <= points.length <= 500`
#   2. `0 <= points[i][0] <= 40000`
#   3. `0 <= points[i][1] <= 40000`
#   4. All points are distinct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
    # Functions and methods for cheme::Cascade.DielectronicCaptureScheme computations

    """
    `Cascade.computeSteps(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested capture & transition amplitudes as well as DielectronicCapture.Line's, AutoIonization.Line's,
            etc. for all pre-specified decay steps of the cascade. When compared with standard computations of these atomic
            processes, however, the amount of output is largely reduced and often just printed into the summary file.
            A set of  data::Cascade.CaptureData  is returned.
    """
    function computeSteps(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        linesA = AutoIonization.Line[];    linesR = PhotoEmission.Line[];    cOrbitals = Dict{Subshell, Orbital}()
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        nt = 0;   st = 0;   previousMeanEn = 0.
        for  step  in  stepList
            st = st + 1
            nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
            sa = "\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc decay lines (without selection rules): "
            println(sa);    if  printSummary   println(iostream, sa)   end

            if      step.process == Basics.Auger()  &&   comp.approach == Cascade.AverageSCA()
                meanEn = 0.;    NoEn = 0
                for  p = 1:length(step.initialMultiplet.levels),  q = 1:length(step.finalMultiplet.levels)
                    en = step.initialMultiplet.levels[p].energy - step.finalMultiplet.levels[q].energy
                    if  en > 0.03    meanEn = meanEn + en;    NoEn = NoEn + 1    end
                end
                if  NoEn > 0     meanEn = meanEn/ NoEn
                else
                    println(">> No transition with positive energy.")
                    continue
                end

                # Generate potential for continuum orbitals for this step
                nrContinuum  = Continuum.gridConsistency(meanEn + 0.1, comp.grid)
                contSettings = Continuum.Settings(false, nrContinuum);
                npot         = Nuclear.nuclearPotential(comp.nuclearModel, comp.grid)
                ## wp1 = compute("radial potential: core-Hartree", grid, wLevel)
                ## wp2 = compute("radial potential: Hartree-Slater", grid, wLevel)
                wp           = Basics.computePotentialDFS(comp.grid, step.finalMultiplet.levels[1])
                pot          = Basics.add(npot, wp)
                #  Generate continuum if not yet available
                generatedKappas = Int64[]
                for  p = 1:length(step.initialMultiplet.levels),  q = 1:length(step.finalMultiplet.levels)
                    en = step.initialMultiplet.levels[p].energy - step.finalMultiplet.levels[q].energy
                    if  en > 0.03
                        symi      = LevelSymmetry(step.initialMultiplet.levels[p].J, step.initialMultiplet.levels[p].parity)
                        symf      = LevelSymmetry(step.finalMultiplet.levels[q].J, step.finalMultiplet.levels[q].parity)
                        kappaList = AngularMomentum.allowedKappaSymmetries(symi, symf)
                        for  kappa in kappaList
                            sh    = Subshell(101, kappa)
                            if      abs(kappa) > step.settings.maxKappa                               ## continuum orbital not needed
                            elseif  kappa in generatedKappas                                          ## already generated for this step
                            elseif  haskey(cOrbitals,sh)   &&   cOrbitals[sh].energy - en / en < 0.15 ## use previous one
                                println(">> No new continum orbital generated for $sh and energy $en ")
                            else                                                                      ## generate new continuum orbital
                                push!(generatedKappas, kappa)
                                cOrbital, phase, normF  = Continuum.generateOrbitalLocalPotential(en, sh, pot, contSettings)
                                cOrbitals[sh]           = cOrbital
                                println(">> New continum orbital generated for $sh and energy $en ")
                            end
                        end
                    end
                end

                newLines = AutoIonization.computeLinesFromOrbitals(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                                   step.settings, cOrbitals, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Auger()
                # Compute continuum orbitals independently for all transitions in the given block.
                newLines = AutoIonization.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                              step.settings, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Radiative()
                newLines = PhotoEmission.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.grid,
                                                             step.settings, output=true, printout=false)
                append!(linesR, newLines);    nt = length(linesR)
            else   error("Unsupported atomic process for cascade computations.")
            end
            sa = "     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                 "to a total of $nt $(string(step.process)) decay lines."
            println(sa);    if  printSummary   println(iostream, sa)   end
        end
        #
        data = Cascade.DecayData(linesR, linesA)
    end

    """
    `Cascade.determineSteps(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, captureList::Array{Cascade.Block,1}, decayList::Array{Cascade.Block,1})`
        ... determines all step::Cascade.Step's that need to be computed for this dielectronic cascade. It considers the autoionization
            between the blocks from the captureList and initialList as well as the radiative stabilization between the blocks from the
            captureList and decayList. It checks that at least on pair of levels supports either an `electron-capture' or
            `radiative stabilization' within the step. A stepList::Array{Cascade.Step,1} is returned, and for which subsequently all
            required transition amplitudes and rates/cross sections are computed.
    """
    function determineSteps(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, capturedList::Array{Cascade.Block,1}, decayList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            for  capturedBlock in capturedList
                for  initialBlock in initialList
                    if  initialBlock.NoElectrons + 1 != capturedBlock.NoElectrons   error("stop a")     end
                    # Check that at least one energy supports autoionization
                    settings = AutoIonization.Settings(AutoIonization.Settings(), maxKappa=7)
                    push!( stepList, Cascade.Step(Basics.Auger(), settings, capturedBlock.confs,     initialBlock.confs,
                                                                            capturedBlock.multiplet, initialBlock.multiplet) )
                end
            end
            #
            for  capturedBlock in capturedList
                for  decayBlock in decayList
                    if  decayBlock.NoElectrons != capturedBlock.NoElectrons   error("stop b")     end
                    # Check that at least one energy supports radiative stabilization
                    if   Basics.determineMeanEnergy(capturedBlock.multiplet) - Basics.determineMeanEnergy(decayBlock.multiplet) > scheme.minPhotonEnergy
                        settings = PhotoEmission.Settings(PhotoEmission.Settings(), gauges=[UseCoulomb, UseBabushkin])
                        push!( stepList, Cascade.Step(Basics.Radiative(), settings, capturedBlock.confs,     decayBlock.confs,
                                                                                capturedBlock.multiplet, decayBlock.multiplet) )
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.generateBlocks(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)`
        ... generate all block::Cascade.Block's, that need to be computed for this electron-capture and subsequent stabilization (DR) cascade,
            and compute also the corresponding multiplets. The different cascade approches realized different strategies how these blocks are
            selected and computed. A blockList::Array{Cascade.Block,1} is returned.
    """
    function generateBlocks(scheme::Cascade.DielectronicCaptureScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)
        blockList = Cascade.Block[];    basis = Basis()
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        if    comp.approach == AverageSCA()
            sa = "\n* Generate blocks for DR plasma rate coefficient computations: \n" *
                 "\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: " *
                 "\n    + orbitals are generated independently for the first block in a Dirac-Fock-Slater potential; " *
                 "\n    + these orbitals are re-used for all other block, together with hydrogenic orbitals for the outer part; " *
                 "\n    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; " *
                 "\n    + only the Coulomb interaction is considered for the electron capture. " *
                 "\n    + only E1 excitations are considered for the stabilization. \n"
            if  printout       println(sa)              end
            if  printSummary   println(iostream, sa)    end
            #
            if  length(confs) > 1
                # Determine a list of hydrogenic orbitals for later use
                relconfList = ConfigurationR[]
                for  confa in confs
                    wa = Basics.generate("configuration list: relativistic", confa)
                    append!( relconfList, wa)
                end
                subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
                ##x @show "aa1", subshellList
                Defaults.setDefaults("relativistic subshell list", subshellList; printout=printout)
                wa                 = Bsplines.generatePrimitives(comp.grid)
                hydrogenicOrbitals = Bsplines.generateOrbitalsHydrogenic(wa, comp.nuclearModel, subshellList; printout=printout)
            end

            for  (ia, confa)  in  enumerate(confs)
                sa = "  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... "
                print(sa);      if  printSummary   println(iostream, sa)   end
                # Now distinguish between the first and all other blocks; for the first block, a SCF is generated and the occupied orbital
                # used also for all other blocks. In addition, a set of hydrogenic orbitals generated for later use
                if  ia == 1
                    basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                else
                    # Generate a list of relativistic configurations and determine an ordered list of subshells for these configurations
                    relconfList  = ConfigurationR[]
                    wa           = Basics.generate("configuration list: relativistic", confa);    append!( relconfList, wa)
                    subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
                    Defaults.setDefaults("relativistic subshell list", subshellList; printout=false)
                    # Generate the relativistic CSF's for the given subshell list
                    csfList = CsfR[]
                    for  relconf in relconfList
                        newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
                        append!( csfList, newCsfs)
                    end
                    # Determine the list of coreSubshells
                    coreSubshellList = Subshell[]
                    for  k in 1:length(subshellList)
                        mocc = Basics.subshell_2j(subshellList[k]) + 1;    is_filled = true
                        for  csf in csfList
                            if  csf.occupation[k] != mocc    is_filled = false;    break   end
                        end
                        if   is_filled    push!( coreSubshellList, subshellList[k])    end
                    end
                    # Add all missing orbitals as hydrogenic
                    orbitals      = copy(basis.orbitals)
                    for  subsh in subshellList
                        if haskey(orbitals, subsh)   ## do nothing
                        else      orbitals[subsh] = hydrogenicOrbitals[subsh];   print("hydrogenic $subsh ...")
                        end
                    end

                    basis         = Basis(true, confa.NoElectrons, subshellList, csfList, coreSubshellList, orbitals)
                end
                multiplet = Basics.perform("computation: mutiplet from orbitals, no CI, CSF diagonal", [confa],  basis.orbitals,
                                            comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        elseif    comp.approach == SCA()
            sa = "\n* Generate blocks for electron-capture & decay cascade computations: \n" *
                 "\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: " *
                 "\n    + each single configuration forms an individual cascade block; " *
                 "\n    + orbitals are generated independently for each block for a Dirac-Fock-Slater potential; " *
                 "\n    + configuration mixing is included for each block, based on H^(DC); " *
                 "\n    + all requested multipoles are considered for the stabilization. \n"
            if  printout       println(sa)              end
            if  printSummary   println(iostream, sa)    end
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")   end
                    basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                    multiplet = Basics.performCI(basis,    comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        else  error("Unsupported cascade approach.")
        end

        return( blockList )
    end

    """
    `Cascade.generateConfigurationsForDielectronicCapture(multiplets::Array{Multiplet,1},  scheme::DielectronicCaptureScheme,
                                                      nm::Nuclear.Model, grid::Radial.Grid)`
        ... generates all possible doubly-excited configurations due to (dielectronic) electron capture into the given multiplets.
            The number and type of such doubly-generated configurations depend on (1) the maximum (electron) energy for capturing an electron
            that is closely related to the (maximum) temperature of the plasma; (2) the fromShells from which (and how many displacements)
            are accepted as well as (3) the maximum principle and orbital angular quantum number of the additional (to-) shellsthe fromShells
            into which electrons excited and/or captured. A Tuple(initialConfList::Array{Configuration,1}, confList::Array{Configuration,1})
            is returned.
    """
    function generateConfigurationsForDielectronicCapture(multiplets::Array{Multiplet,1},  scheme::DielectronicCaptureScheme,
                                                      nm::Nuclear.Model, grid::Radial.Grid)
        # Determine all (reference) configurations from multiplets and generate the 'excited' configurations due to the specificed excitations
        initialConfList = Configuration[]
        for mp  in  multiplets
            confList = Basics.extractNonrelativisticConfigurations(mp.levels[1].basis)
            for  conf in confList   if  conf in initialConfList   nothing   else   push!(initialConfList, conf)      end      end
        end
        coreConfList    = Basics.generateConfigurations(initialConfList, scheme.excitationFromShells, scheme.excitationToShells,
                                                        scheme.NoExcitations)
        captureConfList = Basics.generateConfigurationsWithElectronCapture(coreConfList, scheme.excitationFromShells, scheme.intoShells, 0)
        decayConfList   = Basics.generateConfigurationsWithElectronCapture(coreConfList, scheme.excitationFromShells, scheme.decayShells, 0)
        @show initialConfList
        @show coreConfList
        @show captureConfList
        @show decayConfList
        #
        # Determine first a hydrogenic spectrum for all subshells of the initial and doubly-excited states
        allConfList   = Configuration[];
        append!(allConfList, initialConfList);      append!(allConfList, captureConfList);      append!(allConfList, decayConfList)

        allSubshells  = Basics.extractRelativisticSubshellList(allConfList)
        primitives    = Bsplines.generatePrimitives(grid)
        orbitals      = Bsplines.generateOrbitalsHydrogenic(primitives, nm, allSubshells, printout=true)
        # Exclude configurations with too high mean energies
        en            = Float64[];
        for conf in initialConfList
            wen = Basics.determineMeanEnergy(conf, orbitals, nm, grid)
            push!(en, wen)
        end
        initialMean = sum(en) / length(en)
        println(">>> initial configuration(s) have mena energies  $initialMean  [a.u.].")
        #
        newCaptureConfList = Configuration[];  meanEnergies = Float64[]
        for  conf  in  captureConfList
            confMean = Basics.determineMeanEnergy(conf, orbitals, nm, grid)
            if  confMean - initialMean <= scheme.maxExcitationEnergy     push!(newCaptureConfList, conf)     end
        end

        return( (initialConfList, newCaptureConfList, decayConfList)  )
    end

    """
    `Cascade.perform(scheme::DielectronicCaptureScheme, comp::Cascade.Computation)`
        ... to set-up and perform a dielectronic-recombination (DR) plasma rate coefficient computation that combines the electron capture
            and the subsequent radiative stabilization steps. Such a computation starts from a given set of initial configurations xor
            initial multiplets and (1) generates all doubly-excited configurations due to the capture of an electron with a given maximum
            electron energy; (2) selects all electron capture (inverse Auger) and re-autoionization (Auger) steps and (3) selects
            all steps for radiative stabilization due to given parameters of the scheme::DielectronicCaptureScheme. The results of
            these DR plasma rate computation are comprised into (output) data::ExcitationData, while these data are only printed during
            the generation and nothing is returned.

    `Cascade.perform(scheme::DielectronicCaptureScheme, comp::Cascade.Computation; output=true, outputToFile::Bool=true)`
        ... to perform the same but to return the complete output in a dictionary that is written to disk and can be used in subsequent
            cascade simulation. The particular output depends on the specifications of the cascade.
    """
    function perform(scheme::DielectronicCaptureScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplets if initial configurations are given
        if  comp.initialConfigs != Configuration[]
            basis      = Basics.performSCF(comp.initialConfigs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplet  = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            ##x # Shift the initial level energy by -electronEnergyShift
            ##x if  scheme.electronEnergyShift != 0.
            ##x     multiplet = Basics.shiftTotalEnergies(multiplet, Defaults.convertUnits("energy: to atomic", -scheme.electronEnergyShift))
            ##x     println(">> Shift all initial level energies by -$(scheme.electronEnergyShift) $(Defaults.getDefaults("unit: energy"))")
            ##x end
            multiplets = [Multiplet("initial states", multiplet.levels)]
        else
            #== # Shift the initial level energy by -electronEnergyShift
            if  scheme.electronEnergyShift != 0.
                multiplets = Multiplet
                for  multiplet  in  comp.initialMultiplets
                    newMultiplet = Basics.shiftTotalEnergies(multiplet, Defaults.convertUnits("energy: to atomic", -scheme.electronEnergyShift))
                    push!(multiplets, newMultiplet)
                end
                println(">> Shift all initial level energies by -$(scheme.electronEnergyShift) $(Defaults.getDefaults("unit: energy"))")
            else
                multiplets = comp.initialMultiplets
            end  ==#
            multiplets = comp.initialMultiplets
        end
        # Print out initial configurations and levels
        Cascade.displayLevels(stdout, multiplets, sa="initial ")
        if  printSummary   Cascade.displayLevels(iostream, multiplets, sa="initial ")                            end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        wa  = Cascade.generateConfigurationsForDielectronicCapture(multiplets, comp.scheme, comp.nuclearModel, comp.grid)
        wb1 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[1], sa="initial configurations of the DR cascade ")
        wb2 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[2], sa="doubly-excited capture configurations of the DR cascade ")
        wb3 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[3], sa="decay configurations of the DR cascade ")
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        wc1 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb1)
        wc2 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb2, printout=false)
        wc3 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb3, printout=false)
        # Shift the initial level energy by -electronEnergyShift
        @show scheme.electronEnergyShift
        if  scheme.electronEnergyShift != 0.
            wc1old = wc1;   wc1 = Cascade.Block[]
            for  block in wc1old
                newMultiplet = Basics.shiftTotalEnergies(multiplet, Defaults.convertUnits("energy: to atomic", -scheme.electronEnergyShift))
                push!(wc1, Cascade.Block(block.NoElectrons, block.confs, block.hasMultiplet, newMultiplet))
            end
            println(">> Shift all initial level energies by $(-scheme.electronEnergyShift) $(Defaults.getDefaults("unit: energy"))")
        end
        #
        Cascade.displayBlocks(stdout, wc1, sa="from the initial configurations of the DR cascade ");
        Cascade.displayBlocks(stdout, wc2, sa="from the doubly-excited capture configurations of the DR cascade ")
        Cascade.displayBlocks(stdout, wc3, sa="from the decay configurations of the DR cascade ")
        if  printSummary   Cascade.displayBlocks(iostream, wc1, sa="from the initial configurations of the DR cascade ");
                           Cascade.displayBlocks(iostream, wc2, sa="from the doubly-excited capture configurations of the DR cascade ")
                           Cascade.displayBlocks(iostream, wc3, sa="from the decay configurations of the DR cascade ")            end
        #
        # Determine, modify and compute the transition data for all steps, ie. the PhotoIonization.Line's, etc.
        gMultiplets = Multiplet[];
        for block in wc1  push!(gMultiplets, block.multiplet)    end
        for block in wc2  push!(gMultiplets, block.multiplet)    end
        for block in wc3  push!(gMultiplets, block.multiplet)    end
        we = Cascade.determineSteps(scheme, comp, wc1, wc2, wc3)
        Cascade.displaySteps(stdout, we, sa="electron capture and stabilization ")
        if  printSummary   Cascade.displaySteps(iostream, we, sa="electron capture and stabilization ")    end
        wf   = Cascade.modifySteps(we)
        #
        data = Cascade.computeSteps(scheme, comp, wf)
        if output
            results = Base.merge( results, Dict("name"                          => comp.name) )
            results = Base.merge( results, Dict("cascade scheme"                => comp.scheme) )
            results = Base.merge( results, Dict("initial multiplets:"           => multiplets) )
            results = Base.merge( results, Dict("generated multiplets:"         => gMultiplets) )
            results = Base.merge( results, Dict("decay line data:"              => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "zzz-cascade-dr-rate-computations-" * string(Dates.now())[1:13] * ".jld"
            println("\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                    "\n   results = JLD.load(''$filename'')    ... to load the results back from file.")
            if  printSummary   println(iostream, "\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                                                 "\n   results = JLD.load(''$filename'')    ... to load the results back from file." )      end
            JLD2.@save filename results
        end
        ## return( results )
        return( results )
    end