text:
"""
Estimate a linear model with high dimensional categorical variables / instrumental variables

### Arguments
* `df::AbstractDataFrame`
* `model::Model`: A model created using [`@model`](@ref)
* `save::Union{Bool, Symbol} = false`: Should residuals and eventual estimated fixed effects saved in a dataframe? Use `save = :residuals` to only save residuals. Use `save = :fe` to only save fixed effects.
* `method::Symbol = :lsmr`: Method to deman regressors. `:lsmr` is akin to conjugate gradient descent.  With parallel use `:lsmr_parallel`. To use multi threaded use `lsmr_threads`. Other choices are `:qr` and `:cholesky` (factorization methods)
* `contrasts::Dict = Dict()` An optional Dict of contrast codings for each categorical variable in the `formula`.  Any unspecified variables will have `DummyCoding`.
* `maxiter::Integer = 10000`: Maximum number of iterations
* `tol::Real =1e-8`: Tolerance

### Details
Models with instruments variables are estimated using 2SLS. `reg` tests for weak instruments by computing the Kleibergen-Paap rk Wald F statistic, a generalization of the Cragg-Donald Wald F statistic for non i.i.d. errors. The statistic is similar to the one returned by the Stata command `ivreg2`.

### Examples
```julia
using DataFrames, RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
df.StateC =  categorical(df.State)
df.YearC =  categorical(df.Year)
reg(df, @model(Sales ~ Price, fe = StateC + YearC))
reg(df, @model(Sales ~ NDI, fe = StateC + StateC&Year))
reg(df, @model(Sales ~ NDI, fe = StateC*Year))
reg(df, @model(Sales ~ (Price ~ Pimin)))
reg(df, @model(Sales ~ Price, weights = Pop))
reg(df, @model(Sales ~ NDI, subset = State .< 30))
reg(df, @model(Sales ~ NDI, vcov = robust))
reg(df, @model(Sales ~ NDI, vcov = cluster(StateC)))
reg(df, @model(Sales ~ NDI, vcov = cluster(StateC + YearC)))
reg(df, @model(Sales ~ YearC), contrasts = Dict(:YearC => DummyCoding(base = 80)))
```
"""

function reg(df::AbstractDataFrame, m::Model; kwargs...)
    reg(df, m.f; m.dict..., kwargs...)
end

function reg(df::AbstractDataFrame, f::FormulaTerm;
    fe::Union{Symbol, Expr, Nothing} = nothing,
    vcov::Union{Symbol, Expr, Nothing} = :(simple()),
    weights::Union{Symbol, Expr, Nothing} = nothing,
    subset::Union{Symbol, Expr, Nothing} = nothing,
    maxiter::Integer = 10000, contrasts::Dict = Dict{Symbol, Any}(),
    tol::Real= 1e-8, dof_add::Integer = 0,
    save::Union{Bool, Symbol} = false,  method::Symbol = :lsmr, drop_singletons = true
   )

    if isa(vcov, Symbol)
        vcovformula = VcovFormula(Val{vcov})
    else
        vcovformula = VcovFormula(Val{vcov.args[1]}, (vcov.args[i] for i in 2:length(vcov.args))...)
    end

    ##############################################################################
    ##
    ## Parse formula
    ##
    ##############################################################################
    if  (ConstantTerm(0) ∉ eachterm(f.rhs)) & (ConstantTerm(1) ∉ eachterm(f.rhs))
        f = FormulaTerm(f.lhs, tuple(ConstantTerm(1), eachterm(f.rhs)...))
    end
    formula, formula_endo, formula_iv = decompose_iv(f)
    has_iv = formula_iv != nothing
    has_absorb = fe != nothing
    has_weights = weights != nothing

    ##############################################################################
    ##
    ## Save keyword argument
    ##
    ##############################################################################
    if !isa(save, Bool)
        if save ∉ (:residuals, :fe)
            error("the save keyword argument must be a Bool or a Symbol equal to :residuals or :fe")
        end
    end
    save_residuals = (save == :residuals) | (save == true)
    save_fe = (save == :fe) | ((save == true) & has_absorb)

    ##############################################################################
    ##
    ## Construct new dataframe after removing missing values
    ##
    ##############################################################################

    # create a dataframe without missing values & negative weights
    vars = unique(allvars(formula))
    iv_vars = unique(allvars(formula_iv))
    endo_vars = unique(allvars(formula_endo))
    absorb_vars = unique(allvars(fe))
    vcov_vars = unique(allvars(vcovformula))
    # create a dataframe without missing values & negative weights
    all_vars = unique(vcat(vars, vcov_vars, absorb_vars, endo_vars, iv_vars))

    esample = completecases(df, all_vars)

    if has_weights
        esample .&= isnaorneg(df[!, weights])
    end
    if subset != nothing
        subset = eval(evaluate_subset(df, subset))
        if length(subset) != size(df, 1)
            error("df has $(size(df, 1)) rows but the subset vector has $(length(subset)) elements")
        end
        esample .&= subset
    end

    if has_absorb
        feformula = @eval(@formula(0 ~ $(fe)))
        fes, ids = parse_fixedeffect(df, feformula)
        if drop_singletons
            for fe in fes
                drop_singletons!(esample, fe)
            end
        end
    end

    nobs = sum(esample)
    (nobs > 0) || error("sample is empty")

    # Compute weights
    sqrtw = get_weights(df, esample, weights)
    all(isfinite, sqrtw) || throw("Weights are not finite")

    # Compute pfe, a FixedEffectMatrix
    has_absorb_intercept = false
    if has_absorb
        # in case some FixedEffect does not have interaction, remove the intercept
        if any([isa(fe.interaction, Ones) for fe in fes])
            formula = FormulaTerm(formula.lhs, tuple(ConstantTerm(0), (t for t in eachterm(formula.rhs) if t!= ConstantTerm(1))...))
            has_absorb_intercept = true
        end
        fes = FixedEffect[_subset(fe, esample) for fe in fes]
        pfe = FixedEffectMatrix(fes, sqrtw, Val{method})
    end

    has_intercept = ConstantTerm(1) ∈ eachterm(formula.rhs)

    # Compute data for std errors
    vcov_method_data = VcovMethod(df[esample, vcov_vars], vcovformula)

    ##############################################################################
    ##
    ## Dataframe --> Matrix
    ##
    ##############################################################################
    subdf = columntable(df[esample, vars])
    formula_schema = apply_schema(formula, schema(formula, subdf, contrasts), StatisticalModel)
    # Obtain y
    # for a Vector{Float64}, conver(Vector{Float64}, y) aliases y
    y = convert(Vector{Float64}, response(formula_schema, subdf))
    all(isfinite, y) || throw("Some observations for the dependent variable are infinite")

    y .= y .* sqrtw

    # Obtain X
    Xexo = convert(Matrix{Float64}, modelmatrix(formula_schema, subdf))
    all(isfinite, Xexo) || throw("Some observations for the regressor are infinite")
    Xexo .= Xexo .* sqrtw

    yname, coef_names = coefnames(formula_schema)
    if !isa(coef_names, Vector)
        coef_names = [coef_names]
    end

    yname = Symbol(yname)
    coef_names = Symbol.(coef_names)

    if has_iv
        subdf = columntable(df[esample, endo_vars])
        formula_endo_schema = apply_schema(formula_endo, schema(formula_endo, subdf, contrasts), StatisticalModel)
        Xendo = convert(Matrix{Float64}, modelmatrix(formula_endo_schema, subdf))
        all(isfinite, Xendo) || throw("Some observations for the endogenous variable are infinite")
        Xendo .= Xendo .* sqrtw

        _, coefendo_names = coefnames(formula_endo_schema)
        if !isa(coefendo_names, Vector)
              coefendo_names = [coefendo_names]
          end
        append!(coef_names, Symbol.(coefendo_names))

        subdf = columntable(df[esample, iv_vars])
        formula_iv_schema = apply_schema(formula_iv, schema(formula_iv, subdf, contrasts), StatisticalModel)
        Z = convert(Matrix{Float64}, modelmatrix(formula_iv_schema, subdf))
        all(isfinite, Z) || throw("Some observations for the instrument are infinite")

        Z .= Z .* sqrtw

        if size(Z, 2) < size(Xendo, 2)
            error("Model not identified. There must be at least as many ivs as endogeneneous variables")
        end

        # modify formula to use in predict
        formula = FormulaTerm(formula.lhs, (tuple(eachterm(formula.rhs)..., eachterm(formula_endo.rhs)...)))
        formula_schema = apply_schema(formula, schema(formula, columntable(df), contrasts), StatisticalModel)
    end

    # compute tss now before potentially demeaning y
    tss = compute_tss(y, has_intercept | has_absorb_intercept, sqrtw)

    if has_absorb
        # used to compute tss even without save_fe
        if save_fe
            oldy = deepcopy(y)
            if has_iv
                oldX = hcat(Xexo, Xendo)
            else
                oldX = deepcopy(Xexo)
            end
        end

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

        y, b, c = solve_residuals!(y, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        Xexo, b, c = solve_residuals!(Xexo, pfe; maxiter = maxiter, tol = tol)
        append!(iterations, b)
        append!(convergeds, c)

        if has_iv
            Xendo, b, c = solve_residuals!(Xendo, pfe; maxiter = maxiter, tol = tol)
            append!(iterations, b)
            append!(convergeds, c)

            Z, b, c = solve_residuals!(Z, pfe; maxiter = maxiter, tol = tol)
            append!(iterations, b)
            append!(convergeds, c)
        end

        iterations = maximum(iterations)
        converged = all(convergeds)
        if converged == false
            @warn "convergence not achieved in $(iterations) iterations; try increasing maxiter or decreasing tol."
        end
    end

    ##############################################################################
    ##
    ## Get Linearly Independent Components of Matrix
    ##
    ##############################################################################

    # Compute linearly independent columns + create the Xhat matrix
    if has_iv
        # get linearly independent columns
        # note that I do it after residualizing
        baseall = basecol(Z, Xexo, Xendo)
        basecolXexo = baseall[(size(Z, 2)+1):(size(Z, 2) + size(Xexo, 2))]
        basecolXendo = baseall[(size(Z, 2) + size(Xexo, 2) + 1):end]
        Z = getcols(Z, baseall[1:size(Z, 2)])
        Xexo = getcols(Xexo, basecolXexo)
        Xendo = getcols(Xendo, basecolXendo)
        basecoef = vcat(basecolXexo, basecolXendo)

        # Build
        X = hcat(Xexo, Xendo)
        newZ = hcat(Xexo, Z)
        crossz = cholesky!(Symmetric(newZ' * newZ))
        Pi = crossz \ (newZ' * Xendo)
        Xhat = hcat(Xexo, newZ * Pi)

        # prepare residuals used for first stage F statistic
        ## partial out Xendo in place wrt (Xexo, Z)
        Xendo_res = gemm!('N', 'N', -1.0, newZ, Pi, 1.0, Xendo)
        ## partial out Z in place wrt Xexo
        Pi2 = cholesky!(Symmetric(Xexo' * Xexo)) \ (Xexo' * Z)
        Z_res = gemm!('N', 'N', -1.0, Xexo, Pi2, 1.0, Z)

    else
        # get linearly independent columns
        basecolXexo = basecol(Xexo)
        Xexo = getcols(Xexo, basecolXexo)
        Xhat = Xexo
        X = Xexo
        basecoef = basecolXexo
    end

    ##############################################################################
    ##
    ## Do the regression
    ##
    ##############################################################################

    crossx =  cholesky!(Symmetric(Xhat' * Xhat))
    coef = crossx \ (Xhat' * y)
    residuals = y - X * coef

    ##############################################################################
    ##
    ## Optionally save objects in a new dataframe
    ##
    ##############################################################################

    augmentdf = DataFrame()
    if save_residuals
        augmentdf.residuals =  Vector{Union{Missing, Float64}}(missing, length(esample))
        augmentdf[esample, :residuals] = residuals ./ sqrtw
    end
    if save_fe
        oldX = getcols(oldX, basecoef)
        newfes, b, c = solve_coefficients!(oldy - oldX * coef, pfe; tol = tol, maxiter = maxiter)
        for j in 1:length(fes)
            augmentdf[!, ids[j]] = Vector{Union{Float64, Missing}}(missing, length(esample))
            augmentdf[esample, ids[j]] = newfes[j]
        end
    end

    ##############################################################################
    ##
    ## Test Statistics
    ##
    ##############################################################################

    # Compute degrees of freedom
    dof_absorb = 0
    if has_absorb
        for fe in fes
            # adjust degree of freedom only if fe is not fully nested in a cluster variable:
            if isa(vcovformula, VcovClusterFormula) && any(isnested(fe, v) for v in eachcol(vcov_method_data.clusters))
                    dof_absorb += 1 # if fe is nested you still lose 1 degree of freedom
            else
                #only count groups that exists
                dof_absorb +=  length(Set(fe.refs))
            end
        end
    end
    dof_residual = max(1, nobs - size(X, 2) - dof_absorb - dof_add)

    # Compute rss, tss, r2, r2 adjusted
    rss = sum(abs2, residuals)
    mss = tss - rss
    r2 = 1 - rss / tss
    adjr2 = 1 - rss / tss * (nobs - (has_intercept | has_absorb_intercept)) / dof_residual
    if has_absorb
        r2_within = 1 - rss / compute_tss(y, (has_intercept | has_absorb_intercept), sqrtw)
    end

    # Compute standard error
    vcov_data = VcovData(Xhat, crossx, residuals, dof_residual)
    matrix_vcov = vcov!(vcov_method_data, vcov_data)

    # Compute Fstat
    (F, p) = compute_Fstat(coef, matrix_vcov, nobs, has_intercept, vcov_method_data, vcov_data)

    # Compute Fstat of First Stage
    if has_iv
        Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
        (F_kp, p_kp) = ranktest!(Xendo_res, Z_res, Pip,
                                  vcov_method_data, size(X, 2), dof_absorb)
    end

    # https://github.com/matthieugomez/FixedEffectModels.jl/issues/66
    dof_residual = df_FStat(vcov_method_data, vcov_data, has_intercept)

    ##############################################################################
    ##
    ## Return regression result
    ##
    ##############################################################################

    # add omitted variables
    if !all(basecoef)
        newcoef = fill(zero(Float64), length(basecoef))
        newmatrix_vcov = fill(NaN, (length(basecoef), length(basecoef)))
        newindex = [searchsortedfirst(cumsum(basecoef), i) for i in 1:length(coef)]
        for i in 1:length(coef)
            newcoef[newindex[i]] = coef[i]
            for j in 1:length(coef)
                newmatrix_vcov[newindex[i], newindex[j]] = matrix_vcov[i, j]
            end
        end
        coef = newcoef
        matrix_vcov = newmatrix_vcov
    end

    # return
    if !has_iv && !has_absorb
        return RegressionResult(coef, matrix_vcov, esample, augmentdf,
                                coef_names, yname, f, formula_schema, nobs, dof_residual,
                                rss, tss, r2, adjr2, F, p)
    elseif has_iv && !has_absorb
        return RegressionResultIV(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, formula_schema, nobs, dof_residual,
                                  rss, tss,  r2, adjr2, F, p, F_kp, p_kp)
    elseif !has_iv && has_absorb
        return RegressionResultFE(coef, matrix_vcov, esample, augmentdf,
                                  coef_names, yname, f, formula_schema, fe, nobs, dof_residual,
                                  rss, tss, r2, adjr2, r2_within, F, p, iterations, converged)
    elseif has_iv && has_absorb
        return RegressionResultFEIV(coef, matrix_vcov, esample, augmentdf,
                                   coef_names, yname, f, formula_schema, fe, nobs, dof_residual,
                                   rss, tss, r2, adjr2, r2_within, F, p, F_kp, p_kp,
                                   iterations, converged)
    end
end

##############################################################################
##
## Syntax without keywords
##
##############################################################################
function evaluate_subset(df, ex::Expr)
    if ex.head == :call
        return Expr(ex.head, ex.args[1], (evaluate_subset(df, ex.args[i]) for i in 2:length(ex.args))...)
    else
        return Expr(ex.head, (evaluate_subset(df, ex.args[i]) for i in 1:length(ex.args))...)
    end
end
evaluate_subset(df, ex::Symbol) = df[!, ex]
evaluate_subset(df, ex)  = ex


text:
fib(x::Int64) = x <= 1 ? 1 : fib(x-1) + fib(x - 2)

# added a test case for you to run and benchmark it.
n = 1:20
for i in n
  println(fib(i))
end

ext: jl

max_stars_repo_name: AliChraghi/LanguagesBenchmark

max_issues_repo_name: AliChraghi/LanguagesBenchmark

max_forks_repo_name: AliChraghi/LanguagesBenchmark


text:
include("readData.jl")
include("sparseSCVB0_lda.jl")
include("showTopics.jl")
using MAT
using JLD
using Base.DataFmt.writecsv
# this file shows a demo of sparseSCVB0 algorithm

# hyperparameters of the model
numTopics = 500;
minibatchSize = 20;
numDocumentIterations = 2500; # unnecessary variable. We used this one previously. But this variable is not needed in new implementation
burnIn = 1;
Alpha = 0.1;
Beta = 0.01;
tau = 1000.0;
kappa = 0.9;
tau2 = 10.0;
kappa2 = 0.9;
scale = 100.0;
sampleSize = 5;
sparsity = 1/numTopics

dataFilename = "data/NIPS.txt"; # one line per document, space-separated one-based dictionary indices for each consecutive word in document.
dictionaryFilename = "data/NIPSdict.txt"; # one line per dictionary word

nTrain = 1730; # must be divisible by numDocuments_miniEpoch
nTest = 10;
documents, test_documents, dictionary = readData(dataFilename, dictionaryFilename,nTrain,nTest);
numWords = length(dictionary);
numDocuments = length(documents);
# converting documents to sparse format
documents = wordVectorToSparseCounts(documents);

testTrain, testTest = splitTestSet(test_documents);
testTrain = wordVectorToSparseCounts(testTrain);
testTest = wordVectorToSparseCounts(testTest);
test_numDocuments = length(testTrain);

time_limit = 1.0*3600; #8.0*3600;

miniEpoch = 25;
numDocuments_miniEpoch = 1730;

# running sparseSCVB0 algorithm for LDA
sparseVBparams,saved_topics,saved_Totaltopics,saved_time,saved_iters,iter = sparseSCVB0_lda(documents,numWords,numDocuments,numTopics,
        numDocumentIterations,burnIn,minibatchSize,Alpha,Beta,dictionary,tau,kappa,tau2,kappa2,scale,sampleSize,time_limit,miniEpoch,numDocuments_miniEpoch,sparsity);

numToGet = 10; # number of top words to show as a topic
topWords,topicWordProbs = getImportantWordsInAllTopics(sparseVBparams.wordTopicCounts, dictionary, numTopics,numToGet);

# save & write all SCVB0 parameters in MAT file
file = matopen("sparseSCVB0_saved_topics.mat","w")
write(file,"saved_topics",saved_topics)
close(file)

file = matopen("sparseSCVB0_saved_Totaltopics.mat","w")
write(file,"saved_Totaltopics",saved_Totaltopics)
close(file)

file = matopen("sparseSCVB0_saved_time.mat","w")
write(file,"saved_time",saved_time)
close(file)

file = matopen("sparseSCVB0_saved_iters.mat","w")
write(file,"saved_iters",saved_iters)
close(file)

save("saved_params.jld","saved_topics",saved_topics,"saved_Totaltopics",saved_Totaltopics,"saved_time",saved_time,"saved_iters",saved_iters)

file = matopen("sparseSCVB0_topWords.mat","w")
write(file,"topWords",topWords)
close(file)

writecsv("sparseSCVB0_numIter.csv", iter);
# saving top words of all topics
writecsv("allTopics.csv", topWords);


text:
"""
`module  JAC.SpinAngular`
    ... a submodel of JAC that contains all methods for computing the spin-angular coefficients for one- and two-particle operators
        of given rank and for symmetry-adapted CSF. At the present, these spin-angular can be obtained for all scalar (and symmetric)
        one and two-particle interactions and they are sometimes refered to as 'pure angular coefficients' in order to distinguish
        them from those used in GRASP (and its derivatives) where part of the `physical interaction' were included orginal
        into the angular integrals. The spin-angular coefficients usally appear in the form:

            <CSF_l | operator | CSF_r> = sum_{t} coeff_t (L,abcd) * X^L (abcd)

        where CSF is a standard (jj-coupled) configuration state function, t a summation index, and where X^L(abcd) denotes the
        effective interaction strength.
"""
module SpinAngular

    using  Printf, ..AngularMomentum, ..Basics,  ..Defaults, ..ManyElectron, ..Radial

    export  Xronecker

    """
    `abstract type SpinAngular.AbstractAngularType`
        ... defines an abstract type and a number of data types to work with one- and two-particle operators of given rank, see also:

        + struct OneParticleOperator    ... to represent a single-particle operator with well-defined spherical tensor rank.
        + struct TwoParticleOperator    ... to represent a two-particle operator with well-defined spherical tensor rank.
    """
    abstract type  AbstractAngularType                           end

    """
    `struct  SpinAngular.OneParticleOperator  <:  AbstractAngularType`
        ... a struct for defining the spherial tensor structure of a single-particle operator.

        + rank            ::Int64             ... Rank of the operator.
        + parity          ::Basics.Parity     ... Parity of the operator (if needed ??)
        + sameOrbitalSet  ::Bool              ... True if orbitals for both CSF are taken from the same orbital set (if needed ??)
    """
    struct  OneParticleOperator  <:  AbstractAngularType
        rank              ::Int64
        parity            ::Basics.Parity
        sameOrbitalSet    ::Bool
    end

    """
    `SpinAngular.OneParticleOperator()`  ... constructor for setting the default values.
    """
    function OneParticleOperator()
    	OneParticleOperator( 0, Basics.plus, false)
    end

    # `Base.show(io::IO, op::OneParticleOperator)`  ... prepares a proper printout of the op::OneParticleOperator.
    function Base.show(io::IO, op::OneParticleOperator)
          sa = "One-particle operator O^($(op.rank) $(op.parity))"
          if op.sameOrbitalSet  sa = sa * ", defined for CSF from the same orbital set."
          else                  sa = sa * ", defined for CSF from different orbital sets."
          end
    	  println(io, sa)
    end

    """
    `struct  SpinAngular.TwoParticleOperator`
        ... a struct for defining the spherial tensor structure of a two-particle operator.

        + rank            ::Int64             ... Rank of the operator (if needed ??).
        + parity          ::Basics.Parity     ... Parity of the operator (if needed ??)
        + sameOrbitalSet  ::Bool              ... True if orbitals for both CSF are taken from the same orbital set (if needed ??)
    """
    struct  TwoParticleOperator
        rank              ::Int64
        parity            ::Basics.Parity
        sameOrbitalSet    ::Bool
    end

    """
    `SpinAngular.TwoParticleOperator()`  ... constructor for setting the default values.
    """
    function TwoParticleOperator()
    	TwoParticleOperator( 0, Basics.plus, false)
    end

    # `Base.show(io::IO, op::TwoParticleOperator)`  ... prepares a proper printout of the op::TwoParticleOperator.
    function Base.show(io::IO, op::TwoParticleOperator)
          sa = "Two-particle operator O^($(op.rank) $(op.parity))"
          if op.sameOrbitalSet  sa = sa * ", defined for CSF from the same orbital set."
          else                  sa = sa * ", defined for CSF from different orbital sets."
          end
    	  println(io, sa)
    end

    """
    `struct  SpinAngular.Coefficient1p`
        ... a struct for defining a single spin-angular coefficient for a reduced one-particle matrix element <a || o^(L) || b>.

        + L        ::Int64             ... Rank (L or nu) of the single-particle interaction strength.
        + a        ::Subshell          ... Left-hand subshell (orbital).
        + b        ::Subshell          ... Right-hand subshell (orbital).
        + v        ::Float64           ... (Value of) spin-angular coefficient.
    """
    struct  Coefficient1p
        L          ::Int64
        a          ::Subshell
        b          ::Subshell
        v          ::Float64
    end

    # `Base.show(io::IO, coeff::Coefficient1p)`  ... prepares a proper printout of the coeff::Coefficient1p.
    function Base.show(io::IO, coeff::Coefficient1p)
          sa = "\n V^($(coeff.L)) [$(coeff.a), $(coeff.b)] = $(coeff.v)"
    	  print(io, sa)
    end

    """
    `struct  SpinAngular.Coefficient2p`
        ... a struct for defining a single spin-angular coefficient for a reduced two-particle matrix element <ab || o^(L) || cd>,
            such as the Coulomb interaction strength.

        + L        ::Int64             ... Rank (L or nu) of the single-particle interaction strength.
        + a        ::Subshell          ... Left-hand subshell (orbital).
        + b        ::Subshell          ... Left-hand subshell (orbital).
        + c        ::Subshell          ... Right-hand subshell (orbital).
        + d        ::Subshell          ... Right-hand subshell (orbital).
        + v        ::Float64           ... (Value of) spin-angular coefficient.
    """
    struct  Coefficient2p
        L          ::Int64
        a          ::Subshell
        b          ::Subshell
        c          ::Subshell
        d          ::Subshell
        v          ::Float64
    end

    """
    `struct  SpinAngular.QspaceTerm`
        ... a struct for defining a subshell term/state  |j (nu) alpha Q J> == |j (nu) Q J Nr> for a subshell with well-defined j.

        + j          ::AngularJ64      ... subshell j
        + Q          ::AngularJ64      ... quasi-spin
        + J          ::AngularJ64      ... total J of subshell term
        + Nr         ::Int64           ... Additional quantum number Nr = 0,1,2.
        + index      ::Int64           ... Index in the QspaceTerm-list of given j
    """
    struct  QspaceTerm
        j            ::AngularJ64
        Q            ::AngularJ64
        J            ::AngularJ64
        Nr           ::Int64
        index        ::Int64
    end

    # `Base.show(io::IO, term::QspaceTerm)`  ... prepares a proper printout of term::QspaceTerm.
    function Base.show(io::IO, term::QspaceTerm)
        if  term.Nr == 0   sa = "|$(term.j) $(term.Q) $(term.J); index=$(term.index)>"
        else               sa = "|$(term.j) $(term.Q) $(term.J); Nr=$(term.Nr) index=$(term.index)>"
        end
        println(io, sa)
    end

    # `Base.show(io::IO, coeff::Coefficient2p)`  ... prepares a proper printout of the coeff::Coefficient2p.
    function Base.show(io::IO, coeff::Coefficient2p)
          sa = "\n V^($(coeff.L)) [$(coeff.a), $(coeff.b)| $(coeff.c), $(coeff.d)] = $(coeff.v)"
    	  print(io, sa)
    end

    """
    `@enum   SchemeEta`
        ... defines a enumeration for the allowed values of the Eta scheme to distinguish the coupling and selection of the
            tensor operators.

        + Eta_a         ... a^(qj)_mq
        + Eta_W         ... W^(k_12) = [a x a]
        + Eta_aW        ... W^(k_12, k_2) = [a1 x [a2 x a3]]
        + Eta_Wa        ... W^(k_12, k_2) = [[a1 x a2] x a3]
        + Eta_WW        ... W^(kk0) = [[a1 x a2]^k x [a3 x a4]^k]^0
    """
    @enum   SchemeEta    Eta_none   Eta_a   Eta_W   Eta_aW   Eta_Wa   Eta_WW

    """
    `struct  SpinAngular.SchemeGamma??`
        ... a singleton struct to distinguish between different coupling schemes of the matrix elements.
    """
    struct  SchemeGamma01   end
    struct  SchemeGamma02   end
    struct  SchemeGamma03   end

    """
    `struct  SpinAngular.SubshellTerm`
        ... a struct for defining a subshell term/state  |j^N (nu) alpha Q J> == |j^N (nu) Q J Nr> for a subshell with well-defined j.

        + j          ::AngularJ64      ... subshell j
        + Q          ::AngularJ64      ... quasi-spin
        + occupation ::Int64           ... occupation N
        + seniority  ::Int64           ... seniority
        + J          ::AngularJ64      ... total J of subshell term
        + Nr         ::Int64           ... Additional quantum number Nr = 0,1,2.
    """
    struct  SubshellTerm
        j            ::AngularJ64
        Q            ::AngularJ64
        occupation   ::Int64
        seniority    ::Int64
        J            ::AngularJ64
        Nr           ::Int64
    end

    # `Base.show(io::IO, term::SubshellTerm)`  ... prepares a proper printout of term::SubshellTerm.
    function Base.show(io::IO, term::SubshellTerm)
        if  term.Nr == 0   sa = "|$(term.j)^(term.occupation) ($(term.seniority)) $(term.Q) $(term.J)>"
        else               sa = "|$(term.j)^(term.occupation) ($(term.seniority)) $(term.Q) $(term.J); Nr=$(term.Nr)>"
        end
        println(io, sa)
    end

    """
    `struct  SpinAngular.Lambda`
        ... a struct for defining a subshell term/state  |j (nu) alpha Q J> == |j (nu) Q J Nr> for a subshell with well-defined j.

        + Ji         ::AngularJ64      ... J_i
        + Jj         ::AngularJ64      ... J_j
        + Jip        ::AngularJ64      ... J_i'
        + Jjp        ::AngularJ64      ... J_j'
    """
    struct  Lambda
        Ji           ::AngularJ64
        Jj           ::AngularJ64
        Jip          ::AngularJ64
        Jjp          ::AngularJ64
    end

    # `Base.show(io::IO, lambda::Lambda)`  ... prepares a proper printout of lambda::Lambda.
    function Base.show(io::IO, lambda::Lambda)
        sa = "Lambda($(lambda.Ji), $(lambda.Jj); lambda.Jip), $(lambda.Jjp))"
        println(io, sa)
    end

    include("module-SpinAngular-inc-reducedcoeffs.jl")

    #######################################################################################################################
    #######################################################################################################################

    """
    `SpinAngular.computeCoefficients(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})`
        ... computes the spin-angular coefficients for the reduced one-particle matrix element
            <leftCsf || op^(k) || rightCsf >  if both CSF refer to the same list of subshells; a list of one-particle
            coefficients coeffs::Array{Coefficient1p,1} is returned.
    """
    function  computeCoefficients(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})
        if      op.rank == 0   coeffs = computeCoefficientsScalar(op,    leftCsf, rightCsf, subshells)
        elseif  op.rank >  0   coeffs = computeCoefficientsNonScalar(op, leftCsf, rightCsf, subshells)
        else    error("stop a")
        end

        return (coeffs)
    end

    """
    `SpinAngular.computeCoefficients(op::SpinAngular.TwoParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})`
        ... computes the spin-angular coefficients for the reduced two-particle matrix element
            <leftCsf || op^(k) || rightCsf >   if both CSF refer to the same list of subshells; a Tuple of two lists
            with one- and two-particle coefficients
            tpl::Tuple{coeffs1p::Array{Coefficient1p,1}, coeffs2p::Array{Coefficient2p,1}}  is returned.
    """
    function  computeCoefficients(op::SpinAngular.TwoParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})
        if      op.rank == 0   coeffs = computeCoefficientsScalar(op, leftCsf, rightCsf, subshells)
        else    error("stop a")
        end

        return (coeffs)
    end

    """
    `SpinAngular.computeCoefficientsNonScalar(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})`
        ... computes the spin-angular coefficients for the reduced (scalar) one-particle matrix element
            <leftCsf || op^(0) || rightCsf >  if both CSF refer to the same list of subshells;
            a list of one-particle coefficients coeffs::Array{Coefficient1p,1} is returned.
    """
    function  computeCoefficientsNonScalar(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})
        coeffs1p = Coefficient1p[]
        # Cycle through all subshells of the bra- and ket-CSF
        for  (i,shi)  in  enumerate(subshells)
            for  (j,shj)  in  enumerate(subshells)

                push!(coeffs1p, Coefficient1p(op.rank, shi, shj, 0.))
            end
        end

        return( coeffs1p )
    end

    """
    `SpinAngular.computeCoefficientsScalar(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})`
        ... computes the spin-angular coefficients for the reduced (scalar) one-particle matrix element
            <leftCsf || op^(0) || rightCsf >  if both CSF refer to the same list of subshells;
            a list of one-particle coefficients coeffs::Array{Coefficient1p,1} is returned.
    """
    function  computeCoefficientsScalar(op::SpinAngular.OneParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})
        coeffs1p = Coefficient1p[]
        # Here, I'm following your notations in Eq. (6) in Gaigalas et al. (CPC, 2001)
        # Cycle through all subshells of the bra- and ket-CSF
        for  (i,shi)  in  enumerate(subshells)
            for  (j,shj)  in  enumerate(subshells)
                if  shi.kappa != shj.kappa      continue    end
                Ni  = leftCsf.occupation[i];   Nj  = rightCsf.occupation[j]
                ji  = Basics.subshell_j(shi);  jj  = Basics.subshell_j(shj)
                Ji  = leftCsf.subshellJ[i];    Jj  = rightCsf.subshellJ[j]
                si  = leftCsf.seniorityNr[i];    sj  = rightCsf.seniorityNr[j]
                Qi  = AngularJ64(1);           Qj  = AngularJ64(1)                          # !! This is yet not correct
                Nri = 0;                       Nrj = 0                                      # !! This is yet not correct
                iTerm     = SpinAngular.SubshellTerm(ji, Qi, si, Ji, Nri)
                jTerm     = SpinAngular.SubshellTerm(jj, Qj, sj, Jj, Nrj)
                braLambda = SpinAngular.Lambda(ji, ji, ji, ji)                              # !! This is yet not correct
                ketLambda = SpinAngular.Lambda(ji, ji, ji, ji)                              # !! This is yet not correct
                wa        = (-1)^SpinAngular.phaseDelta(op, Ni, Nj) * sqrt(Basics.twice(ji)+1) *
                            SpinAngular.recouplingMatrix(op, ji, jj, braLambda, ketLambda)
                #
                if  shi.n == shj.n    wb = SpinAngular.irreducibleTensor(iTerm, Ni, AngularM64(1//2), AngularM64(1//2), iTerm, Ni)
                else                  wb = 0.
                end
                if  shi.n != shj.n    wb = wb + SpinAngular.irreducibleTensor(iTerm, Ni, AngularM64(1//2),  jTerm, Nj) *
                                                SpinAngular.irreducibleTensor(iTerm, Ni, AngularM64(-1//2), jTerm, Nj)      end
                push!(coeffs1p, Coefficient1p(op.rank, shi, shj, wa*wb))
            end
        end

        return( coeffs1p )
    end

    """
    `SpinAngular.computeCoefficientsScalar(op::SpinAngular.TwoParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})`
        ... computes the spin-angular coefficients for the reduced (scalar) two-particle matrix element
            <leftCsf || op^(0) || rightCsf >; a Tuple of two lists with one- and two-particle coefficients
            tpl::Tuple{coeffs1p::Array{Coefficient1p,1}, coeffs2p::Array{Coefficient2p,1}}  is returned.
    """
    function  computeCoefficientsScalar(op::SpinAngular.TwoParticleOperator, leftCsf::CsfR, rightCsf::CsfR, subshells::Array{Subshell,1})
        coeffs1p   = Coefficient1p[];       coeffs2p    = Coefficient2p[]
        # Here, I'm following your notations in Eq. (13) in Gaigalas et al. (CPC, 2001)
        # Cycle through all subshells of the bra- and ket-CSF
        for  (i,shi)  in  enumerate(subshells)
            for  (j,shj)  in  enumerate(subshells)
                for  (ip,ship)  in  enumerate(subshells)
                    for  (jp,shjp)  in  enumerate(subshells)
                        schemeGamma = SpinAngular.determineSchemeGamma(shi, shj, ship, shjp)
                        schemeEta   = SpinAngular.determineSchemeEta(shi, shj, ship, shjp)

                        Ni  = leftCsf.occupation[i];   Nj  = leftCsf.occupation[j];   Nip = rightCsf.occupation[ip];  Njp  = rightCsf.occupation[jp]
                        ji  = Basics.subshell_j(shi);  jj  = Basics.subshell_j(shj);  jip = Basics.subshell_j(ship);  jjp  = Basics.subshell_j(shjp)
                        Ji  = leftCsf.subshellJ[i];    Jj  = leftCsf.subshellJ[j];    Jip = rightCsf.subshellJ[ip];   Jjp  = rightCsf.subshellJ[jp]
                        si  = leftCsf.seniorityNr[i];    sj  = leftCsf.seniorityNr[j];    sip = leftCsf.seniorityNr[ip];    sjp  = leftCsf.seniorityNr[jp]
                        Qi  = AngularJ64(1);           Qj  = AngularJ64(1);           Qip = AngularJ64(1);            Qjp  = AngularJ64(1)                          # !! This is yet not
                        Nri = 0;                       Nrj = 0;                       Nrip= 0;                        Nrjp= 0
                        # !! This is yet not correct
                        iTerm     = SpinAngular.SubshellTerm(ji, Qi, si, Ji, Nri)
                        jTerm     = SpinAngular.SubshellTerm(jj, Qj, sj, Jj, Nrj)
                        ipTerm    = SpinAngular.SubshellTerm(jip, Qip, sip, Jip, Nrip)
                        jpTerm    = SpinAngular.SubshellTerm(jjp, Qjp, sjp, Jjp, Nrjp)
                        #
                        wa = 0.
                        for  k12 = 0:10                                                                 # !! This is yet not correct
                            braLambda = SpinAngular.Lambda(ji, ji, ji, ji)                              # !! This is yet not correct
                            ketLambda = SpinAngular.Lambda(ji, ji, ji, ji)                              # !! This is yet not correct
                            k         = 1                                                               # !! This is yet not correct
                            wa        = wa + (-1)^SpinAngular.phaseDelta(op, Ni, Nj, Nip, Njp) *
                                        SpinAngular.coefficientThetap(schemeEta, k, shi, shj, ship, shjp) *
                                        SpinAngular.irreducibleTensorT(schemeEta, ji, jj, jip, jjp, braLambda, ketLambda, schemeGamma) *
                                        SpinAngular.recouplingMatrix(op, ji, jj, jip, jjp, braLambda, ketLambda, schemeGamma)
                        end
                        #
                        if  length(coeffs2p) < 7    push!(coeffs2p, Coefficient2p(op.rank, shi, shj, ship, shjp, wa))   end
                        #  This condition is of course arbitrary; from where do we get the one-particle coefficients
                        #  from kinetic energy, like for a Hamiltonian ? !!
                    end
                end
            end
        end

        return( coeffs2p )
    end

    """
    `SpinAngular.coefficientThetap(eta::SchemeEta, k::Int64, shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)`
        ... computes the coefficient Theta' (ni li ji, nj lj jj; ni' li' ji', nj lj jj', Eta) apart from the effective
            interaction strength that is set X^k (shi, shj, ship, shjp) = 1 in the procedure; a value::Float64 is returned.
    """
    function  coefficientThetap(eta::SchemeEta, k::Int64, shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)
        # Here likely call: AngularMomentum.
        if      eta == Eta_none
        elseif  eta == Eta_a
        elseif  eta == Eta_W
        else    error("stop a")
        end

        return( 1.0 )
    end

    """
    `SpinAngular.determineSchemeEta(shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)`
        ... determines the coupling scheme Eta for a set of subshells (shi, shj, ship, shjp);
            a scheme::SchemeEta is returned.
            !! Here, I do not really understand how you read off Eta = Eta_none, ... from the shells and or their occupation;
               this may require to provide further input to this procedure. !!
    """
    function  determineSchemeEta(shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)
        # Determine the correct scheme here
        return( Eta_none )
    end

    """
    `SpinAngular.determineSchemeGamma(shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)`
        ... determines the coupling scheme Gamma for a set of subshells (shi, shj, ship, shjp);
            a scheme::SchemeGamma?? is returned.
            !! Here, I do not really understand how you read off Gamma = 1, ..., 42  from the shells and or their occupation;
               this may require to provide further input to this procedure. !!
    """
    function  determineSchemeGamma(shi::Subshell, shj::Subshell, ship::Subshell, shjp::Subshell)
        # Determine the correct scheme here
        return( SchemeGamma01() )
    end

    """
    `SpinAngular.irreducibleTensor(term::SubshellTerm, N::Int64, mq::AngularM64, termp::SubshellTerm, Np::Int64)`
        ... computes the submatrix elements (j^N alpha Q J || a^(qj)_{m_q} || j'^N' alpha' Q' J')
            as defined by Eq. (7) in Gaigalas et al. (CPC, 2001); a value::Float64 is returned.
            !! We should check of whether we better include the occupation N into SubshellTerm; this would make the communication
               simpler. !!
    """
    function  irreducibleTensor(term::SubshellTerm, N::Int64, mq::AngularM64, termp::SubshellTerm, Np::Int64)
        # Here likely call: AngularMomentum.ClebschGordan, AngularMomentum.Wigner_3j

        return( 1.0 )
    end

    """
    `SpinAngular.irreducibleTensor(term::SubshellTerm, N::Int64, mLeft::AngularM64, mRight::AngularM64, termp::SubshellTerm, Np::Int64)`
        ... computes the submatrix elements (j^N alpha Q J || [a^(qj)_{m_q1} x a^(qj)_{m_q2}]^kj || j'^N' alpha' Q' J')
            as defined by Eq. (8) in Gaigalas et al. (CPC, 2001); a value::Float64 is returned.
            !! We should check of whether we better include the occupation N into SubshellTerm; this would make the communication
               simpler. !!
    """
    function  irreducibleTensor(term::SubshellTerm, N::Int64, mLeft::AngularM64, mRight::AngularM64, termp::SubshellTerm, Np::Int64)
        # Here likely call: AngularMomentum.ClebschGordan, AngularMomentum.Wigner_3j

        return( 1.0 )
    end

    """
    `SpinAngular.irreducibleTensorT(eta::SchemeEta, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                                    braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma01)`
        ... computes matrix elements T(Eta, ja, jb, ja', jb', Lambda^bra, Lambda^ket, Gamma=1) recouplings matrix for a two-particle
            operator; a value::Float64 is returned.
    """
    function  irreducibleTensorT(eta::SchemeEta, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                                 braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma01)
        # Here likely call functions:  AngularMomentum.
        if      eta == Eta_none
        elseif  eta == Eta_a
        elseif  eta == Eta_W
        else    error("stop a")
        end

        return( 1.0 )
    end

    """
    `SpinAngular.phaseDelta(op::SpinAngular.OneParticleOperator, Ni::Int64, Nj::Int64)`
        ... computes the phase Delta for a one-particle operator from the occupation number of the shells involved;
            a phase::Int64 is returned
    """
    function  phaseDelta(op::SpinAngular.OneParticleOperator, Ni::Int64, Nj::Int64)
        phase = 0                           # !! This is yet not correct
        return( phase )
    end

    """
    `SpinAngular.phaseDelta(op::SpinAngular.TwoParticleOperator, Ni::Int64, Nj::Int64, Nip::Int64, Njp::Int64)`
        ... computes the phase Delta for a two-particle operator from the occupation number of the shells involved;
            a phase::Int64 is returned
    """
    function  phaseDelta(op::SpinAngular.TwoParticleOperator, Ni::Int64, Nj::Int64, Nip::Int64, Njp::Int64)
        phase = 0                           # !! This is yet not correct
        return( phase )
    end

    """
    `SpinAngular.recouplingMatrix(op::SpinAngular.OneParticleOperator, ja::AngularJ64, jb::AngularJ64,
                                  braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda)`
        ... computes the recoupling matrix R(ja, jb, Lambda^bra, Lambda^ket) recouplings matrix for a one-particle
            operator; a value::Float64 is returned.
    """
    function  recouplingMatrix(op::SpinAngular.OneParticleOperator, ja::AngularJ64, jb::AngularJ64,
                               braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda)
        return( 1.0 )
    end

    """
    `SpinAngular.recouplingMatrix(op::SpinAngular.TwoParticleOperator, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                                  braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma01)`
        ... computes the recoupling matrix R(ja, jb, ja', jb', Lambda^bra, Lambda^ket, Gamma=1) recouplings matrix for a two-particle
            operator; a value::Float64 is returned.
    """
    function  recouplingMatrix(op::SpinAngular.TwoParticleOperator, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                               braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma01)
        return( 1.0 )
    end

    """
    `SpinAngular.recouplingMatrix(op::SpinAngular.TwoParticleOperator, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                                  braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma02)`
        ... computes the recoupling matrix R(ja, jb, ja', jb', Lambda^bra, Lambda^ket, Gamma=2) recouplings matrix for a two-particle
            operator; a value::Float64 is returned.
    """
    function  recouplingMatrix(op::SpinAngular.TwoParticleOperator, ja::AngularJ64, jb::AngularJ64, jap::AngularJ64, jbp::AngularJ64,
                               braLambda::SpinAngular.Lambda, ketLambda::SpinAngular.Lambda, scheme::SpinAngular.SchemeGamma02)
        return( 1.0 )
    end

end # module

#========================================================
Comments & thoughts on the implementation of this module
--------------------------------------------------------

++  I looked briefly through your ANO/RATIP implementation to remember the basic ingredients of such an implementation
    and to suggest some first functions which you likely need in 'this' or a rather similar form.
++  Overall, however, I belief that the overall code should be a factor 3-5 shorter and more easy to check.
++  Above you find the basic data structures which will allow me/JAC to communicate with this SpinAngular module.

++  In the first instance, I'm interested only in the spin-angular coefficients that belong to a single reduced
    many-electron matrix element ... I do currently not plan to benefit from the simultaneous computation of a whole
    raw or matrix of (matrix) elements.

++  I will be happy (and grateful) if we can agree about proper names of functions and if we make use of the
    standard Julia features; names like  anco_case_6 or similar are no longer appropriate and will make you own
    life difficult.
    Let's introduce a AbstractAngularCase if really necessary.

++  Please, use ? AngularMomentum.Wigner_6j(a, b, c, d, e, f) and similar for Wigner_3j and _9j symbols to get the Wigner symbols
    if needed.

++  Please, use CamelCase notation an start all variables with a lowercase letter; arrays are usually given the same name
    with an additional s:  coeff::Coefficient1p  vs.  coeffs::Array{Coefficient1p,1}
++  All angular momenta should be of type j::AngularJ64 (for j-values) and of type m::AngularM64 (for m-values,
    and if they appear at all in your implementation), please.
++  For a given subhsell sh::Subshell, get n and kappa simply by: sh.n, sh.kappa
++  For a given subhsell sh::Subshell, get l and j simply by:     Basics.subshell_l(sh), Basics.subshell_j(sh)

++  Look for ?CsfR ... to get the underlying definition of a CSF in JAC; it allows a direct access to J, parity,
    the occupation, seniority and all angular momenta from the coupling.

++  if csf.useStandardSubshells  subshells = Defaults.GBL_STANDARD_SUBSHELL_LIST
    else                         subshells = csf.subshells
    end
    In the latter case, these subshells could be different for the leftCsf and rightCsf.

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


text:
struct InvalidSignatureError <: Exception
end

struct MalformedJWTError <: Exception
    msg::String
end

struct NotSupportedJWTError <: Exception
    msg::String
end

Base.showerror(io::IO, e::InvalidSignatureError) = print(io, "Signature verification failed.")
Base.showerror(io::IO, e::MalformedJWTError) = print(io, "Malformed JWT: $(e.msg).")
Base.showerror(io::IO, e::NotSupportedJWTError) = print(io, "JWT format not supported: $(e.msg).")


text:
@testset "black box variational inference" begin

    Random.seed!(1)

    @gen function model()
        slope = @trace(normal(-1, exp(0.5)), :slope)
        intercept = @trace(normal(1, exp(2.0)), :intercept)
    end

    @gen function approx()
        @param slope_mu::Float64
        @param slope_log_std::Float64
        @param intercept_mu::Float64
        @param intercept_log_std::Float64
        @trace(normal(slope_mu, exp(slope_log_std)), :slope)
        @trace(normal(intercept_mu, exp(intercept_log_std)), :intercept)
    end

    # to regular black box variational inference
    init_param!(approx, :slope_mu, 0.)
    init_param!(approx, :slope_log_std, 0.)
    init_param!(approx, :intercept_mu, 0.)
    init_param!(approx, :intercept_log_std, 0.)

    observations = choicemap()
    update = ParamUpdate(GradientDescent(1, 100000), approx)
    update = ParamUpdate(GradientDescent(1., 1000), approx)
    black_box_vi!(model, (), observations, approx, (), update;
        iters=2000, samples_per_iter=100, verbose=false)
    slope_mu = get_param(approx, :slope_mu)
    slope_log_std = get_param(approx, :slope_log_std)
    intercept_mu = get_param(approx, :intercept_mu)
    intercept_log_std = get_param(approx, :intercept_log_std)
    @test isapprox(slope_mu, -1., atol=0.001)
    @test isapprox(slope_log_std, 0.5, atol=0.001)
    @test isapprox(intercept_mu, 1., atol=0.001)
    @test isapprox(intercept_log_std, 2.0, atol=0.001)

    # smoke test for black box variational inference with Monte Carlo objectives
    init_param!(approx, :slope_mu, 0.)
    init_param!(approx, :slope_log_std, 0.)
    init_param!(approx, :intercept_mu, 0.)
    init_param!(approx, :intercept_log_std, 0.)
    black_box_vimco!(model, (), observations, approx, (), update, 20;
        iters=50, samples_per_iter=100, verbose=false, geometric=false)

    init_param!(approx, :slope_mu, 0.)
    init_param!(approx, :slope_log_std, 0.)
    init_param!(approx, :intercept_mu, 0.)
    init_param!(approx, :intercept_log_std, 0.)
    black_box_vimco!(model, (), observations, approx, (), update, 20;
        iters=50, samples_per_iter=100, verbose=false, geometric=true)

end

@testset "minimal vae" begin

    Random.seed!(1)

    n = 10

    @gen function model()
        @param theta::Float64
        for i in 1:n
            z = ({(:z, i)} ~ normal(theta, 1))
            {(:x, i)} ~ normal(z, 1)
        end
    end

    @gen function approx(xs)
        @param mu_coeffs::Vector{Float64} # 2 x 1; should be [opt_theta / 2, 0.5]
        @param log_std::Float64
        for i in 1:n
            mu = mu_coeffs[1] + mu_coeffs[2] * xs[i]
            {(:z, i)} ~ normal(mu, exp(log_std))
        end
    end

    observations = choicemap()
    xs = Float64[]
    for i in 1:n
        z = normal(-1.0, 1)
        x = normal(z, 1)
        push!(xs, x)
        observations[(:x, i)] = x
    end

    # the optimum maximum marginal likelihood is at theta = mean(xs)
    opt_theta = sum(xs) / length(xs)
    println("opt_theta: $opt_theta")
    # see https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf
    posterior_precisions = 1.0 + 1 * 1.0 # n = 1
    posterior_means = (xs .+ opt_theta) ./ 2.0
    posterior_precisions = 2.0
    @gen function optimum_approx()
        for i in 1:n
            {(:z, i)} ~ normal(posterior_means[i], sqrt(1.0 / posterior_precisions))
        end
    end
    init_param!(model, :theta, opt_theta)
    approx_trace = simulate(optimum_approx, ())
    (model_trace, _) = generate(model, (), merge(get_choices(approx_trace), observations))
    # note that p(z1..zn, x1..xn) / p(z1..zn | x1..xn) = p(x1...xn) - for all z1..zn
    log_marginal_likelihood = get_score(model_trace) - get_score(approx_trace)
    println("true optimum log_marginal_likelihood: $log_marginal_likelihood")

    # using BBVI with score function estimator
    init_param!(model, :theta, 0.0)
    init_param!(approx, :mu_coeffs, zeros(2))
    init_param!(approx, :log_std, 0.0)
    approx_update = ParamUpdate(FixedStepGradientDescent(0.0001), approx)
    model_update = ParamUpdate(FixedStepGradientDescent(0.002), model)
    @time (_, _, elbo_history, _) =
        black_box_vi!(model, (), model_update, observations,
                      approx, (xs,), approx_update;
                      iters=3000, samples_per_iter=20, verbose=false)
    @test isapprox(get_param(model, :theta), opt_theta, atol=1e-1)
    println("final theta: $(get_param(model, :theta))")
    println("final elbo estimate: $(elbo_history[end])")
    @test isapprox(elbo_history[end], log_marginal_likelihood, rtol=0.1)

    # using VIMCO
    init_param!(model, :theta, 0.0)
    init_param!(approx, :mu_coeffs, zeros(2))
    init_param!(approx, :log_std, 0.0)
    approx_update = ParamUpdate(FixedStepGradientDescent(0.001), approx)
    model_update = ParamUpdate(FixedStepGradientDescent(0.01), model)
    @time (_, _, elbo_history, _) =
        black_box_vimco!(model, (), model_update, observations,
                         approx, (xs,), approx_update, 10;
                         iters=1000, samples_per_iter=10, verbose=false)
    println("final theta: $(get_param(model, :theta))")
    println("final elbo estimate: $(elbo_history[end])")
    @test isapprox(elbo_history[end], log_marginal_likelihood, rtol=0.1)
end


text:
"""
    ADVector{T}

A cache of vectors for making it easier to ForwardDiff. A vector of the desired type
    is extracted by indexing:

    vec = adv[eltype(x)]

Since ForwardDiff Duals are parameterized on the chunk size, this type will store any
number of vectors for different chunk sizes. A new vector of duals will be automatically
created (once) if it doesn't exist yet.

For best performance (by avoiding type-instability), annotate the indexing operation:

    vec = adv[T]::Vector{T}

While the annotation could have been included in the indexing operation itself, some
benchmark results found it significantly faster to apply the annotation at the call site.

# Constructor

    ADVector{T}(n)

where `n` is the size of the vectors, and `T` is the numeric type (typically `Float64`).
"""
struct ADVector{T}
    v::Vector{T}
    d::Dict{DataType,Vector{D} where {D<:ForwardDiff.Dual{Nothing,T}}}
    function ADVector{T}(n::Integer) where {T}
        new{T}(zeros(T, n), Dict{DataType,Vector{D} where D}())
    end
end
Base.getindex(adv::ADVector, ::Type{<:Number}) = adv.v
function Base.getindex(adv::ADVector, T::Type{<:ForwardDiff.Dual})
    if haskey(adv.d, T)
        return adv.d[T]
    else
        d = zeros(T, length(adv.v))
        adv.d[T] = d
        return d
    end
end

########################################
# Explicit Methods
########################################

struct Euler <: Explicit
    Euler(n::Integer, m::Integer) = new()
end

function integrate(::Euler, model, x, u, t, h)
    xdot = dynamics(model, x, u, t)
    return x + h * xdot
end

function integrate!(int::Euler, model, xn, x, u, t, h)
    dynamics!(model, xn, x, u, t)
    xn .*= h
    xn .+= x
    return nothing
end

function jacobian!(int::Euler, sig::FunctionSignature, model::ContinuousDynamics, J, xn, x, u, t, h)
    # Call the user-defined Continuous-time Jacobian
    jacobian!(model, J, xn, x, u, t)
    J .*= h
    for i = 1:state_dim(model)
        J[i, i] += 1.0
    end
    return nothing
end

struct RK3 <: Explicit
    k1::ADVector{Float64}
    k2::ADVector{Float64}
    k3::ADVector{Float64}
    A::Vector{Matrix{Float64}}
    B::Vector{Matrix{Float64}}
    dA::Vector{Matrix{Float64}}
    dB::Vector{Matrix{Float64}}
    function RK3(n::Integer, m::Integer)
        k1, k2 = ADVector{Float64}(n), ADVector{Float64}(n)
        k3 = ADVector{Float64}(n)
        A = [zeros(n, n) for i = 1:3]
        B = [zeros(n, m) for i = 1:3]
        dA = [zeros(n, n) for i = 1:3]
        dB = [zeros(n, m) for i = 1:3]
        new(k1, k2, k3, A, B, dA, dB)
    end
end
getks(int::RK3, ::Type{T}) where {T} =
    int.k1[T]::Vector{T}, int.k2[T]::Vector{T}, int.k3[T]::Vector{T}

function integrate(::RK3, model, x, u, t, h)
    k1 = dynamics(model, x,            u, t      ) * h
    k2 = dynamics(model, x + k1 / 2,   u, t + h/2) * h
    k3 = dynamics(model, x - k1 + 2k2, u, t + h  ) * h
    return x + (k1 + 4k2 + k3) / 6
end

function integrate!(int::RK3, model, xn, x, u, t, h)
    T = eltype(xn)
    k1, k2, k3 = getks(int, T)
    dynamics!(model, k1, x, u, t)
    @. xn = x + k1 * h / 2
    dynamics!(model, k2, xn, u, t + h / 2)
    @. xn = x - k1 * h + 2 * k2 * h
    dynamics!(model, k3, xn, u, t + h)
    @. xn = x + h * (k1 + 4k2 + k3) / 6
    return nothing
end

function jacobian!(int::RK3, sig::StaticReturn, model, J, xn, x, u, t, h)
    n, m = size(model)
    ix = SVector{n}(1:n)
    iu = SVector{m}(n+1:n+m)
    k1 = dynamics(model, x,            u, t      ) * h
    k2 = dynamics(model, x + k1 / 2,   u, t + h/2) * h

    jacobian!(model, J, xn, x, u, t)
    A1, B1 = J[ix, ix], J[ix, iu]

    jacobian!(model, J, xn, x + k1 / 2, u, t + h / 2)
    A2, B2 = J[ix, ix], J[ix, iu]

    jacobian!(model, J, xn, x - k1 + 2k2, u, t + h)
    A3, B3 = J[ix, ix], J[ix, iu]

    dA1 = A1 * h
    dA2 = A2 * (I + 0.5 * dA1) * h
    dA3 = A3 * (I - dA1 + 2 * dA2) * h

    dB1 = B1 * h
    dB2 = B2 * h + 0.5 * A2 * dB1 * h
    dB3 = B3 * h + A3 * (2dB2 - dB1) * h

    J[ix, ix] .= I + (dA1 + 4dA2 + dA3) / 6
    J[ix, iu] .= (dB1 + 4dB2 + dB3) / 6

    return nothing
end

function jacobian!(int::RK3, sig::InPlace, model, J, xn, x, u, t, h)
    # x,u,t,h = state(z), control(z), time(z), timestep(z)
    k1, k2, k3 = getks(int, Float64)
    A1, A2, A3 = int.A[1], int.A[2], int.A[3]
    B1, B2, B3 = int.B[1], int.B[2], int.B[3]
    dA1, dA2, dA3 = int.dA[1], int.dA[2], int.dA[3]
    dB1, dB2, dB3 = int.dB[1], int.dB[2], int.dB[3]
    n, m = size(model)
    ix, iu = 1:n, n+1:n+m

    jacobian!(model, J, k1, x, u, t)
    dynamics!(model,    k1, x, u, t)
    A1 .= @view J[ix, ix]
    B1 .= @view J[ix, iu]

    @. xn = x + k1 * h / 2
    jacobian!(model, J, k2, xn, u, t + h / 2)
    dynamics!(model,    k2, xn, u, t + h / 2)
    A2 .= @view J[ix, ix]
    B2 .= @view J[ix, iu]

    @. xn = x - k1 * h + 2 * k2 * h
    jacobian!(model, J, k3, xn, u, t + h)
    dynamics!(model,    k3, xn, u, t + h)
    A3 .= @view J[ix, ix]
    B3 .= @view J[ix, iu]

    # dA = A1 * h
    dA1 .= A1 .* h

    # dA2 = A2 * (I + 0.5 * dA1) * h
    mul!(dA2, A2, dA1, 0.5, 0.0)
    dA2 .+= A2
    dA2 .*= h

    # dA3 = A3 * (I - dA1 + 2 * dA2) * h
    mul!(dA3, A3, dA2, 2.0, 0.0)
    mul!(dA3, A3, dA1,-1.0, 1.0)
    dA3 .+= A3
    dA3 .*= h

    # dB1 = B1 * h
    dB1 .= B1 .* h

    # dB2 = B2 * h + 0.5 * A2 * dB1 * h
    dB2 .= B2
    mul!(dB2, A2, dB1, 0.5, 1.0)
    dB2 .*= h

    # dB3 = B3 * h + A3 * (2dB2 - dB1) * h
    dB3 .= B3
    mul!(dB3, A3, dB2, 2.0, 1.0)
    mul!(dB3, A3, dB1,-1.0, 1.0)
    dB3 .*= h

    @. J[ix, ix] = (dA1 + 4dA2 + dA3) / 6
    for i = 1:n
        J[i, i] += 1.0
    end
    @. J[ix, iu] = (dB1 + 4dB2 + dB3) / 6

    return nothing
end

"Fourth-order Runge-Kutta method with zero-order-old on the controls"
struct RK4 <: Explicit
    k1::ADVector{Float64}
    k2::ADVector{Float64}
    k3::ADVector{Float64}
    k4::ADVector{Float64}
    A::Vector{Matrix{Float64}}
    B::Vector{Matrix{Float64}}
    dA::Vector{Matrix{Float64}}
    dB::Vector{Matrix{Float64}}
    function RK4(n::Integer, m::Integer)
        k1, k2 = ADVector{Float64}(n), ADVector{Float64}(n)
        k3, k4 = ADVector{Float64}(n), ADVector{Float64}(n)
        A = [zeros(n, n) for i = 1:4]
        B = [zeros(n, m) for i = 1:4]
        dA = [zeros(n, n) for i = 1:4]
        dB = [zeros(n, m) for i = 1:4]
        new(k1, k2, k3, k4, A, B, dA, dB)
    end
end
getks(int::RK4, ::Type{T}) where {T} =
    int.k1[T]::Vector{T}, int.k2[T]::Vector{T}, int.k3[T]::Vector{T}, int.k4[T]::Vector{T}

function integrate(::RK4, model, x, u, t, h)
    k1 = dynamics(model, x, u, t) * h
    k2 = dynamics(model, x + k1 / 2, u, t + h / 2) * h
    k3 = dynamics(model, x + k2 / 2, u, t + h / 2) * h
    k4 = dynamics(model, x + k3, u, t + h) * h
    x + (k1 + 2k2 + 2k3 + k4) / 6
end

function integrate!(int::RK4, model, xn, x, u, t, h)
    T = eltype(xn)
    k1, k2, k3, k4 = getks(int, T)
    dynamics!(model, k1, x, u, t)
    @. xn = x + k1 * h / 2
    dynamics!(model, k2, xn, u, t + h / 2)
    @. xn = x + k2 * h / 2
    dynamics!(model, k3, xn, u, t + h / 2)
    @. xn = x + k3 * h
    dynamics!(model, k4, xn, u, t + h)
    @. xn = x + h * (k1 + 2k2 + 2k3 + k4) / 6
    return nothing
end

function jacobian!(int::RK4, sig::StaticReturn, model, J, xn, x, u, t, h)
    n, m = size(model)
    ix = SVector{n}(1:n)
    iu = SVector{m}(n+1:n+m)
    k1 = dynamics(model, x, u, t) * h
    k2 = dynamics(model, x + k1 / 2, u, t + h / 2) * h
    k3 = dynamics(model, x + k2 / 2, u, t + h / 2) * h

    jacobian!(model, J, xn, x, u, t)
    A1, B1 = J[ix, ix], J[ix, iu]

    jacobian!(model, J, xn, x + k1 / 2, u, t + h / 2)
    A2, B2 = J[ix, ix], J[ix, iu]

    jacobian!(model, J, xn, x + k2 / 2, u, t + h / 2)
    A3, B3 = J[ix, ix], J[ix, iu]

    jacobian!(model, J, xn, x + k3, u, t + h)
    A4, B4 = J[ix, ix], J[ix, iu]

    dA1 = A1 * h
    dA2 = A2 * (I + 0.5 * dA1) * h
    dA3 = A3 * (I + 0.5 * dA2) * h
    dA4 = A4 * (I + dA3) * h

    dB1 = B1 * h
    dB2 = B2 * h + 0.5 * A2 * dB1 * h
    dB3 = B3 * h + 0.5 * A3 * dB2 * h
    dB4 = B4 * h + A4 * dB3 * h

    J[ix, ix] .= I + (dA1 + 2dA2 + 2dA3 + dA4) / 6
    J[ix, iu] .= (dB1 + 2dB2 + 2dB3 + dB4) / 6

    return nothing
end

function jacobian!(int::RK4, sig::InPlace, model, J, xn, x, u, t, h)
    # x,u,t,h = state(z), control(z), time(z), timestep(z)
    k1, k2, k3, k4 = getks(int, Float64)
    A1, A2, A3, A4 = int.A[1], int.A[2], int.A[3], int.A[4]
    B1, B2, B3, B4 = int.B[1], int.B[2], int.B[3], int.B[4]
    dA1, dA2, dA3, dA4 = int.dA[1], int.dA[2], int.dA[3], int.dA[4]
    dB1, dB2, dB3, dB4 = int.dB[1], int.dB[2], int.dB[3], int.dB[4]
    n, m = size(model)
    ix, iu = 1:n, n+1:n+m

    jacobian!(model, J, k1, x, u, t)
    dynamics!(model,    k1, x, u, t)
    A1 .= @view J[ix, ix]
    B1 .= @view J[ix, iu]

    @. xn = x + k1 * h / 2
    jacobian!(model, J, k2, xn, u, t + h / 2)
    dynamics!(model,    k2, xn, u, t + h / 2)
    A2 .= @view J[ix, ix]
    B2 .= @view J[ix, iu]

    @. xn = x + k2 * h / 2
    jacobian!(model, J, k3, xn, u, t + h / 2)
    dynamics!(model,    k3, xn, u, t + h / 2)
    A3 .= @view J[ix, ix]
    B3 .= @view J[ix, iu]

    @. xn = x + k3 * h
    jacobian!(model, J, k4, xn, u, t + h)
    A4 .= @view J[ix, ix]
    B4 .= @view J[ix, iu]

    # dA = A1 * h
    dA1 .= A1 .* h

    # dA2 = A2 * (I + 0.5 * dA1) * h
    mul!(dA2, A2, dA1, 0.5, 0.0)
    dA2 .+= A2
    dA2 .*= h

    # dA3 = A3 * (I + 0.5 * dA2) * h
    mul!(dA3, A3, dA2, 0.5, 0.0)
    dA3 .+= A3
    dA3 .*= h

    # dA4 = A4 * (I + dA3) * h
    mul!(dA4, A4, dA3, 1.0, 0.0)
    dA4 .+= A4
    dA4 .*= h

    # dB1 = B1 * h
    dB1 .= B1 .* h

    # dB2 = B2 * h + 0.5 * A2 * dB1 * h
    dB2 .= B2
    mul!(dB2, A2, dB1, 0.5, 1.0)
    dB2 .*= h

    # dB3 = B3 * h + 0.5 * A3 * dB2 * h
    dB3 .= B3
    mul!(dB3, A3, dB2, 0.5, 1.0)
    dB3 .*= h

    # dB4 = B4 * h + A4 * dB3 * h
    mul!(dB4, A4, dB3)
    dB4 .+= B4
    dB4 .*= h

    @. J[ix, ix] = (dA1 + 2dA2 + 2dA3 + dA4) / 6
    for i = 1:n
        J[i, i] += 1.0
    end
    @. J[ix, iu] = (dB1 + 2dB2 + 2dB3 + dB4) / 6

    return nothing
end

########################################
# Implicit Methods
########################################

struct ImplicitMidpoint <: Implicit
    xmid::ADVector{Float64}
    function ImplicitMidpoint(n::Integer, m::Integer)
        new(ADVector{Float64}(n))
    end
end

function dynamics_error(
    ::ImplicitMidpoint,
    model::ContinuousDynamics,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    t, h = time(z1), timestep(z1)
    x1, u1 = state(z1), control(z1)
    x2, u2 = state(z2), control(z2)
    xmid = (x1 + x2) / 2
    fmid = dynamics(model, xmid, u1, t + h / 2)
    x1 + h * fmid - x2
end

function dynamics_error!(
    int::ImplicitMidpoint,
    model::ContinuousDynamics,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    T = eltype(y2)
    t, h = time(z1), timestep(z1)
    x1, u1 = state(z1), control(z1)
    x2, u2 = state(z2), control(z2)
    xmid = int.xmid[T]::Vector{T}

    @. xmid = (x1 + x2) / 2
    dynamics!(model, y2, xmid, u1, t + h / 2)
    @. y2 = x1 + h * y2 - x2
    return nothing
end

function dynamics_error_jacobian!(
    int::ImplicitMidpoint,
    ::StaticReturn,
    model::ContinuousDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    n, m = size(model)
    ix = SVector{n}(1:n)
    iu = SVector{m}(n+1:n+m)
    t, h = time(z1), timestep(z1)
    x1, u1 = state(z1), control(z1)
    x2, u2 = state(z2), control(z2)

    xmid = (x1 + x2) / 2
    jacobian!(model, J1, y1, xmid, u1, t)

    J1 .*= h
    J1[ix,ix] ./= 2
    J2[ix,ix] .= J1[ix,ix]
    J2[ix,iu] .= 0
    for i = 1:n
        J1[i,i] += 1.0
        J2[i,i] -= 1.0
    end
    return nothing
end

function dynamics_error_jacobian!(
    int::ImplicitMidpoint,
    ::InPlace,
    model::ContinuousDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    n, m = size(model)
    ix = 1:n
    iu = n+1:n+m
    t, h = time(z1), timestep(z1)
    x1, u1 = state(z1), control(z1)
    x2, u2 = state(z2), control(z2)
    xmid = int.xmid[Float64]

    @. xmid = (x1 + x2) / 2
    jacobian!(model, J1, y1, xmid, u1, t)
    J1 .*= h
    A1 = view(J1, ix, ix)
    A2 = view(J2, ix, ix)
    A1 ./= 2
    A2 .= A1
    view(J2, ix, iu) .= 0
    for i = 1:n
        J1[i,i] += 1.0
        J2[i,i] -= 1.0
    end
    return nothing
end

a = 1


text:
using TensorOperations
using Lints
using Fermi.DIIS
using Fermi.Integrals: IntegralHelper
using Fermi.Orbitals: CanonicalOrbital, CanonicalOrbitals, OrbDict

# Define Algorithims
abstract type RHFAlgorithm end
struct ConventionalRHF <: RHFAlgorithm end
struct DFRHF <: RHFAlgorithm end

# Define Guesses
abstract type RHFGuess end
struct CoreGuess   <: RHFGuess end
struct GWHGuess    <: RHFGuess end
struct HuckelGuess <: RHFGuess end

"""
    Fermi.HartreeFock.RHF

Wave function object for Restricted Hartree-Fock methods

# High Level Interface
    RHF()

Computes RHF using information from Fermi.CurrentOptions.

# Fields:

    molecule    Molecule object used to compute the RHF wave function
    energy      RHF Energy
    ndocc       Number of doubly occupied spatial orbitals
    nvir        Number of virtual spatial orbitals
    C           Array with MO coefficients
    eps         Array with MO energies

# Relevant options

These options can be set with `@set <option> <value>`

| Option         | What it does                      | Type      | choices [default]     |
|----------------|-----------------------------------|-----------|-----------------------|
| `scf_alg`      | picks SCF algorithm               | `String`  | "df" ["conventional"] |
| `scf_max_rms`  | RMS density convergence criterion |`Float64`  | [10^-10]              |
| `scf_max_iter` | Max number of iterations          | `Int`     | [50]                  |
| `basis`        | What basis set to use             | `String`  | ["sto-3g"]            |
| `jkfit`        | What aux. basis set to use for JK | `String`  | ["auto"]              |
| `oda`          | Whether to use ODA                | `Bool`    | [`True`]              |
| `oda_cutoff`   | When to turn ODA off (RMS)        | `Float64` | [1E-1]                |
| `oda_shutoff`  | When to turn ODA off (iter)       | `Int`     | [20]                  |
| `scf_guess`    | Which guess density to use        |           | "core" ["gwh"]        |

# Lower level interfaces

    RHF(molecule::Molecule, aoint::IntegralHelper, C::Array{Float64,2}, ERI::Array{Float64,N}, Λ::Array{Float64,2}) where N

The RHF kernel. Computes RHF on the given `molecule` with integral information defined in `aoint`. Starts from
the given C matrix.

_struct tree:_

**RHF** <: AbstractHFWavefunction <: AbstractReferenceWavefunction <: AbstractWavefunction
"""
struct RHF <: AbstractHFWavefunction
    molecule::Molecule
    energy::Float64
    ndocc::Int
    nvir::Int
    eps::Array{Float64,1}
    ints::IntegralHelper
end

function select_alg(A::String)
    implemented = Dict{String,Any}(
        "conventional" => (ConventionalRHF()),
        "df"           => (DFRHF())
       )

    try
        return implemented[A]
    catch KeyError
        throw(Fermi.InvalidFermiOption("Invalid RHF algorithm: $(A)"))
    end
end

function select_guess(A::String)
    implemented = Dict{String,Any}(
        "gwh"  => GWHGuess(),
        "core" => CoreGuess()
    )
    try
        return implemented[A]
    catch KeyError
        throw(Fermi.InvalidFermiOption("Invalid RHF guess: $(A)"))
    end
end

function RHF()
    molecule = Molecule()
    RHF(molecule)
end

function RHF(molecule::Molecule)
    Alg = select_alg(Fermi.CurrentOptions["scf_alg"])
    ints = Fermi.Integrals.IntegralHelper()
    guess = select_guess(Fermi.CurrentOptions["scf_guess"])
    RHF(molecule, ints, Alg, guess)
end

function RHF(molecule::Molecule, aoint::IntegralHelper, Alg::B, guess::GWHGuess) where B <: RHFAlgorithm

    #form GWH guess
    @output "Using GWH Guess\n"
    S = Array(Hermitian(aoint["S"]))
    F = eigen(S,sortby=x->1/abs(x))
    U = F.vectors
    d = F.values
    Λ = Array(S^(-1/2))*U#*diagm(abs.(d))^(-1/2)
    idxs = [abs(d[i]) > 1E-7 for i=1:size(S,1)]
    @output "Found {} linear dependencies. Projected them out.\n" size(S,1) - sum(idxs)
    Λ = convert(Array{Float64},real.(Λ[:,idxs]))
    H = Hermitian(aoint["T"] + aoint["V"])
    ndocc = molecule.Nα#size(S,1)
    nvir = size(S,1) - ndocc
    F = Array{Float64,2}(undef, ndocc+nvir, ndocc+nvir)
    Hmax = maximum(diag(H))
    Hmin = minimum(diag(H))
    for i = 1:ndocc+nvir
        F[i,i] = H[i,i]
        for j = i+1:ndocc+nvir
            F[i,j] = 0.875*S[i,j]*(H[i,i] + H[j,j])
            F[j,i] = F[i,j]
        end
    end
    F = F#[idxs,idxs]
    Ft = Λ'*F*Λ

    # Get orbital energies and transformed coefficients
    eps,Ct = eigen(Hermitian(Ft))

    # Reverse transformation to get MO coefficients
    C = Λ*Ct
    Co = C[:,1:ndocc]
    D = Fermi.contract(Co,Co,"um","vm")
    Eguess = RHFEnergy(D,Array(H),F)
    @output "Guess energy: {}\n" Eguess

    RHF(molecule, aoint, C, Λ, Alg)
end

function RHF(molecule::Molecule, aoint::IntegralHelper, Alg::B, guess::CoreGuess) where B <: RHFAlgorithm

    #Form core guess
    @output "Using Core Guess\n"
    S = Hermitian(aoint["S"])
    Λ = S^(-1/2)
    H = Hermitian(aoint["T"] + aoint["V"])
    F = Array{Float64,2}(undef, ndocc+nvir, ndocc+nvir)
    F .= H
    Ft = Λ*F*transpose(Λ)

    # Get orbital energies and transformed coefficients
    eps,Ct = eigen(Hermitian(Ft))

    # Reverse transformation to get MO coefficients
    C = Λ*Ct

    RHF(molecule, aoint, C, Λ, Alg)
end

function RHF(wfn::RHF)
    Alg = select_alg(Fermi.CurrentOptions["scf_alg"])
    ints = Fermi.Integrals.IntegralHelper()
    RHF(wfn, ints, Alg)
end

function RHF(wfn::RHF, aoint::IntegralHelper, Alg::B) where B <: RHFAlgorithm

    # Projection of A→B done using equations described in Werner 2004
    # https://doi.org/10.1080/0026897042000274801
    @output "Using {} wave function as initial guess\n" wfn.basis
    Ca = wfn.ints.orbs["FU"]
    Sbb = aoint["S"]
    S = Hermitian(aoint["S"])
    Λ = S^(-1/2)
    Sab = Lints.projector(wfn.LintsBasis, aoint.LintsBasis)
    T = transpose(Ca)*Sab*(Sbb^-1)*transpose(Sab)*Ca
    Cb = (Sbb^-1)*transpose(Sab)*Ca*T^(-1/2)
    Cb = real.(Cb)
    RHF(Fermi.Geometry.Molecule(), aoint, Cb, Alg)
end

function RHFEnergy(D::Array{Float64,2}, H::Array{Float64,2},F::Array{Float64,2})
    return sum(D .* (H .+ F))
end

#actual HF routine is in here
include("SCF.jl")


text:
function TriMesh_Search( x, eConn, eAdjacency, point )
#=
   TriMesh_Search - Locates an element that contains a given point.
   The function provides only one such element and makes an
   explicit assumption that the triangulation is Delauney.

   Usage:
     e, r, s = TriMesh_Search(x,eConn,eAdjacency,point)

   where:
     x    is an array of the nodal coordinates.
          Dimension: [ nNodes, 2 ]

     eConn  is the element connectivity.  The standard orientation is assumed
            so the first three nodes describe the vertices.
            Dimension: [ nElements, nLocalNoes ]

     eAdjacency  an [nElements, 3] array containing adjacent elements or -1
                 where -1 indicates a boundary (no adjacent element).

     point  coordinates of a point of interest in the mesh.

     e      the index of the element where "point" is located.

     r,s    the finite element coordinate of "point" in the reference element.

   Licensing:
     This code is distributed under the MIT license.

   Author: Jeff Borggaard, Virginia Tech
           part of FEMfunctions.jl
=#

  nElements = size(eConn,1)

  element = 1   # initial guess of the element containing the point
  elementOld = 1
  count = 0
  found = false
  global iso1 = 0.0
  global iso2 = 0.0
  global iso3 = 0.0
  while ~found
    xLocal = x[ eConn[element,1:3], : ]    # view?

#    plot( [ xLocal[1,1] xLocal[2,1] ],[ xLocal[1,2] xLocal[2,2]],...
#          [ xLocal[2,1] xLocal[3,1] ],[ xLocal[2,2] xLocal[3,2]],...
#          [ xLocal[3,1] xLocal[1,1] ],[ xLocal[3,2] xLocal[1,2]] )

    delta = ( xLocal[3,1]-xLocal[2,1] )*( xLocal[1,2]-xLocal[2,2] ) -
            ( xLocal[3,2]-xLocal[2,2] )*( xLocal[1,1]-xLocal[2,1] )

    global iso1  = ( ( xLocal[2,1]-point[1] )*( xLocal[3,2]-point[2] ) -
              ( xLocal[3,1]-point[1] )*( xLocal[2,2]-point[2] ) ) / delta;

    global iso2  = ( ( xLocal[3,1]-point[1] )*( xLocal[1,2]-point[2] ) -
              ( xLocal[1,1]-point[1] )*( xLocal[3,2]-point[2] ) ) / delta;

    global iso3  = ( ( xLocal[1,1]-point[1] )*( xLocal[2,2]-point[2] ) -
              ( xLocal[2,1]-point[1] )*( xLocal[1,2]-point[2] ) ) / delta;

    if iso1>=0 && iso2>=0 && iso3>=0
      found = true

    else
      count = count + 1;
      index = sortperm([iso1, iso2, iso3])

      if eAdjacency[element,index[1]]>0
        element_old = element;
        element = eAdjacency[element,index[1]]
        choice = 1
      elseif eAdjacency[element,index[2]]>0
        element_old = element
        element = eAdjacency[element,index[2]]
        choice = 2
      elseif eAdjacency[element,index[3]]>0
        element_old = element
        element = eAdjacency[element,index[3]]
        choice = 3;
      else
         # iso1
         # iso2
         # iso3
         # eAdjacency[element,:]
         # index
         error("no next element found, Delaunay mesh?")
      end
#       if ( index(choice)==1 )
#         plot( (xLocal(2,1)+xLocal(3,1))/2, (xLocal(2,2)+xLocal(3,2))/2, '*')
#       elseif ( index(choice)==2 )
#         plot( (xLocal(1,1)+xLocal(3,1))/2, (xLocal(1,2)+xLocal(3,2))/2, '*')
#       else
#         plot( (xLocal(1,1)+xLocal(2,1))/2, (xLocal(1,2)+xLocal(2,2))/2, '*')
#       end
       if count>nElements   # if we are somehow in a loop, get out...
         element = nElements
         count   = 0
       end
    end

  end

  return element, iso2, iso3

end


text:
function cosine_similarity(x, y)
    return dot(x, y) / (norm(x) * norm(y))
end

function train_and_find_support_vectors(model, data, solver)

    SVDD.set_data!(model, data)
    try
        SVDD.fit!(model, solver)
    catch e
        # For small data sets fall back to returning all points as SVs
        if size(data, 2) <= 5
            return collect(1:size(data, 2))
        end
        throw(SamplingException("Failed fitting model on data set of size $(size(data, 2)) due to $e."))
    end
    return SVDD.get_support_vectors(model)
end

function mcc(classification, ground_truth; pos_class = :outlier, neg_class = :inlier)::Float64
    @assert length(classification) == length(ground_truth)
    tp = sum((classification .== pos_class) .& (ground_truth .== pos_class))
    fp = sum((classification .== pos_class) .& (ground_truth .== neg_class))
    tn = sum((classification .== neg_class) .& (ground_truth .== neg_class))
    fn = sum((classification .== neg_class) .& (ground_truth .== pos_class))

    d = sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn))
    return d == 0 ? 0.0 : (tp * tn - fp * fn) / d
end


text:
# https://www.ibm.com/support/knowledgecenter/SSSA5P_12.10.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/EpOpt.html
const CPLEX_REDUCED_COST_TOLERANCE = 1.0e-6

function _check_weight_matrix(weight_matrix)
    # Return true if there may be a new path in this weight matrix, false if it is hopeless.

    # If all weights are zero, impossible to get a negative reduced cost.
    # Very efficient test, unlike the next one.
    if nnz(weight_matrix) == 0
        return false
    end

    # If all weights are negative, there is a negative-cost cycle.
    return ! all(weight_matrix .<= CPLEX_REDUCED_COST_TOLERANCE)
end

function _build_path(state::LightGraphs.AbstractPathState, demand::Edge)
    # TODO: use LightGraphs.enumerate_paths? It returns lists of vertices and not edges...
    reverse_nodes = [dst(demand)]
    current_node = dst(demand)
    while current_node != src(demand)
        current_node = state.parents[current_node]
        push!(reverse_nodes, current_node)
    end
    return [Edge(reverse_nodes[i], reverse_nodes[i - 1]) for i in length(reverse_nodes):-1:2]
end

function _oblivious_iterative_solve_pricing_problem_master(rd::RoutingData,
                                                           dual_values_matrices,
                                                           dual_value_convexity,
                                                           matrices::Vector{Dict{Edge{Int}, Float64}})
    paths = Dict{Edge{Int}, Vector{Edge}}()

    for demand in demands(rd)
        # Determine the weight matrix to use for the computations:
        #     for e: \sum_{D\in\Delta} dual_{e, d_k} d_k
        rd.logmessage("Pricing for $demand")
        weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
        for e in edges(rd)
            for matrix in matrices
                if e in keys(dual_values_matrices[matrix])
                    # Ignore positive dual values, they are mostly numerical errors.
                    if dual_values_matrices[matrix][e] >= - CPLEX_REDUCED_COST_TOLERANCE
                        # rd.logmessage("[$matrix][$e] $(dual_values_matrices[matrix][e])")
                        continue
                    end

                    weight_matrix[src(e), dst(e)] -= matrix[demand] * dual_values_matrices[matrix][e]
                end
            end
        end

        # Check that at least some weights are nonzero.
        if ! _check_weight_matrix(weight_matrix)
            continue
        end

        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
            continue
        end
        reduced_cost = state.dists[dst(demand)] - dual_value_convexity[demand]
        if reduced_cost >= - CPLEX_REDUCED_COST_TOLERANCE
            continue
        end
        rd.logmessage(reduced_cost)

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
            rd.logmessage("New path $(path) is already in the formulation!")
        else
            paths[demand] = path
        end
    end

    return paths
end

function _oblivious_iterative_solve_pricing_problem_subproblem(rd::RoutingData, dual_values_capacity)
    paths = Dict{Edge{Int}, Vector{Edge}}()

    # Determine the weight matrix to use for the computations.
    weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
    for e in keys(dual_values_capacity)
        # Ignore negative dual values, they are mostly numerical errors.
        if dual_values_capacity[e] <= CPLEX_REDUCED_COST_TOLERANCE
            # rd.logmessage("[$e] Ignoring dual value: $(dual_values_capacity[e])")
            continue
        end

        weight_matrix[src(e), dst(e)] -= dual_values_capacity[e]
    end

    if ! _check_weight_matrix(weight_matrix)
        return paths
    end

    # Perform the actual pricing.
    for demand in demands(rd)
        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
            continue
        end
        # TODO: state.dists[dst(demand)] - dual_value_convexity[demand] to avoid cycling!?
        if state.dists[dst(demand)] >= - CPLEX_REDUCED_COST_TOLERANCE
            continue
        end
        rd.logmessage(state.dists[dst(demand)])

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
            rd.logmessage("New path $(path) is already in the formulation!")
        else
            paths[demand] = path
        end
    end

    return paths
end

function _add_path_to_master_formulation(rd::RoutingData, rm::RoutingModel, demand::Edge, path_id::Int)
    # Update RoutingModel. Create a new variable for the new path, add it in the corresponding constraints.

    # New variable and convexity constraint.
    add_routing_var!(rm, demand, path_id)

    # Traffic constraints.
    for matrix in keys(rm.constraints_matrices)
        if iszero(matrix[demand])
            continue
        end

        for e in keys(rm.constraints_matrices[matrix])
            if e in rd.paths_edges[path_id]
                set_normalized_coefficient(rm.constraints_matrices[matrix][e], rm.routing[demand, path_id], matrix[demand])
            end
        end
    end
end

function _add_path_to_subproblem_formulation(rd::RoutingData, s_rm::RoutingModel, demand::Edge, path_id::Int)
    # Update RoutingModel. Create a new variable for the new path, add it in the corresponding constraints.

    add_routing_var!(s_rm, demand, path_id)

    # New variable.
    s_rm.routing[demand, path_id] = @variable(s_rm.model, lower_bound=0, upper_bound=1)
    set_name(s_rm.routing[demand, path_id], "routing[$demand, $path_id]")

    # Capacity constraints.
    for e in edges(rd)
        if e in rd.paths_edges[path_id]
            set_normalized_coefficient(s_rm.constraints_capacity[e], s_rm.routing[demand, path_id], 1)
        end
    end
end

function solve_master_problem(rd::RoutingData, rm::RoutingModel, ::Load, ::MinimumMaximum, formulation::FormulationType, ::Val{true}, ::CuttingPlane, ::ObliviousUncertainty, ::UncertainDemand)
    n_new_paths = 0
    @assert rm.mu !== nothing

    result = nothing
    current_routing = nothing

    while true
        # Solve the current master problem.
        m = rm.model
        optimize!(m)
        result = termination_status(m)
        current_routing = Routing(rd, value.(rm.routing)) # TODO: remember all generated routings.

        # Check if there are still columns to add.
        dual_values_matrices = Dict(tm => Dict(e => dual(constraint) for (e, constraint) in d) for (tm, d) in rm.constraints_matrices)
        dual_value_convexity = dual.(rm.constraints_convexity)
        matrices = collect(keys(rm.constraints_matrices))

        new_paths = _oblivious_iterative_solve_pricing_problem_master(rd, dual_values_matrices, dual_value_convexity, matrices)
        n_new_paths += length(new_paths)

        if length(new_paths) == 0
            break
        end

        # Add the new paths to the formulation.
        for (d, path) in new_paths
            path_id = add_path!(rd, d, path)
            _add_path_to_master_formulation(rd, rm, d, path_id)
        end

        # TODO: propose to output problems at each iteration (_export_lp_if_allowed).
    end

    return result, current_routing, n_new_paths
end

function solve_subproblem(rd::RoutingData, rm::RoutingModel, s_rm::RoutingModel, e_bar::Edge, ::Load, ::MinimumMaximum, ::FormulationType, ::Val{true}, ::CuttingPlane, ::ObliviousUncertainty, ::UncertainDemand)
    n_new_paths = 0
    result = nothing

    while true
        # Solve the current subproblem.
        optimize!(s_rm.model)
        result = termination_status(s_rm.model)

        # Check if there are still columns to add.
        dual_values_capacity = Dict(e => dual(constraint) for (e, constraint) in s_rm.constraints_capacity)

        rd.logmessage("Subpricing for edge $e_bar")
        new_paths = _oblivious_iterative_solve_pricing_problem_subproblem(rd, dual_values_capacity)
        n_new_paths += length(new_paths)

        if length(new_paths) == 0
            break
        end

        # Add the new paths to the formulation.
        for (d, path) in new_paths
            path_id = _add_path_to_data(rd, d, path)
            _add_path_to_master_formulation(rd, rm, d, path_id)
            _add_path_to_subproblem_formulation(rd, s_rm, d, path_id)
        end

        # TODO: propose to output problems at each iteration (_export_lp_if_allowed).
    end

    routing_v = value.(s_rm.routing)
    demand = Dict(d => sum(routing_v[d, p_id] for p_id in rd.demand_to_path_ids[d]) for d in demands(rd))
    candidate_matrix = ConsideredTrafficMatrix(objective_value(s_rm.model), e_bar, demand)

    return result, n_new_paths, candidate_matrix
end


text:
module TestAqua

using Aqua
using Maybe
using Test

Aqua.test_all(
    Maybe;
    project_extras = true,
    stale_deps = (; ignore = [:Compat]),
    deps_compat = true,
    project_toml_formatting = true,
)

@testset "Compare test/Project.toml and test/environments/main/Project.toml" begin
    @test Text(read(joinpath(@__DIR__, "Project.toml"), String)) ==
          Text(read(joinpath(@__DIR__, "environments", "main", "Project.toml"), String))
end

end  # module