text:
using Distributions, Random, StaticArrays, LinearAlgebra
using Test

import Distributions: GenericMvTDist
import PDMats: PDMat

# Set location vector mu and scale matrix Sigma as in
# Hofert M. On Sampling from the Multivariate t Distribution. The R Journal
mu = [1., 2]
Sigma = [4. 2; 2 3]

# LogPDF evaluation for varying degrees of freedom df
# Julia's output is compared to R's corresponding values obtained via R's mvtnorm package
# R code exemplifying how the R values (rvalues) were obtained:
# options(digits=20)
# library("mvtnorm")
# mu <- 1:2
# Sigma <- matrix(c(4, 2, 2, 3), ncol=2)
# dmvt(c(-2., 3.), delta=mu, sigma=Sigma, df=1)
rvalues = [-5.6561739738159975133,
           -5.4874952805811396672,
           -5.4441948098568158088,
           -5.432461875138580254,
           -5.4585441614404803801]
df = [1., 2, 3, 5, 10]
for i = 1:length(df)
    d = MvTDist(df[i], mu, Sigma)
    @test isapprox(logpdf(d, [-2., 3]), rvalues[i], atol=1.0e-8)
    dd = typeof(d)(params(d)...)
    @test d.df == dd.df
    @test Vector(d.μ) == Vector(dd.μ)
    @test Matrix(d.Σ) == Matrix(dd.Σ)
end

# test constructors for mixed inputs:
@test typeof(GenericMvTDist(1, Vector{Float32}(mu), PDMat(Sigma))) == typeof(GenericMvTDist(1., mu, PDMat(Sigma)))

@test typeof(GenericMvTDist(1, mu, PDMat(Array{Float32}(Sigma)))) == typeof(GenericMvTDist(1., mu, PDMat(Sigma)))

d = GenericMvTDist(1, Array{Float32}(mu), PDMat(Array{Float32}(Sigma)))
@test typeof(convert(GenericMvTDist{Float64}, d)) == typeof(GenericMvTDist(1, mu, PDMat(Sigma)))
@test typeof(convert(GenericMvTDist{Float64}, d.df, d.dim, d.zeromean, d.μ, d.Σ)) == typeof(GenericMvTDist(1, mu, PDMat(Sigma)))
@test partype(d) == Float32
@test d == deepcopy(d)

@test size(rand(MvTDist(1., mu, Sigma))) == (2,)
@test size(rand(MvTDist(1., mu, Sigma), 10)) == (2,10)
@test size(rand(MersenneTwister(123), MvTDist(1., mu, Sigma))) == (2,)
@test size(rand(MersenneTwister(123), MvTDist(1., mu, Sigma), 10)) == (2,10)

# static array for mean/variance
mu_static = @SVector [1., 2]
# depends on PDMats#101 (merged but not released)
# Sigma_static = @SMatrix [4. 2; 2 3]

for i in 1:length(df)
    d = GenericMvTDist(df[i], mu_static, PDMat(Sigma))
    @test d.μ isa SVector
    @test isapprox(logpdf(d, [-2., 3]), rvalues[i], atol=1.0e-8)
    dd = typeof(d)(params(d)...)
    @test d.df == dd.df
    @test d.μ == dd.μ
    @test Matrix(d.Σ) == Matrix(dd.Σ)
end


text:
using StructJuMP

NS = 3;                        # number of scenarios
probability = [1/3, 1/3, 1/3]; # probability

CROPS  = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
PURCH  = 1:2 # set of crops to purchase (wheat and corn, resp.)
SELL   = 1:4 # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)

Cost     = [150., 230., 260.]    # cost of planting crops
Budget   = 500.                # budget capacity
Purchase = [238., 210.];        # purchase price
Sell     = [170., 150., 36., 10.] # selling price
Yield    = [3.0 3.6 24.0;
            2.5 3.0 20.0;
            2.0 2.4 16.0]
Minreq   = [200., 240., 0.]      # minimum crop requirement

function farmer_stochastic()
    m = StructuredModel(num_scenarios = NS)

    @variable(m, x[i=CROPS] >= 0, Int)
    @objective(m, Min, sum(Cost[i] * x[i] for i=CROPS))
    @constraint(m, sum(x[i] for i=CROPS) <= Budget)

    for s in 1:NS
        blk = StructuredModel(parent = m, id = s, prob = probability[s])
        @variable(blk, y[j=PURCH] >= 0)
        @variable(blk, w[k=SELL] >= 0)
        @objective(blk, Min,
            sum(Purchase[j] * y[j] for j=PURCH)
            - sum(Sell[k] * w[k] for k=SELL))
        @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
        @constraint(blk, const_minreq_beets, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
        @constraint(blk, const_aux, w[3] <= 6000)
    end

    return m
end

function farmer_blocks()
    m = StructuredModel(num_scenarios = NS)

    @variable(m, x[i=CROPS,s=1:NS] >= 0, Int)
    @variable(m, y[j=PURCH,s=1:NS] >= 0)
    @variable(m, w[k=SELL,s=1:NS] >= 0)
    @objective(m, Min,
        sum(probability[s] * Cost[i] * x[i,s] for i=CROPS for s=1:NS)
        + sum(probability[s] * Purchase[j] * y[j,s] for j=PURCH for s=1:NS)
        - sum(probability[s] * Sell[k] * w[k,s] for k=SELL for s=1:NS))
    @constraint(m, nonant[i=CROPS,s=2:NS], x[i,s-1] - x[i,s] == 0)

    for s in 1:NS
        blk = StructuredModel(parent = m, id = s)
        @objective(blk, Min, 0.)
        @constraint(blk, sum(x[i,s] for i=CROPS) <= Budget)
        @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j,s] + y[j,s] - w[j,s] >= Minreq[j])
        @constraint(blk, const_minreq_beets, Yield[s,3] * x[3,s] - w[3,s] - w[4,s] >= Minreq[3])
        @constraint(blk, const_aux, w[3,s] <= 6000)
    end
    return m
end

ext: jl


text:
"Type-safe wrapper for `Int64`, references a variable."
struct VarRef
    val::Int64
end

"Type-safe wrapper for `Int64`, references a constraint."
struct ConsRef
    val::Int64
end

#to be moved to MOI_wrapper
"""
SCIPData holds pointers to SCIP data.
It does not perform memory management and should not be created directly.
"""
mutable struct SCIPData
    scip::Ref{Ptr{SCIP_}}
    vars::Dict{VarRef, Ref{Ptr{SCIP_VAR}}}
    conss::Dict{ConsRef, Ref{Ptr{SCIP_CONS}}}
    var_count::Int64
    cons_count::Int64

    # Map from user-defined Julia types (keys are <: AbstractConstraintHandler
    # or <: AbstractConstraint, respectively) to the corresponding SCIP objects.
    # The reverse mapping is handled by SCIP itself.
    # This also serves to prevent premature GC.
    conshdlrs::Dict{Any, Ptr{SCIP_CONSHDLR}}
    conshdlrconss::Dict{Any, Ptr{SCIP_CONS}}

    # Map from user-defined types (keys are <: AbstractSeparator) to the
    # corresponding SCIP objects.
    sepas::Dict{Any, Ptr{SCIP_SEPA}}
end

# Protect SCIPData from GC for ccall with Ptr{SCIP_} argument.
Base.unsafe_convert(::Type{Ptr{SCIP_}}, scipd::SCIPData) = scipd.scip[]

"Release references and free memory."
function free_scip(scipd::SCIPData)
    # Avoid double-free (SCIP will set the pointers to NULL).
    if scipd.scip[] != C_NULL
        for c in values(scipd.conss)
            @SCIP_CALL SCIPreleaseCons(scipd, c)
        end
        for v in values(scipd.vars)
            @SCIP_CALL SCIPreleaseVar(scipd, v)
        end
        # only scipd.scip is GC-protected during ccall!
        GC.@preserve scipd begin
            @SCIP_CALL SCIPfree(scipd.scip)
        end
    end
    @assert scipd.scip[] == C_NULL
end

"Set a parameter's current value."
function get_parameter(scipd::SCIPData, name::AbstractString)
    param = SCIPgetParam(scipd, name)
    if param == C_NULL
        error("Unrecognized parameter: $name")
    end
    paramtype = SCIPparamGetType(param)
    if paramtype === SCIP_PARAMTYPE_BOOL
        value = Ref{SCIP_Bool}()
        @SCIP_CALL SCIPgetBoolParam(scipd, name, value)
        if value[] == TRUE
            return true
        elseif value[] == FALSE
            return false
        else
            error("encountered invalid value for a boolean: $(value[])")
        end
    elseif paramtype === SCIP_PARAMTYPE_INT
        value = Ref{Cint}()
        @SCIP_CALL SCIPgetIntParam(scipd, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_LONGINT
        value = Ref{Clonglong}()
        @SCIP_CALL SCIPgetLongintParam(scipd, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_REAL
        value = Ref{Cdouble}()
        @SCIP_CALL SCIPgetRealParam(scipd, name, value)
        return value[]
    elseif paramtype === SCIP_PARAMTYPE_CHAR
        value = Ref{Cchar}()
        @SCIP_CALL SCIPgetCharParam(scipd, name, value)
        return Char(value[])
    elseif paramtype === SCIP_PARAMTYPE_STRING
        value = Ref{Ptr{Cchar}}()
        @SCIP_CALL SCIPgetStringParam(scipd, name, value)
        return unsafe_string(value[])
    else
        error("Unexpected parameter type: $paramtype")
    end
end

"Set a parameter."
function set_parameter(scipd::SCIPData, name::AbstractString, value)
    param = SCIPgetParam(scipd, name)
    if param == C_NULL
        error("Unrecognized parameter: $name")
    end
    paramtype = SCIPparamGetType(param)
    if paramtype === SCIP_PARAMTYPE_BOOL
        @SCIP_CALL SCIPsetBoolParam(scipd, name, value)
    elseif paramtype === SCIP_PARAMTYPE_INT
        @SCIP_CALL SCIPsetIntParam(scipd, name, value)
    elseif paramtype === SCIP_PARAMTYPE_LONGINT
        @SCIP_CALL SCIPsetLongintParam(scipd, name, value)
    elseif paramtype === SCIP_PARAMTYPE_REAL
        @SCIP_CALL SCIPsetRealParam(scipd, name, value)
    elseif paramtype === SCIP_PARAMTYPE_CHAR
        @SCIP_CALL SCIPsetCharParam(scipd, name, value)
    elseif paramtype === SCIP_PARAMTYPE_STRING
        @SCIP_CALL SCIPsetStringParam(scipd, name, value)
    else
        error("Unexpected parameter type: $paramtype")
    end
    return nothing
end

"Return pointer to SCIP variable."
function var(scipd::SCIPData, vr::VarRef)::Ptr{SCIP_VAR}
    return scipd.vars[vr][]
end

"Return pointer to SCIP constraint."
function cons(scipd::SCIPData, cr::ConsRef)::Ptr{SCIP_CONS}
    return scipd.conss[cr][]
end

"Store reference to variable, return VarRef"
function store_var!(scipd::SCIPData, var__::Ref{Ptr{SCIP_VAR}})
    scipd.var_count += 1
    vr = VarRef(scipd.var_count)
    scipd.vars[vr] = var__
    return vr
end

"Store reference to constraint, return ConsRef"
function store_cons!(scipd::SCIPData, cons__::Ref{Ptr{SCIP_CONS}})
    scipd.cons_count += 1
    cr = ConsRef(scipd.cons_count)
    scipd.conss[cr] = cons__
    return cr
end

"Add variable to problem (continuous, no bounds), return var ref."
function add_variable(scipd::SCIPData)
    var__ = Ref{Ptr{SCIP_VAR}}(C_NULL)
    @SCIP_CALL SCIPcreateVarBasic(scipd, var__, "", -SCIPinfinity(scipd), SCIPinfinity(scipd),
                           0.0, SCIP_VARTYPE_CONTINUOUS)
    @SCIP_CALL SCIPaddVar(scipd, var__[])
    return store_var!(scipd, var__)
end

"Delete variable from problem."
function delete(scipd::SCIPData, vr::VarRef)
    # delete variable from SCIP problem
    deleted = Ref{SCIP_Bool}()
    @SCIP_CALL SCIPdelVar(scipd, var(scipd, vr), deleted)
    deleted[] == TRUE || error("Variable at $(vr.val) could not be deleted!")

    # release memory and remove reference
    @SCIP_CALL SCIPreleaseVar(scipd, scipd.vars[vr])
    delete!(scipd.vars, vr)
    return nothing
end

"Delete constraint from problem."
function delete(scipd::SCIPData, cr::ConsRef)
    @SCIP_CALL SCIPdelCons(scipd, cons(scipd, cr))

    # release memory and remove reference
    @SCIP_CALL SCIPreleaseCons(scipd, scipd.conss[cr])
    delete!(scipd.conss, cr)
    return nothing
end

"""
Add (ranged) linear constraint to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references for affine terms.
- `coefs::AbstractArray{Float64}`: coefficients for affine terms.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_linear_constraint(scipd::SCIPData, varrefs, coefs, lhs, rhs)
    @assert length(varrefs) == length(coefs)
    vars = [var(scipd, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicLinear(
        scipd, cons__, "", length(vars), vars, coefs, lhs, rhs)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add (ranged) quadratic constraint to problem, return cons ref.

# Arguments
- `linrefs::AbstractArray{VarRef}`: variable references for affine terms.
- `lincoefs::AbstractArray{Float64}`: coefficients for affine terms.
- `quadrefs1::AbstractArray{VarRef}`: first variable references for quadratic terms.
- `quadrefs2::AbstractArray{VarRef}`: second variable references for quadratic terms.
- `quadcoefs::AbstractArray{Float64}`: coefficients for quadratic terms.
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_quadratic_constraint(scipd::SCIPData, linrefs, lincoefs,
                                  quadrefs1, quadrefs2, quadcoefs, lhs, rhs)
    @assert length(linrefs) == length(lincoefs)
    @assert length(quadrefs1) == length(quadrefs2)
    @assert length(quadrefs1) == length(quadcoefs)

    linvars = [var(scipd, vr) for vr in linrefs]
    quadvars1 = [var(scipd, vr) for vr in quadrefs1]
    quadvars2 = [var(scipd, vr) for vr in quadrefs2]

    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicQuadratic(
        scipd, cons__, "", length(linvars), linvars, lincoefs,
        length(quadvars1), quadvars1, quadvars2, quadcoefs, lhs, rhs)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add second-order-cone constraint to problem, return cons ref.

Does not support the full generality of SCIP's constraint (offsets and
coefficients). The first entry in `varrefs` is used for the special variable on
the right-hand side.
"""
function add_second_order_cone_constraint(scipd::SCIPData, varrefs)
    vars = [var(scipd, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicSOC(scipd, cons__, "", length(vars) - 1,
                               vars[2:end], C_NULL, C_NULL, 0.0, vars[1], 1.0, 0.0)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add special-ordered-set of type 1 to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references
- `weights::AbstractArray{Float64}`: numeric weights
"""
function add_special_ordered_set_type1(scipd::SCIPData, varrefs, weights)
    @assert length(varrefs) == length(weights)
    vars = [var(scipd, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicSOS1(scipd, cons__, "", length(vars), vars, weights)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add special-ordered-set of type 2 to problem, return cons ref.

# Arguments
- `varrefs::AbstractArray{VarRef}`: variable references
- `weights::AbstractArray{Float64}`: numeric weights
"""
function add_special_ordered_set_type2(scipd::SCIPData, varrefs, weights)
    @assert length(varrefs) == length(weights)
    vars = [var(scipd, vr) for vr in varrefs]
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicSOS2(scipd, cons__, "", length(vars), vars, weights)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add abspower constraint to problem, return cons ref.

    lhs ≤ sign(x + a) * abs(x + a)^n + c*z ≤ rhs

# Arguments
- `x::VarRef`: reference for power variable
- `a::Float64`: offset for power variable
- `n::Float64`: exponent for power variable, n >= 1
- `z::VarRef`: reference for linear variable
- `c::Float64`: coefficient for linear variable
- `lhs::Float64`: left-hand side for ranged constraint
- `rhs::Float64`: right-hand side for ranged constraint

Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
"""
function add_abspower_constraint(scipd::SCIPData, x, a, n, z, c, lhs, rhs)
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SCIP_CALL SCIPcreateConsBasicAbspower(
        scipd, cons__, "", var(scipd, x), var(scipd, z), n, a, c, lhs, rhs)
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

"""
Add indicator constraint to problem, return cons ref.

    y = 1 ==> a^T x ≤ rhs

y has to be a binary variable, or SCIP will error.

# Arguments
- `y::VarRef`: reference for binary indicator variable
- `x::Vector{VarRef}`: reference vector for variables
- `a::Float64`: coefficients for x variable
- `rhs::Float64`: right-hand side for linear constraint
"""
function add_indicator_constraint(scipd::SCIPData, y, x, a, rhs)
    SCIPvarIsBinary(var(scipd, y)) > 0 || error("indicator variable must be binary.")
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    xref = [var(scipd, x[i]) for i in eachindex(x)]
    @SCIP_CALL SCIPcreateConsBasicIndicator(
        scipd, cons__, "", var(scipd, y), length(x), xref, a, rhs
    )
    @SCIP_CALL SCIPaddCons(scipd, cons__[])
    return store_cons!(scipd, cons__)
end

# Transform SCIP C function name as follows:
# 1. Remove leading "SCIP" part (drop the first four characters).
# 2. Convert camel case to snake case.
# For example, `SCIPprintStatusStatistics` becomes `print_status_statistics`.
const STATISTICS_FUNCS = map(x -> Symbol(camel_case_to_snake_case(string(x)[5 : end])), SCIP_STATISTICS_FUNCS)

for (scip_statistics_func, statistics_func) in zip(SCIP_STATISTICS_FUNCS, STATISTICS_FUNCS)
    @eval begin
        """
            $($statistics_func)(scipd::SCIPData)

        Print statistics (calls `$($scip_statistics_func)`) to standard output.
        """
        function $statistics_func end

        function $statistics_func(scipd::SCIPData)
            ret = $scip_statistics_func(scipd, C_NULL)
            ret !== nothing && @assert ret == SCIP_OKAY
            return nothing
        end
    end
end


text:
# asymptotic Gaussian likelihood function of statistic
function logL(θ, m, n, η, ϵ, withdet=true)
    S = size(η,2)
    k = size(m,1)
    ms = zeros(S, k)
    Threads.@threads for s = 1:S
        y, junk = SVmodel(θ, n, η[:,s], ϵ[:,s])
        y = min.(y, 100.0)
        y = max.(y,-100.0)
        ms[s,:] = sqrt(n)*aux_stat(y)
    end
    mbar = mean(ms,dims=1)[:]
    if ~any(isnan.(mbar))
        Σ = cov(ms)
        x = (m .- mbar)
        lnL = try
            if withdet
                lnL = -0.5*log(det(Σ)) - 0.5*x'*inv(Σ)*x # for Bayesian
            else
                lnL = 0.5*x'*inv(Σ)*x # for classic indirect inference (note sign change)
            end
        catch
            lnL = -Inf
        end
     else
         lnL = -Inf
     end
     return lnL
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
"""
Count the number of individuals in each state at each event time
"""
function state_timeseries(event_db::edb)
	susceptible = fill(0, length(event_db.event_times[event_db.event_times .< Inf]))
	infectious = fill(0, length(event_db.event_times[event_db.event_times .< Inf]))
	if size(event_db.events)[2] == 4
		recovered = fill(0, length(event_db.event_times[event_db.event_times .< Inf]))
		for t = 1:length(event_db.event_times[event_db.event_times .< Inf])
			susceptible[t] = sum(find_state(event_db, event_db.event_times[t], "S"))
			infectious[t] = sum(find_state(event_db, event_db.event_times[t], "I"))
			recovered[t] = sum(find_state(event_db, event_db.event_times[t], "R"))
		end
		return DataFrame(time = event_db.event_times[event_db.event_times .< Inf], S = susceptible, I = infectious, R = recovered)
	end
  if size(event_db.events)[2] == 3
    for t = 1:length(event_db.event_times[event_db.event_times .< Inf])
        susceptible[t] = sum(find_state(event_db, event_db.event_times[t], "S"))
        infectious[t] = sum(find_state(event_db, event_db.event_times[t], "I"))
    end
    return DataFrame(time = event_db.event_times[event_db.event_times .< Inf], S = susceptible, I = infectious)
  end
end

"""
Determine infection status for each individual at a specific time,
and create a dataframe which includes their coordinates for easy
spatial plotting of status through time
"""
function spatial_states(pop_db::DataFrame, event_db::edb, time::Real)
  i = find_state(event_db, time, "I")
  state = fill("S", length(i))
  state[i] = "I"
  if size(event_db.events)[2] == 4
      state[find_state(event_db, time, "R")] = "R"
  end
  return append!(DataFrame(x=fill(NaN, 3), y=fill(NaN, 3), state=["S", "I", "R"]), DataFrame(x=pop_db[:,2], y=pop_db[:,3], state=state))
end

"""
Make a spatial animation/slider in iJulia
"""
function state_animation(framewindow, event_db, pop_db)
  @manipulate for time = 1:framewindow:ceil(maximum(event_db.event_times[event_db.event_times .< Inf]))
      plot(spatial_states(pop_db, event_db, time), x="x", y="y", Geom.point, color="state", Scale.discrete_color_manual("black", "orange", "blue"))
  end
end


text:
# PF00501 in Pfam 30.0 has 3560 columns, 423 sin inserts
@benchgroup "filtercolumns" ["IO", "Pfam", "MSA"] begin

    seq = replace("""
    ........................................................................................
    ........................................................................................
    .m-----.----...--.---..--..........................-------....-..-..-.....-....-....-...
    ......-....-....-....-.....-.......-........-..........-........-..............-....-...
    ............-........G..............V......E.............K.............G.........D.I....
    ...............I.G.......L....K..........G...................R....N........V............
    .P......E..............W.....L.........I..A......D.L.......G..V...Q.........M.....A.....
    ..........G...G......C....S.....L.....N......L.P.Y......................................
    ........Q.Q..K.E..................E...I.....M...V...D.......L.....L....H......E.....I...
    ....G.......T.......-....-....--..........................--..........................-.
    -.--..-..-.......--..-....-.....-.....-.....-......-......-.........-......-........-...
    ....-.....-........-......-........-.......-.......-.......-.........-......-......-....
    ....-....-........-..........-......-........-....-...-..-..-.-.-..-.--.---........---..
    ........................................................................................
    .........................................----.-.--.-......-.-.-.-...-....-.....-...-....
    ..-............-...............-...........-.......-...........-......-.........-.......
    ...-............-...............-.................-.......-.....-.....-...---...........
    ...........--.-.--.-.--.........----.......--..-.--.....................................
    ...--...-.-....--..-....-..-...-...-....-..-....-..-..................-....-....-.-...-.
    ...-....-............-..........-.....-..........-.........-......-............-......-.
    ...-.....-.......-............-.......-...........-......-......-....--..........----...
    ..........................-..-.-..-.........-...........-.........-...........-.........
    -......-.......-.........-......-...-...-..........-...-....-..--..--...................
    ..........--..-.-.---..-.............-----..............................................
    .....................................................................................---
    ...................................-...--..-............................................
    ........................................................................................
    .............................................-.....-...........................-.......-
    -..-...-...-....-....-......-.........-..-...-...-..-..-..-.----........................
    ...................--.-..-..-.-....-..........-.....-.-.-----........--.................
    ...............................................................---.---.-...-...-.-.-....
    .-...-.....-........................--.....-....................................-....-..
    -...-...-..-......-....-..............-.....-.......-........-.......-...-..-.-......-..
    ......................................................--.--..........-.-.............-..
    -..........-....-.....-..-.............-.................-..-............-.............-
    ......-....-.........-...-..................-...--...................-...-..-..-...--..-
    .......................................................-....--..-.-...-...-...-..-..---.
    ...-..-..............-..........-...-................................-...-....-.-..-..-.
    -.-...-..-..-----vvys...................................................................
    ........................................................................................
    ........................................
    """, '\n' => "")

    annot = Annotations()
    setannotresidue!(annot, "K1PKS6_CRAGI/1-58", "SEQ", seq)

    mask = convert(BitArray, Bool[isuppercase(char) || char == '-' for char in seq])
    indexes = collect(1:length(seq))[mask]

    @bench "boolean mask" filtercolumns!(copy($annot), $mask)
    @bench "index array"  filtercolumns!(copy($annot), $indexes)
end


text:
# # Translational Oscillations by a Rotational Actuator (TORA)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Sherlock-Benchmark-9-TORA.ipynb)

module TORA_Sigmoid  #jl

using ClosedLoopReachability, MAT
using ClosedLoopReachability: UniformAdditivePostprocessing, NoSplitter

# The following option determines whether the verification settings should be
# used or not. The verification settings are chosen to show that the safety
# property is satisfied. Concretely we split the initial states into small
# chunks and run many analyses.
const verification = false;

# This model consists of a cart attached to a wall with a spring. The cart is
# free to move on a friction-less surface. The car has a weight attached to an
# arm, which is free to rotate about an axis. This serves as the control input
# to stabilize the cart at $x = 0$.
#
# ## Model
#
# The model is four dimensional, given by the following equations:
#
# ```math
# \left\{ \begin{array}{lcl}
#       \dot{x}_1 &=& x_2 \\
#       \dot{x}_2 &=& -x_1 + 0.1 \sin(x_3) \\
#       \dot{x}_3 &=& x_4  \\
#       \dot{x}_4 &=& u
# \end{array} \right.
# ```
#
# A neural network controller was trained for this system. The trained network
# has 3 hidden layers, with 100 neurons in each layer (i.e., a total of 300
# neurons). Note that the output of the neural network $f(x)$ needs to be
# normalized in order to obtain $u$, namely $u = f(x) - 10$. The sampling time
# for this controller is 1s.

@taylorize function TORA!(dx, x, p, t)
    x₁, x₂, x₃, x₄, u = x

    aux = 0.1 * sin(x₃)
    dx[1] = x₂
    dx[2] = -x₁ + aux
    dx[3] = x₄
    dx[4] = u
    dx[5] = zero(u)
    return dx
end

path = @modelpath("Sherlock-Benchmark-9-TORA", "controllerToraSigmoid.mat")
controller = read_nnet_mat(path, act_key="act_fcns");

# ## Specification

# The verification problem is safety. For an initial set of $x_1 ∈ [0.6, 0.7]$,
# $x_2 ∈ [−0.7, −0.6]$, $x_3 ∈ [−0.4, −0.3]$, and $x_4 ∈ [0.5, 0.6]$, the
# system has to stay within the box $x ∈ [−2, 2]^4$ for a time window of 20s.

X₀ = Hyperrectangle(low=[0.6, -0.7, -0.4, 0.5], high=[0.7, -0.6, -0.3, 0.6])
U = ZeroSet(1)

vars_idx = Dict(:state_vars=>1:4, :control_vars=>5)
ivp = @ivp(x' = TORA!(x), dim: 5, x(0) ∈ X₀ × U)

period = 1.0  # control period
control_postprocessing = UniformAdditivePostprocessing(-10.0)  # control postprocessing

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       postprocessing=control_postprocessing)

## Safety specification: x[1], x[2], x[3], x[4] ∈ [-2, 2] for all t ≤ T
T = 20.0  # time horizon
T_warmup = 2 * period  # shorter time horizon for dry run
T_reach = verification ? T : 5.0  # shorter time horizon if not verifying

safe_states = cartesian_product(BallInf(zeros(4), 2.0), Universe(1))
predicate = X -> X ⊆ safe_states;

# ## Results

alg = TMJets(abstol=1e-10, orderT=8, orderQ=3)
alg_nn = DeepZ()
if verification
    splitter = BoxSplitter([4, 4, 3, 5])
else
    splitter = NoSplitter()
end

function benchmark(; T=T, silent::Bool=false)
    ## We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg,
                                 splitter=splitter)
    sol = res_sol.value
    silent || print_timed(res_sol)

    ## Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    solz = overapproximate(sol, Zonotope)
    res_pred = @timed predicate(solz)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is satisfied.")
    else
        silent || println("The property may be violated.")
    end
    return solz
end

benchmark(T=T_warmup, silent=true)  # warm-up
res = @timed benchmark(T=T_reach)  # benchmark
sol = res.value
println("total analysis time")
print_timed(res);

# We also compute some simulations:

import DifferentialEquations

println("simulation")
res = @timed simulate(prob, T=T; trajectories=10, include_vertices=true)
sim = res.value
print_timed(res);

# Finally we plot the results
using Plots
import DisplayAs

function plot_helper(fig, vars)
    if vars[1] == 0
        safe_states_projected = project(safe_states, [vars[2]])
        time = Interval(0, T)
        safe_states_projected = cartesian_product(time, safe_states_projected)
    else
        safe_states_projected = project(safe_states, vars)
    end
    plot!(fig, safe_states_projected, color=:lightgreen, lab="safe states")
    if !verification && 0 ∉ vars
        plot!(fig, project(X₀, vars), lab="X₀")
    end
    plot!(fig, sol, vars=vars, color=:yellow, lab="")
    plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
    fig = DisplayAs.Text(DisplayAs.PNG(fig))
end

vars = (1, 2)
fig = plot(xlab="x₁", ylab="x₂")
fig = plot_helper(fig, vars)
## savefig("TORA-x1-x2.png")
fig

#-

vars = (1, 3)
fig = plot(xlab="x₁", ylab="x₃")
plot_helper(fig, vars)

#-

vars = (1, 4)
fig = plot(xlab="x₁", ylab="x₄")
plot_helper(fig, vars)

#-

vars=(2, 3)
fig = plot(xlab="x₂", ylab="x₃")
plot_helper(fig, vars)

#-

vars=(2, 4)
fig = plot(xlab="x₂", ylab="x₄")
plot_helper(fig, vars)

#-

vars=(3, 4)
fig = plot(xlab="x₃", ylab="x₄")
fig = plot_helper(fig, vars)
## savefig("TORA-x3-x4.png")
fig

#-

vars = (0, 1)
fig = plot(xlab="t", ylab="x₁")
plot_helper(fig, vars)

#-

vars=(0, 2)
fig = plot(xlab="t", ylab="x₂")
plot_helper(fig, vars)

#-

vars=(0, 3)
fig = plot(xlab="t", ylab="x₃")
plot_helper(fig, vars)

#-

vars=(0, 4)
fig = plot(xlab="t", ylab="x₄")
plot_helper(fig, vars)

#-

# Here we plot the control functions for each run:
tdom = range(0, 20, length=length(controls(sim, 1)))
fig = plot(xlab="t", ylab="u")
[plot!(fig, tdom, [c[1] for c in controls(sim, i)], lab="") for i in 1:length(sim)]
fig = DisplayAs.Text(DisplayAs.PNG(fig))

#-

end  #jl
nothing  #jl


text:
"""
DFTK --- The density-functional toolkit. Provides functionality for experimenting
with plane-wave density-functional theory algorithms.
"""
module DFTK

using Printf
using Markdown
using LinearAlgebra
using Requires
using TimerOutputs
using spglib_jll
using Unitful
using UnitfulAtomic

export Vec3
export Mat3
export mpi_nprocs
export mpi_master
export setup_threading, disable_threading
include("common/timer.jl")
include("common/asserting.jl")
include("common/constants.jl")
include("common/ortho.jl")
include("common/types.jl")
include("common/spherical_harmonics.jl")
include("common/split_evenly.jl")
include("common/mpi.jl")
include("common/threading.jl")
include("common/printing.jl")

export PspHgh
include("pseudo/NormConservingPsp.jl")
include("pseudo/PspHgh.jl")

export ElementPsp
export ElementCohenBergstresser
export ElementCoulomb
export charge_nuclear
export charge_ionic
export n_elec_valence
export n_elec_core
include("elements.jl")

export Smearing
export Model
export PlaneWaveBasis
export compute_fft_size
export G_vectors
export G_vectors_cart
export r_vectors
export r_vectors_cart
export Kpoint
export G_to_r
export G_to_r!
export r_to_G
export r_to_G!
include("Smearing.jl")
include("Model.jl")
include("structure.jl")
include("PlaneWaveBasis.jl")
include("fft.jl")
include("orbitals.jl")
include("show.jl")

export Energies
include("Energies.jl")

export Hamiltonian
export HamiltonianBlock
export energy_hamiltonian
export Kinetic
export ExternalFromFourier
export ExternalFromReal
export AtomicLocal
export PowerNonlinearity
export Hartree
export Xc
export AtomicNonlocal
export Ewald
export PspCorrection
export Entropy
export Magnetic
export Anyonic
export apply_kernel
export compute_kernel
include("terms/terms.jl")

include("occupation.jl")
export compute_density
export total_density
export spin_density
export ρ_from_total_and_spin
include("densities.jl")
include("interpolation_transfer.jl")
export compute_transfer_matrix
export transfer_blochwave
export transfer_blochwave_kpt

export PreconditionerTPA
export PreconditionerNone
export lobpcg_hyper
export diag_full
export diagonalize_all_kblocks
include("eigen/preconditioners.jl")
include("eigen/diag.jl")

export model_atomic
export model_DFT
export model_PBE
export model_LDA
include("standard_models.jl")

export KerkerMixing, KerkerDosMixing, SimpleMixing, DielectricMixing
export LdosMixing, HybridMixing, χ0Mixing
export scf_nlsolve_solver
export scf_damping_solver
export scf_anderson_solver
export scf_CROP_solver
export self_consistent_field
export direct_minimization
export newton
export load_scfres, save_scfres
include("scf/chi0models.jl")
include("scf/mixing.jl")
include("scf/scf_solvers.jl")
include("scf/self_consistent_field.jl")
include("scf/direct_minimization.jl")
include("scf/newton.jl")
include("scf/scfres.jl")
include("scf/potential_mixing.jl")

export symmetry_operations
export standardize_atoms
export bzmesh_uniform
export bzmesh_ir_wedge
export kgrid_from_minimal_spacing, kgrid_from_minimal_n_kpoints
include("symmetry.jl")
include("bzmesh.jl")

export guess_density
export random_density
export load_psp
export list_psp
include("guess_density.jl")
include("pseudo/load_psp.jl")
include("pseudo/list_psp.jl")

export pymatgen_structure
export ase_atoms
export load_lattice
export load_basis
export load_model
export load_density
export load_atoms
export load_magnetic_moments
include("external/abinit.jl")
include("external/load_from_python.jl")
include("external/load_from_file.jl")
include("external/ase.jl")
include("external/pymatgen.jl")

export compute_bands
export high_symmetry_kpath
export plot_bandstructure
include("postprocess/band_structure.jl")

export compute_forces
export compute_forces_cart
include("postprocess/forces.jl")
export compute_stresses
include("postprocess/stresses.jl")
export compute_dos
export compute_ldos
export compute_nos
export plot_dos
include("postprocess/dos.jl")
export compute_χ0
export apply_χ0
include("postprocess/chi0.jl")
export compute_current
include("postprocess/current.jl")

# ForwardDiff workarounds
include("workarounds/dummy_inplace_fft.jl")
include("workarounds/forwarddiff_rules.jl")

function __init__()
    # Use "@require" to only include fft_generic.jl once IntervalArithmetic or
    # DoubleFloats has been loaded (via a "using" or an "import").
    # See https://github.com/JuliaPackaging/Requires.jl for details.
    #
    # The global variable GENERIC_FFT_LOADED makes sure that things are
    # only included once.
    @require IntervalArithmetic="d1acc4aa-44c8-5952-acd4-ba5d80a2a253" begin
        include("workarounds/intervals.jl")
        !isdefined(DFTK, :GENERIC_FFT_LOADED) && include("workarounds/fft_generic.jl")
    end
    @require DoubleFloats="497a8b3b-efae-58df-a0af-a86822472b78" begin
        !isdefined(DFTK, :GENERIC_FFT_LOADED) && include("workarounds/fft_generic.jl")
    end
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plotting.jl")
    @require JLD2="033835bb-8acc-5ee8-8aae-3f567f8a3819"  include("external/jld2io.jl")
    @require WriteVTK="64499a7a-5c06-52f2-abe2-ccb03c286192" include("external/vtkio.jl")
    @require NCDatasets="85f8d34a-cbdd-5861-8df4-14fed0d494ab" begin
        include("external/etsf_nanoquanta.jl")
    end
end

end # module DFTK


text:
using Test
using Base.Iterators: Iterators
using DataStructures: OrderedDict
using QuantumLattices.Essentials.QuantumNumbers
using QuantumLattices.Interfaces: ⊕, ⊗, dimension, expand, permute, decompose
using QuantumLattices.Prerequisites.Traits: contentnames, getcontent

import QuantumLattices.Essentials.QuantumNumbers: periods

@abeliannumber "CN" (:N,) (Inf,)
@abeliannumber "Z4" (:Z,) (4,)
@abeliannumber "CNZ4" (:N, :Z) (Inf, 4)

@testset "regularize" begin
    @test regularize(CNZ4, [1.5, 5.0]) == [1.5, 1.0]
    @test regularize(CNZ4, [4.0, -1.0]) == [4.0, 3.0]
    @test regularize(CNZ4, [1.5 4.0; 5.0 -1.0]) == [1.5 4.0; 1.0 3.0]
end

@testset "AbelianNumbers" begin
    @test contentnames(AbelianNumbers) == (:form, :table, :indptr)

    qn₁, qn₂ = CNZ4(+1.0, +3.0), CNZ4(-1.0, +1.0)
    qns₁ = AbelianNumbers('U', [qn₁, qn₂], [0, 2, 5], :indptr)
    qns₂ = AbelianNumbers('U', [qn₁, qn₂], [2, 3], :counts)
    @test isequal(qns₁, qns₂)
    @test periods(qn₁) == periods(qn₂) == (Inf, 4)
    qns₃ = AbelianNumbers(OrderedDict(qn₁=>2, qn₂=>3))
    qns₄ = AbelianNumbers(OrderedDict(qn₁=>1:2, qn₂=>3:5))
    @test qns₁ == qns₂ == qns₃ == qns₄

    qns = AbelianNumbers('U', [qn₁, qn₂])
    @test qns == AbelianNumbers('U', [qn₁, qn₂], [1, 1], :counts)

    qns = AbelianNumbers('U', [qn₁, qn₂], [0, 3, 5], :indptr)
    @test getcontent(qns, Val(:table)) == qns.contents
    @test qns|>dimension == 5
    @test qns|>string == "QNS(2, 5)"
    @test qns|>size == (2,)
    @test qns|>eltype == CNZ4
    @test qns|>typeof|>eltype == CNZ4
    @test qns|>issorted == false
    @test repr(qns) == "QNS(CNZ4(1.0, 3.0)=>1:3, CNZ4(-1.0, 1.0)=>4:5)"
    @test qns[1] == qn₁
    @test qns[2] == qn₂
    @test qns[1:2] == AbelianNumbers('U', [qn₁, qn₂], [0, 3, 5], :indptr)
    @test qns[[2, 1]] == AbelianNumbers('G', [qn₂, qn₁], [0, 2, 5], :indptr)
    @test qns|>collect == [qn₁, qn₂]
    @test qns|>Iterators.reverse|>collect == [qn₂, qn₁]
    @test qns|>keys|>collect == [qn₁, qn₂]
    @test values(qns, :indptr)|>collect == [1:3, 4:5]
    @test values(qns, :counts)|>collect == [3, 2]
    @test pairs(qns, :indptr)|>collect == [qn₁=>1:3, qn₂=>4:5]
    @test pairs(qns, :counts)|>collect == [qn₁=>3, qn₂=>2]
end

@testset "OrderedDict" begin
    qn₁, qn₂ = CNZ4(1.0, 2.0), CNZ4(2.0, 3.0)
    qns = AbelianNumbers('U', [qn₁, qn₂], [2, 3], :counts)
    @test OrderedDict(qns, :indptr) == OrderedDict(qn₁=>1:2, qn₂=>3:5)
    @test OrderedDict(qns, :counts) == OrderedDict(qn₁=>2, qn₂=>3)
end

@testset "arithmetic" begin
    qn = CNZ4(1.0, 3.0)
    @test qn|>dimension == 1
    @test qn|>typeof|>dimension == 1
    @test +qn == CNZ4(+1.0, +3.0)
    @test -qn == CNZ4(-1.0, +1.0)
    @test qn*4 == 4*qn == CNZ4(4.0, 0.0)
    @test qn^3 == CNZ4(3.0, 1.0)

    qn₁, qn₂ = CNZ4(1.0, 2.0), CNZ4(2.0, 3.0)
    @test qn₁+qn₂ == CNZ4(3.0, 1.0)
    @test qn₁-qn₂ == CNZ4(-1.0, 3.0)
    @test kron(qn₁, qn₂, signs=(+1, +1)) == CNZ4(+3.0, 1.0)
    @test kron(qn₁, qn₂, signs=(+1, -1)) == CNZ4(-1.0, 3.0)
    @test kron(qn₁, qn₂, signs=(-1, +1)) == CNZ4(+1.0, 1.0)
    @test kron(qn₁, qn₂, signs=(-1, -1)) == CNZ4(-3.0, 3.0)
    @test union(qn₁, qn₂, signs=(+1, +1)) == AbelianNumbers('G', [+qn₁, +qn₂], [0, 1, 2], :indptr)
    @test union(qn₁, qn₂, signs=(+1, -1)) == AbelianNumbers('G', [+qn₁, -qn₂], [0, 1, 2], :indptr)
    @test union(qn₁, qn₂, signs=(-1, +1)) == AbelianNumbers('G', [-qn₁, +qn₂], [0, 1, 2], :indptr)
    @test union(qn₁, qn₂, signs=(-1, -1)) == AbelianNumbers('G', [-qn₁, -qn₂], [0, 1, 2], :indptr)
    @test ⊗(qn₁, qn₂) == kron(qn₁, qn₂)
    @test ⊕(qn₁, qn₂) == union(qn₁, qn₂)

    qns = AbelianNumbers('U', [qn₁, qn₂], [0, 2, 4], :indptr)
    @test +qns == AbelianNumbers('U', [+qn₁, +qn₂], [0, 2, 4], :indptr)
    @test -qns == AbelianNumbers('U', [-qn₁, -qn₂], [0, 2, 4], :indptr)
    @test qns*3 == 3*qns == AbelianNumbers('U', [qn₁*3, qn₂*3], [0, 2, 4], :indptr)
    @test qns^2 == AbelianNumbers('G', [qn₁+qn₁, qn₁+qn₂, qn₁+qn₁, qn₁+qn₂, qn₂+qn₁, qn₂+qn₂, qn₂+qn₁, qn₂+qn₂], [2, 2, 2, 2, 2, 2, 2, 2], :counts)

    @test qns+qn == qn+qns == AbelianNumbers('U', [qn₁+qn, qn₂+qn], [0, 2, 4], :indptr)
    @test qns-qn == AbelianNumbers('U', [qn₁-qn, qn₂-qn], [0, 2, 4], :indptr)
    @test qn-qns == AbelianNumbers('U', [qn-qn₁, qn-qn₂], [0, 2, 4], :indptr)

    qns₁ = AbelianNumbers('U', [+qn₁, -qn₂], [0, 2, 4], :indptr)
    qns₂ = AbelianNumbers('U', [-qn₁, +qn₂], [0, 3, 4], :indptr)
    @test union(qns₁, qns₂, signs=(+1, +1)) == AbelianNumbers('G', [+qn₁, -qn₂, -qn₁, +qn₂], [2, 2, 3, 1], :counts)
    @test union(qns₁, qns₂, signs=(+1, -1)) == AbelianNumbers('G', [+qn₁, -qn₂, +qn₁, -qn₂], [2, 2, 3, 1], :counts)
    @test union(qns₁, qns₂, signs=(-1, +1)) == AbelianNumbers('G', [-qn₁, +qn₂, -qn₁, +qn₂], [2, 2, 3, 1], :counts)
    @test union(qns₁, qns₂, signs=(-1, -1)) == AbelianNumbers('G', [-qn₁, +qn₂, +qn₁, -qn₂], [2, 2, 3, 1], :counts)
    @test kron(qns₁, qns₂, signs=(+1, +1)) == AbelianNumbers('G', [+qn₁-qn₁, +qn₁+qn₂, +qn₁-qn₁, +qn₁+qn₂, -qn₂-qn₁, -qn₂+qn₂, -qn₂-qn₁, -qn₂+qn₂], [3, 1, 3, 1, 3, 1, 3, 1], :counts)
    @test kron(qns₁, qns₂, signs=(+1, -1)) == AbelianNumbers('G', [+qn₁+qn₁, +qn₁-qn₂, +qn₁+qn₁, +qn₁-qn₂, -qn₂+qn₁, -qn₂-qn₂, -qn₂+qn₁, -qn₂-qn₂], [3, 1, 3, 1, 3, 1, 3, 1], :counts)
    @test kron(qns₁, qns₂, signs=(-1, +1)) == AbelianNumbers('G', [-qn₁-qn₁, -qn₁+qn₂, -qn₁-qn₁, -qn₁+qn₂, +qn₂-qn₁, +qn₂+qn₂, +qn₂-qn₁, +qn₂+qn₂], [3, 1, 3, 1, 3, 1, 3, 1], :counts)
    @test kron(qns₁, qns₂, signs=(-1, -1)) == AbelianNumbers('G', [-qn₁+qn₁, -qn₁-qn₂, -qn₁+qn₁, -qn₁-qn₂, +qn₂+qn₁, +qn₂-qn₂, +qn₂+qn₁, +qn₂-qn₂], [3, 1, 3, 1, 3, 1, 3, 1], :counts)
    @test ⊕(qns₁, qns₂) == union(qns₁, qns₂)
    @test ⊗(qns₁, qns₂) == kron(qns₁, qns₂)
end

@testset "sort" begin
    qn₁, qn₂, qn₃, qn₄ = CNZ4(3.0, 2.0), CNZ4(4.0, 1.0), CNZ4(4.0, 2.0), CNZ4(3.0, 3.0)
    oldqns = AbelianNumbers('G', [qn₁, qn₂, qn₃, qn₄, qn₂, qn₃, qn₁, qn₄], [2, 1, 3, 4, 1, 1, 2, 3], :counts)
    newqns, permutation = sort(oldqns)
    @test newqns == AbelianNumbers('C', [qn₁, qn₄, qn₂, qn₃], [4, 7, 2, 4], :counts)
    @test permutation == [1, 2, 13, 14, 7, 8, 9, 10, 15, 16, 17, 3, 11, 4, 5, 6, 12]
end

@testset "findall" begin
    qn₁, qn₂, qn₃ = CNZ4(2.0, 3.0), CNZ4(1.0, 2.0), CNZ4(0.0, 0.0)
    qns = AbelianNumbers('C', [qn₂, qn₁], [2, 3], :counts)
    @test findall(qn₁, qns, :compression) == [2]
    @test findall(qn₂, qns, :compression) == [1]
    @test findall(qn₃, qns, :compression) == []
    @test findall(qn₁, qns, :expansion) == [3, 4, 5]
    @test findall(qn₂, qns, :expansion) == [1, 2]
    @test findall(qn₃, qns, :expansion) == []

    qns = AbelianNumbers('G', [qn₁, qn₂, qn₁], [2, 2, 1], :counts)
    @test findall(qn₁, qns, :compression) == [1, 3]
    @test findall(qn₂, qns, :compression) == [2]
    @test findall(qn₃, qns, :compression) == []
    @test findall(qn₁, qns, :expansion) == [1, 2, 5]
    @test findall(qn₂, qns, :expansion) == [3, 4]
    @test findall(qn₃, qns, :expansion) == []
end

@testset "filter" begin
    qn₁, qn₂ = CNZ4(1.0, 2.0), CNZ4(3.0, 0.0)
    qns = AbelianNumbers('G', [qn₁, qn₂, qn₁, qn₂], [1, 2, 3, 4], :counts)
    @test filter((qn₁, qn₂), qns) == AbelianNumbers('G', [qn₁, qn₂, qn₁, qn₂], [1, 2, 3, 4], :counts)
    @test filter(qn₂, qns) == AbelianNumbers('G', [qn₂, qn₂], [2, 4], :counts)
end

@testset "prod" begin
    qn₁, qn₂ = CNZ4(1.0, 1.0), CNZ4(2.0, 3.0)
    qns₁ = AbelianNumbers('U', [qn₂, qn₁], [4, 5], :counts)
    qns₂ = AbelianNumbers('U', [qn₁, qn₂], [2, 3], :counts)
    qns, records = prod(qns₁, qns₂, signs=(+1, -1))
    @test qns == AbelianNumbers('C', [qn₁-qn₂, zero(CNZ4), qn₂-qn₁], [15, 22, 8], :counts)
    @test records == Dict(  (qn₁-qn₂) =>  Dict((qn₁, qn₂)=>1:15),
                            zero(CNZ4) => Dict((qn₂, qn₂)=>1:12, (qn₁, qn₁)=>13:22),
                            (qn₂-qn₁) =>  Dict((qn₂, qn₁)=>1:8)
                            )
end

@testset "expand" begin
    qn = CNZ4(3.0, 2.0)
    qns = AbelianNumbers(qn, 4)
    @test expand(qns, :indices) == [1, 1, 1, 1]
    @test expand(qns, :contents) == [qn, qn, qn, qn]
end

@testset "decompose" begin
    qn₁, qn₂ = CNZ4(0.0, 0.0), CNZ4(1.0, 1.0)
    qns = AbelianNumbers('U', [qn₁, qn₂], [0, 1, 2], :indptr)
    target = CNZ4(2.0, 2.0)
    result = Set(((1, 1, 2, 2), (1, 2, 1, 2), (1, 2, 2, 1), (2, 1, 1, 2), (2, 1, 2, 1), (2, 2, 1, 1)))
    @test ⊆(Set(decompose(target, qns, qns, qns, qns; signs=(1, 1, 1, 1), nmax=10, method=:bruteforce)), result)
    @test ⊆(Set(decompose(target, qns, qns, qns, qns; signs=(1, 1, 1, 1), nmax=10, method=:montecarlo)), result)
end

@testset "permute" begin
    qn₁, qn₂, qn₃ = CNZ4(1.0, 2.0), CNZ4(3.0, 0.0), CNZ4(4.0, 1.0)
    qns = AbelianNumbers('G', [qn₁, qn₂, qn₃], [2, 3, 4], :counts)
    @test permute(qns, [3, 2, 1], :compression) == AbelianNumbers('G', [qn₃, qn₂, qn₁], [4, 3, 2], :counts)
    @test permute(qns, [4, 6, 9, 8], :expansion) == AbelianNumbers('G', [qn₂, qn₃, qn₃, qn₃], [1, 1, 1, 1], :counts)
end

@testset "ConcreteAbelianNumbers" begin
    @test SpinZ|>fieldnames == (:Sz,)
    @test SpinZ|>periods == (Inf,)
    @test ParticleNumber|>fieldnames == (:N,)
    @test ParticleNumber|>periods == (Inf,)
    @test SpinfulParticle|>fieldnames == (:N, :Sz)
    @test SpinfulParticle|>periods == (Inf, Inf)
    @test spinzs(0.5) == AbelianNumbers('C', [SpinZ(-0.5), SpinZ(0.5)], [0, 1, 2], :indptr)
    @test particlenumbers(1.0) == AbelianNumbers('C', [ParticleNumber(0.0), ParticleNumber(1.0)], [0, 1, 2], :indptr)
    @test spinfulparticles(0.5) == AbelianNumbers('C', [SpinfulParticle(0.0, 0.0), SpinfulParticle(1.0, -0.5), SpinfulParticle(1.0, 0.5), SpinfulParticle(2.0, 0.0)], [0, 1, 2, 3, 4], :indptr)

    @test Momentum₁{10} |> periods == (10,)
    @test Momentum₂{10, 15} |> periods == (10, 15)
    @test Momentum₃{10, 15, 20} |> periods == (10, 15, 20)

    @test Momentum₁{10}(1) == Momentum₁{10}(11) == Momentum₁{10}(-9)
    @test Momentum₂{10}(1, 1) == Momentum₂{10}(11, 11) == Momentum₂{10}(-9, -9)
    @test Momentum₃{10}(1, 1, 1) == Momentum₃{10}(11, 11, 11) == Momentum₃{10}(-9, -9, -9)
    @test Momentum₂{10, 20}(1, 1) == Momentum₂{10, 20}(11, 21) == Momentum₂{10, 20}(-9, -19)
    @test Momentum₃{10, 20, 30}(1, 1, 1) == Momentum₃{10, 20, 30}(11, 21, 31) == Momentum₃{10, 20, 30}(-9, -19, -29)

    @test Int(Momentum₁{10}(2)) == 3
    @test Int(Momentum₂{10, 20}(2, 3)) == 44
    @test Int(Momentum₃{10, 20, 30}(2, 3, 4)) == 1295

    @test expand(Momentum₁{10}(1), [[1.0, 0.0, 0.0]]) == [0.1, 0.0, 0.0]
    @test expand(Momentum₂{10, 100}(1, 1), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]) == [0.1, 0.01, 0.0]
    @test expand(Momentum₃{10, 100, 1000}(1, 1, 1), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) == [0.1, 0.01, 0.001]
end


text:
module Test1
    function test()
        println("Test1: test")
    end
    function test2()
        println("Test1: test2")
    end
end

max_stars_repo_name: srz-zumix/wandbox-api

max_issues_repo_name: srz-zumix/wandbox-api

max_forks_repo_name: srz-zumix/wandbox-api


text:
using InvertibleNetworks, NNlib, LinearAlgebra, Test

# Data
nx = 128
ny = 128
nz = 128
n_in = 3
batchsize = 4
k = 3   # kernel size
s = 1   # stride
p = 1   # padding

###################################################################################################
# Hyperbolic layer

# Data
X_prev = zeros(Float32, nx, ny, nz, n_in, batchsize)
X_curr = randn(Float32, nx, ny, nz, n_in, batchsize)

# Layer
HL = HyperbolicLayer(nx, ny, nz, n_in, batchsize, k, s, p; action="down", α=1f-1, hidden_factor=2)

Y_curr, Y_new = HL.forward(X_prev, X_curr)
X_prev_, X_curr_ = HL.inverse(Y_curr, Y_new)

# Test invertibility
@test isapprox(norm(X_prev_)/prod(size(X_prev_)), 0f0; atol=1f-6)
@test isapprox(norm(X_curr - X_curr_)/norm(X_curr_), 0f0; atol=1f-6)


text:
# ------------------------------------------------------------------
# During time evolution, we may be concerned that though we do have
# convergence of the function in the L^2 norm, we may have a failure
# of the function's derivative in the sparse basis to correctly
# converge
#
# Consequently, we develop a simple filtering method which we can
# apply at every
# ------------------------------------------------------------------

# Efficiency criticality: LOW

# Accuracy criticality: MEDIUM

function filter(D::Int, k::Int, n::Int; scheme="sparse", p=5, alpha=15)

    cutoff    = get_cutoff(scheme, D, n)
    modes    = ntuple(i-> k, D)
    ls        = ntuple(i->(n+1),D)

    V = Float64[]
    for level in CartesianIndices(ls) #This really goes from 0 to l_i for each i
        cutoff(level) && continue

        cells = ntuple(i -> 1<<max(0, level[i]-2), D)
        depth = sum
        depth = (1<<sum((level - 1).I))/(1<<(n+1))
        for cell in CartesianIndices(cells)
            for mode in CartesianIndices(modes)
                val = filter(depth; p=p, alpha=alpha)
                push!(V, val)
            end
        end
    end
    return spdiagm(0 => V)
end


text:
const DATADIR = joinpath(dirname(pathof(@__MODULE__)), "..", "data")

# US Equities from A to E (so far)

for i ∈ [:AAPL, :BA, :CAT, :DELL, :EBAY, :F, :GE, :TX]
  ex = :(
    const $i = readtimearray(joinpath(DATADIR, $("$i.csv")), meta = $("$i")) |> cleanup_colname!)
  eval(ex)
end

# smaller datasets used for testing time-related packages
const dates = Date(2000,1,1):Day(1):Date(2001,12,31)
const cl        = AAPL[:Close][dates]
const op        = AAPL[:Open][dates]
const ohlc      = AAPL[:Open, :High, :Low, :Close][dates]
const ohlcv     = AAPL[:Open, :High, :Low, :Close, :Volume][dates]
const datetime1 = readtimearray(joinpath(DATADIR, "datetime1.csv"))
const datetime2 = readtimearray(joinpath(DATADIR, "datetime2.csv"),
                                format="yyyy-mm-dd HH:MM:SS")
const mdata     = TimeArray(timestamp(cl), values(cl), colnames(cl), "Apple") # data with meta field occupied


text:
#=
- Julia version: 1.1.0
- Author: Jirik
=#
# # using Revise
using Test
using Logging
using LarSurf

# @testset "Init and deinit" begin
# # Dataset can be downloaded from here:
# # https://www.ircad.fr/research/3d-ircadb-01/
# # this is path where data on all my computers are stored
# 	dicomdir = homedir() * "/data/medical/orig/3Dircadb1.1/MASKS_DICOM/portalvein/"
#
# 	datap = LarSurf.read_dicom_dir(dicomdir, "image_")
# 	data3d = datap["data3d"]
# 	voxelsize_mm = datap["voxelsize_mm"]
# 	@test size(data3d, 2) == 512
# 	@test size(data3d, 3) == 512
# 	@test voxelsize_mm[2] == voxelsize_mm[3]
#
# end


text:
function createSonumList(::Type{T},nbit::Int) where {T<:AbstractFloat}
    return zeros(T,2^(nbit-1)+1)
end

function createSonumTable(::Type{T},nbit::Int) where {T<:AbstractFloat}
    n = 2^(nbit-1)+1
    A = Array{T,2}(undef,n,n)
    return Symmetric(A),A
end


text:
module CSDD

using LogicCircuits
# using ..ProbabilisticCircuits
using ..Utils

export

#CredalSDDs
estimate_credal_parameters, estimate_credal_parameters_cached,
CredalΔNode, CredalΔ, CredalLeafNode, CredalInnerNode,
CredalLiteral, Credal⋀, Credal⋁, CredalCache,
# , variable, num_parameters, compute_log_likelihood,
log_marginal_upper,
log_marginal_lower,
log_prob_lower,
log_prob_upper,
conditional_lower,
conditional_upper
# log_proba_lower
# ProbCircuits
# ProbΔNode, ProbΔ, ProbΔ, ProbLeafNode, ProbInnerNode,
# ProbLiteral, Prob⋀, Prob⋁, ProbCache, variable, num_parameters, compute_log_likelihood,
# log_proba,
# log_likelihood, estimate_parameters, log_likelihood_per_instance, marginal_log_likelihood_per_instance,
# initial_mixture_model, estimate_parameters_from_aggregates, compute_ensemble_log_likelihood,
# expectation_step, maximization_step, expectation_step_batch, train_mixture_with_structure, check_parameter_integrity,
# ll_per_instance_per_component, ll_per_instance_for_ensemble,estimate_parameters_cached,
# sample,
# MPE, MAP,prob_origin, copy_node, conjoin_like, disjoin_like, literal_like, normalize, replace_node,

# ProbFlowCircuits
# marginal_pass_up, marginal_pass_down, marginal_pass_up_down,

# Mixtures
# Mixture, AbstractFlatMixture, FlatMixture, FlatMixtureWithFlow,component_weights,FlatMixtureWithFlows,
# log_likelihood, log_likelihood_per_instance, log_likelihood_per_instance_component,
# init_mixture_with_flows, reset_mixture_aggregate_flows, aggregate_flows, estimate_parameters,
# AbstractMetaMixture, MetaMixture,AbstractFlatMixture,AbstractMixture, components, num_components,

# EM Learner
# train_mixture,

# Bagging
# bootstrap_samples_ids, learn_mixture_bagging, learn_mixture_bagging2,
# init_bagging_samples, train_bagging,

# # VtreeLearner
# MetisContext, metis_top_down, BlossomContext, blossom_bottom_up!,
# test_top_down, test_bottom_up!,learn_vtree_bottom_up,

# MutualInformation
# mutual_information, DisCache, conditional_entropy, sum_entropy_given_x,

# Clustering
# clustering,

# Queries
# pr_constraint, psdd_entropy, psdd_kl_divergence

include("CredalCircuits.jl")
include("CredalFlowCircuits.jl")
include("VtreeLearner.jl")
include("Queries.jl")

end


text:
using Random

"""
A pile of tiles.
An abstraction for a collection of Mahjong tiles, similar to a deque.

Some properties:
- ordered and unsorted
- size only strictly decreases once created
"""
struct TilePile
    tiles::Vector{SuitedTile}
end

@doc """
    TilePile(tiles)
    TilePile(name)

Create a particular predefined/standard `TilePile` given the `name`, or
create the `TilePile` from a specific list of tiles.

## Supported Names

### `:standard`

4 of each numbered tile from the 3 suits, plus 4 of each wind, 4 of each dragon, and the 8 flowers/seasons

"""
function TilePile(type::Symbol)
    if type == :standard
        tiles = [
            SuitedTile(suit, i)
            for suit in all_suits
            for i in 1:max_number(suit)
            for _ in 1:typical_duplication(suit)
        ]
    else
        @warn "Not a recognized type of Mahjong game (don't know how to create a pile for it)"
        tiles = []
    end

    return TilePile(tiles)
end

Base.length(tp::TilePile) = length(tp.tiles)
Base.first(tp::TilePile) = first(tp.tiles)
Base.first(tp::TilePile, n::Integer) = first(tp.tiles, n)
Base.last(tp::TilePile) = last(tp.tiles)
Base.last(tp::TilePile, n::Integer) = last(tp.tiles, n)
Base.iterate(tilepile::TilePile) = iterate(tilepile.tiles)
Base.iterate(tilepile::TilePile, state) = iterate(tilepile.tiles, state)

"""
    popfirst!(tilepile)

Remove the "first" tile from the tile pile.
Removing from this end is intended to be used for flower/four-of-a-kind replacement
(this is to improve average performance).
"""
Base.popfirst!(tp::TilePile) = popfirst!(tp.tiles)

"""
    popfirst!(tilepile, num_tiles)

Remove the "first" `num_tiles` tiles from the tile pile as an array.
Removing from this end is intended to be used for flower/four-of-a-kind replacement
(this is to improve average performance).
"""
function Base.popfirst!(tp::TilePile, n::Integer)
    n < 0 && throw("cannot pop fewer than 0 tiles")

    return splice!(tp.tiles, 1:min(n, lastindex(tp.tiles)))
end

"""
    pop!(tilepile)

Remove a tile from the end of the tile pile.
Removing from this end is used for the usual pickup of tiles from the deck
(not for flower/four-of-a-kind replacement).
"""
Base.pop!(tp::TilePile) = pop!(tp.tiles)

"""
    pop!(tilepile, num_tiles)

Remove `num_tiles` tiles from the end of the tile pile as an array.
Removing from this end is used for the usual pickup of tiles from the deck
(not for flower/four-of-a-kind replacement).
"""
function Base.pop!(tp::TilePile, n::Integer)
    n < 0 && throw("cannot pop fewer than 0 tiles")

    return splice!(tp.tiles, max(1, lastindex(tp.tiles) - n + 1):lastindex(tp.tiles))
end
Random.shuffle!(tp::TilePile) = (shuffle!(tp.tiles); tp)
Random.shuffle!(rng, tp::TilePile) = (shuffle!(rng, tp.tiles); tp)

Base.circshift(tp::TilePile, shift) = TilePile(circshift(tp.tiles, shift))

function Base.show(io::IO, mime::MIME"text/plain", tilepile::TilePile)
    for tile in tilepile.tiles
        show(io, mime, tile)
    end
end
function Base.show(io::IO, tilepile::TilePile)
    print(io, "TilePile([")
    for (index, tile) in enumerate(tilepile)
        print(io, tile)
        index != length(tilepile) && print(io, ", ")
    end
    print(io, "])")
end

function deal_tiles(tp::TilePile; quadruple_rounds=3, single_rounds=1)
    player_hands = [SuitedTile[] for _ in 1:4]
    for _ in 1:quadruple_rounds, player in 1:4
        append!(player_hands[player], pop!(tp, 4))
    end
    for _ in 1:single_rounds, player in 1:4
        push!(player_hands[player], pop!(tp))
    end
    push!(first(player_hands), pop!(tp))

    return Hand.(player_hands)
end


text:
using eeg

function main
# julia code
end

function side : main
# julia code
end


text:
"""
`module  JAC.ReducedDensityMatrix`
    ... a submodel of JAC that contains all methods for computing reduced density matrices, natural orbitals, radial density
        distributions and related information for some level(s).
"""
module ReducedDensityMatrix

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

    """
    `struct  ReducedDensityMatrix.Settings  <:  AbstractPropertySettings`
        ... defines a type for providing and parameters for reduced density matrices and other parameters.

        + calcNatural              ::Bool             ... True if natural orbitals need to be calculated, and false otherwise.
        + calcDensity              ::Bool             ... True if the radial density need to be calculated, and false otherwise.
        + calcIpq                  ::Bool             ... True if orbital interaction I_pq need to be calculated, and false otherwise.
        + calc2pRDM                ::Bool             ... True if the 2p RDM need to be calculated, and false otherwise;
                                                          the 1pRDM is calculated by default.
        + printBefore              ::Bool             ... True if a list of selected levels is printed before the actual computations start.
        + levelSelection           ::LevelSelection   ... Specifies the selected levels, if any.
    """
    struct Settings  <:  AbstractPropertySettings
        calcNatural                ::Bool
        calcDensity                ::Bool
        calcIpq                    ::Bool
        calc2pRDM                  ::Bool
        printBefore                ::Bool
        levelSelection             ::LevelSelection
    end

    """
    `ReducedDensityMatrix.Settings(; calcNatural::Bool=true, calcDensity::Bool=true, calcIpq::Bool=false,
                                     printBefore::Bool=true, levelSelection::LevelSelection=LevelSelection())
        ... keyword constructor to overwrite selected value of reduced-density matrix computations.
    """
    function Settings(; calcNatural::Bool=true, calcDensity::Bool=true, calcIpq::Bool=false,
                                     printBefore::Bool=true, levelSelection::LevelSelection=LevelSelection())
        Settings(calcNatural, calcDensity, calcIpq, printBefore, levelSelection)
    end

    # `Base.show(io::IO, settings::ReducedDensityMatrix.Settings)`
    #       ... prepares a proper printout of the variable settings::ReducedDensityMatrix.Settings.
    function Base.show(io::IO, settings::ReducedDensityMatrix.Settings)
        println(io, "calcNatural:              $(settings.calcNatural)  ")
        println(io, "calcDensity:              $(settings.calcDensity)  ")
        println(io, "calcIpq:                  $(settings.calcIpq)  ")
        println(io, "calc2pRDM:                $(settings.calc2pRDM)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `struct  ReducedDensityMatrix.Outcome`
        ... defines a type to keep the outcome of a reduced density-matrix and natural-orbital computation.

        + level                     ::Level                 ... Atomic level to which the outcome refers to.
        + naturalSubshells          ::Array{Subshell,1}     ... List of natural orbitals (subshells).
        + naturalOccupation         ::Array{Float64,1}      ... Occupation numbers of natural orbitals.
        + naturalOrbitalExpansion   ::Dict{Subshell, Array{Float64,1}}
            ... Dictionary of the expansion of (one-electron) natural orbitals in terms of the standard orbitals
                from the given basis.
        + naturalOrbitals           ::Dict{Subshell, Orbital} . Dictionary of (one-electron) natural orbitals.
        + orbitalInteraction        ::Array{Float64,2}      ... Orbital interaction I_pq = I[ip,iq].
        + rho1p                     ::Array{Float64,2}      ... One-particle RDM rho^(1p) [ip,iq].
        + rho2p                     ::Array{Float64,4}      ... Two-particle RDM rho^(1p) [ip,iq].
        + electronDensity           ::Radial.Density        ... Radial density distribution.
    """
    struct Outcome
        level                       ::Level
        naturalSubshells            ::Array{Subshell,1}
        naturalOccupation           ::Array{Float64,1}
        naturalOrbitalExpansion     ::Dict{Subshell, Array{Float64,1}}
        naturalOrbitals             ::Dict{Subshell, Orbital}
        orbitalInteraction          ::Array{Float64,2}
        rho1p                       ::Array{Float64,2}
        rho2p                       ::Array{Float64,4}
        electronDensity             ::Radial.Density
    end

    """
    `ReducedDensityMatrix.Outcome()`  ... constructor for an `empty` instance of Hfs.Outcome for the computation of isotope-shift properties.
    """
    function Outcome()
        Outcome(Level(), Subshell[], Float64[], Dict{Subshell, Array{Float64,1}}, Dict{Subshell, Orbital}(),
                zeros(1,1), zeros(1,1), zeros(1,1,1,1), Radial.Density() )
    end

    # `Base.show(io::IO, outcome::ReducedDensityMatrix.Outcome)`  ... prepares a proper printout of the variable outcome::ReducedDensityMatrix.Outcome.
    function Base.show(io::IO, outcome::ReducedDensityMatrix.Outcome)
        println(io, "level:                   $(outcome.level)  ")
        println(io, "naturalSubshells:        $(outcome.naturalSubshells)  ")
        println(io, "naturalOccupation:       $(outcome.naturalOccupation)  ")
        println(io, "naturalOrbitalExpansion: $(outcome.naturalOrbitalExpansion)  ")
        println(io, "naturalOrbitals:         $(outcome.naturalOrbitals)  ")
        println(io, "orbitalInteraction:      $(outcome.orbitalInteraction)  ")
        println(io, "rho1p:                   $(outcome.rho1p)  ")
        println(io, "rho2p:                   $(outcome.rho2p)  ")
        println(io, "electronDensity:         $(outcome.electronDensity)  ")
    end

    """
    `ReducedDensityMatrix.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                          settings::ReducedDensityMatrix.Settings; output=true)`
        ... to compute (as selected) the 1-particle and 2-particle reduced density matrices, natural orbitals or other
            requested information for the levels of interest.
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                             settings::ReducedDensityMatrix.Settings; output=true)
        println("")
        printstyled("ReducedDensityMatrix.computeOutcomes(): The computation of the reduced density matrices starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------------ \n", color=:light_green)
        #
        outcomes = ReducedDensityMatrix.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    ReducedDensityMatrix.displayOutcomes(outcomes, settings)    end
        # Calculate all requested properties for each outcome
        newOutcomes = ReducedDensityMatrix.Outcome[]
        for  outcome in outcomes
            newOutcome = ReducedDensityMatrix.computeProperties(outcome, nm, grid, settings)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        ReducedDensityMatrix.displayResults(stdout, newOutcomes, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    ReducedDensityMatrix.displayResults(iostream, newOutcomes, settings)   end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `ReducedDensityMatrix.computeNaturalOrbitalExpansion(naturalSubshells::Array{Subshell,1}, rho1p::Array{Subshell,2}, level::Level)`
        ... to perform the expansion of the natural orbitals for the given level in terms of the standard orbitals
            as defined by the basis. A tuple (naturalOcc::Array{Float64,1}, naturalExp::Dict{Subshell, Array{Float64,1}})
            is returned which provides the natural occpuation numbers and expansion coefficients with regard
            to the given list of natural orbitals. This method makes use of angular coefficients for a zero-rank operator
            to calculate and diagonalize the expansion matrix.
    """
    function  computeNaturalOrbitalExpansion(naturalSubshells::Array{Subshell,1}, rho1p::Array{Subshell,2}, level::Level)
        naturalOcc = Float64[];      naturalExp = Dict{Subshell, Array{Float64,1}}();       lenNO = length(naturalSubshells)
        @show naturalSubshells
        @show rho_1p
        # Diagonalize the 1p RDM
        #
        for i = 1:lenNO   push!(naturalOcc, rho_pq[i,i])
                        naturalExp[ naturalSubshells[i] ] = vector
        end

        return( naturalOcc, naturalExp )
    end

    """
    `ReducedDensityMatrix.computeNaturalOrbitals(naturalSubshells::Array{Subshell,1}, naturalExp::Dict{Subshell, Array{Float64,1}}, level::Level)`
        ... to compute the natural orbitals as superposition of the standard orbitals; for each subshell in naturalExp,
            an orbital is computed, and naturalOrbitals::Dict{Subshell, Orbital} returned
    """
    function  computeNaturalOrbitals(naturalSubshells::Array{Subshell,1}, naturalExp::Dict{Subshell, Array{Float64,1}}, level::Level)
        naturalOrbitals = Dict{Subshell, Orbital}()
        println(">> computeNaturalOrbitals() ... not yet implemented !!")

        return( naturalOrbitals )
    end

    """
    `ReducedDensityMatrix.computeOrbitalInteractions(naturalOrbitals::Dict{Subshell, Orbital}, level::Level)`
        ... to compute the natural orbitals as superposition of the standard orbitals; for each subshell in naturalExp,
            an orbital is computed, and naturalOrbitals::Dict{Subshell, Orbital} returned
    """
    function  computeOrbitalInteractions(naturalOrbitals::Dict{Subshell, Orbital}, level::Level)
        orbitalInteraction = zeros(2,2)
        println(">> computeOrbitalInteractions() ... not yet implemented !!")

        return( orbitalInteraction )
    end

    """
    `ReducedDensityMatrix.computeProperties(outcome::ReducedDensityMatrix.Outcome, nm::Nuclear.Model, grid::Radial.Grid,
                                            settings::ReducedDensityMatrix.Settings)
        ... to compute all properties for a given level; an outcome::ReducedDensityMatrix.Outcome is returned for
            which the properties are now evaluated explicitly.
    """
    function  computeProperties(outcome::ReducedDensityMatrix.Outcome, nm::Nuclear.Model, grid::Radial.Grid,
                                settings::ReducedDensityMatrix.Settings)
        naturalOccupation = outcome.naturalOccupation;        naturalOrbitalExp  = outcome.naturalOrbitalExpansion;
        naturalOrbitals   = outcome.naturalOrbitals;          orbitalInteraction = outcome.orbitalInteraction;
        electronDensity   = outcome.electronDensity;          rho1p              = outcome.rho1p;             rho2p = outcome.rho2p

        naturalSubshells   = copy(outcome.level.basis.subshells)

        if  settings.calcNatural
            rho1p                                = compute1pRDM(naturalSubshells, outcome.level)
            naturalOccupation, naturalOrbitalExp = computeNaturalOrbitalExpansion(naturalSubshells, rho1p, outcome.level)
            naturalOrbitals                      = computeNaturalOrbitals(naturalSubshells, naturalOrbitalExp, outcome.level)
        end

        if  settings.calcDensity
            electronDensity = computeRadialDistribution(naturalOrbitalExp, grid, outcome.level)
        end

        if  settings.calcIpq
            orbitalInteraction = computeOrbitalInteractions(naturalOrbitals, outcome.level)
        end

        if  settings.calc2pRDM
            rho2p = compute2pRDM(naturalSubshells, outcome.level)
        end

        newOutcome = ReducedDensityMatrix.Outcome( outcome.level, naturalSubshells,  naturalOccupation, naturalOrbitalExp,
                                                   naturalOrbitals, orbitalInteraction, rho1p, rho2p, electronDensity )
        return( newOutcome )
    end

    """
    `ReducedDensityMatrix.computeRadialDistribution(naturalExp::Dict{Subshell, Array{Float64,1}}, grid::Radial.Grid, level::Level)`
        ... to compute the radial distribution function from the given standard orbitals and natural expansion
            coefficients. An radial disribution D::Array{Float64,1} is returned that specifies the radial distribution
            w.r.t. the given grid.
    """
    function  computeRadialDistribution(naturalExp::Dict{Subshell, Array{Float64,1}}, grid::Radial.Grid, level::Level)
        D = Radial.Density()
        println(">> computeRadialDistribution() ... not yet implemented !!")

        return( D )
    end

    """
    `ReducedDensityMatrix.compute1pRDM(naturalSubshells::Array{Subshell,1}, level::Level)`
        ... to compute 1p RDM; a rdm::Array{Float64,2} is returned.
    """
    function  compute1pRDM(naturalSubshells::Array{Subshell,1}, level::Level)
        lenNO = length(naturalSubshells);    rho_pq     = zeros(lenNO,lenNO)
        # Cycle over all matrix elements of the CSF basis
        subshellList = level.basis.subshells
        opa          = SpinAngular.OneParticleOperator(0, plus, true)
        for  (ir, rcsf) in enumerate(level.basis.csfs)
            for  (is, scsf) in enumerate(level.basis.csfs)
                # Calculate angular coefficient for rank-0 operator
                wa = SpinAngular.computeCoefficients(opa, rcsf, scsf, subshellList)
                # Cycle over the pair of natural subshells in rho_pq
                for (ip,p)  in  enumerate(naturalSubshells)
                    for (iq,q)  in  enumerate(naturalSubshells)
                        for  coeff in wa
                            if  (p == coeff.a   &&  q == coeff.b)  ||  (p == coeff.b   &&  q == coeff.a)
                                jj = Basics.subshell_2j(level.basis.orbitals[coeff.a].subshell)
                                rho_pq[ip,iq] = rho_pq[iq,ip] = rho_pq[ip,iq] + level.mc[ir] * coeff.T * sqrt( jj + 1) * level.mc[is]
                            end
                        end
                    end
                end
            end
        end
        @show rho_pq

        return( rho_pq )
    end

    """
    `ReducedDensityMatrix.compute2pRDM(p::Subshell, q::Subshell, r::Subshell, s::Subshell, naturalSubshells::Array{Subshell,1},
                                       naturalExp::Dict{Subshell, Array{Float64,1}}, level::Level)`
        ... to compute 2p RDM for the given pair (p,q;r,s) of subshells; a rdm::Array{Float64,4} is returned.
    """
    function  compute2pRDM(p::Subshell, q::Subshell, r::Subshell, s::Subshell, naturalSubshells::Array{Subshell,1},
                           naturalExp::Dict{Subshell, Array{Float64,1}}, level::Level)
        ns = length(naturalSubshells);  rdm = zeros(ns,ns,ns,ns)
        println(">> compute2pRDM() ... not yet implemented !!")

        return( rdm )
    end

    """
    `ReducedDensityMatrix.determineOutcomes(multiplet::Multiplet, settings::ReducedDensityMatrix.Settings)`
        ... to determine a list of Outcomes's for the computation of reduced density matrices and other information
            for the given multiplet. It takes into account the particular selections and settings.
            An Array{ReducedDensityMatrix.Outcome,1} is returned. Apart from the level specification, all physical
            properties are set to zero during the initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::ReducedDensityMatrix.Settings)
        outcomes = ReducedDensityMatrix.Outcome[]
        @show settings.levelSelection
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, ReducedDensityMatrix.Outcome(level, Subshell[], Float64[], Dict{Subshell, Array{Float64,1}}(),
                                        Dict{Subshell, Orbital}(), zeros(1,1), zeros(1,1), zeros(1,1,1,1), Radial.Density()) )
            end
        end
        return( outcomes )
    end

    """
    `ReducedDensityMatrix.displayOutcomes(outcomes::Array{ReducedDensityMatrix.Outcome,1}, settings::ReducedDensityMatrix.Settings)`
        ... to display a list of levels that have been selected for the computations. A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise. Moreover, the
            selected properties and information is printed as well, though not yet calculated.
    """
    function  displayOutcomes(outcomes::Array{ReducedDensityMatrix.Outcome,1}, settings::ReducedDensityMatrix.Settings)
        nx = 43
        println(" ")
        println("  Results for the natural occupation number, natural orbitals and kp RDM are printed in turn for the following levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = Basics.LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        println("\n\n  The following properties are calculated: \n")
        println("  + Natural occupations numbers")
        println("  + Natural orbital expansion in terms of the standard orbitals.")
        if  settings.calcNatural   println("  + Representation of the natural orbitals.")  end
        if  settings.calcDensity   println("  + Radial electron density distribution.")    end
        if  settings.calcIpq       println("  + Orbital interaction I_pq.")                end

        return( nothing )
    end

    """
    `ReducedDensityMatrix.displayResults(stream::IO, outcomes::Array{ReducedDensityMatrix.Outcome,1}, settings::ReducedDensityMatrix.Settings)`
        ... to display a list of levels that have been selected for the computations. A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise. Moreover, the
            selected properties and information is printed as well, though not yet calculated.
    """
    function  displayResults(stream::IO, outcomes::Array{ReducedDensityMatrix.Outcome,1}, settings::ReducedDensityMatrix.Settings)
        nx = 43
        println(stream, " ")
        println(stream, "  Results for the natural occupation number, natural orbitals and kp RDM are printed in turn for the following levels:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = Basics.LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            println(stream,  sa )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        # Now print all selected result for each outcome
        for  outcome in outcomes
            nx = 43
            println(stream, " ")
            println(stream, "  Natural occupation numbers:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
            # Use table.strings ... to split a vector in strings of 10 subshells/numbers ... and which are returned
            # as a tuple.
            # Split  and generate the header of the table
            stuple = ( "aa", "bb")
            for  (i, st)  in  enumerate(stuple)
                if       i == 1  sa = "1****" * st
                else             sa = "     " * st      end
                println(stream, sa)
            end
            println(stream, "  ", TableStrings.hLine(nx))
            # Split and generate the occupation numbers of the table
            stuple = ( "cc", "dd")
            for  (i, st)  in  enumerate(stuple)
                if       i == 1  sa = "1****" * st
                else             sa = "     " * st      end
                println(sa)
            end
            println(stream, "  ", TableStrings.hLine(nx))
            #
            nx = 43
            println(stream, " ")
            println(stream, "  Natural orbital expansion:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(nx))
            # Use table.strings ... to split a vector in strings of 10 subshells/numbers ... and which are returned
            # as a tuple.
            # Split  and generate the header of the table
            stuple = ( "aa", "bb")
            for  (i, st)  in  enumerate(stuple)
                if       i == 1  sa = "1****" * st
                else             sa = "     " * st      end
                println(stream, sa)
            end
            println(stream, "  ", TableStrings.hLine(nx))
            #
            for subsh in outcome.naturalSubshells
                # Split and generate the occupation numbers of the table
                stuple = ( "cc", "dd")
                for  (i, st)  in  enumerate(stuple)
                    if       i == 1  sa = "1****" * st
                    else             sa = "     " * st      end
                end
            end
            println(stream, "  ", TableStrings.hLine(nx))
            #
            #
            if  settings.calcNatural
                println("\n  Natural orbitals are calculated for the following subshells and are available by outcome.naturalOrbitals:\n")
                sa = "  ";      for  subsh in outcome.naturalSubshells    sa = sa * "   $subsh"     end
                println(stream, sa)
            end
            #
            #
            if  settings.calcDensity
                println(stream, "\n  Radial electron distribution:\n")
                for  i = 1:20:outcome.electronDensity.grid.NoPoints
                    println(stream, "   " * ("      "*string(i))[end-6:end] * ")    " *
                                    @sprintf("%.4e", grid.r[i]) * @sprintf("%.4e", outcome.electronDensity[i])  )
                end
            end
            #
            #
            if  settings.calcIpq
                nx = 43
                println(stream, "\n  Orbital interactions I_pq:\n")
                println(stream, "  ", TableStrings.hLine(nx))
                # Use table.strings ... to split a vector in strings of 10 subshells/numbers ... and which are returned
                # as a tuple.
                # Split  and generate the header of the table
                stuple = ( "aa", "bb")
                for  (i, st)  in  enumerate(stuple)
                    if       i == 1  sa = "NO / NO    " * st
                    else             sa = "           " * st      end
                    println(stream, stream, sa)
                end
                println("  ", TableStrings.hLine(nx))
                # Split and generate the occupation numbers of the table
                for  subsh in outcome.naturalSubshells
                    stuple = ( "cc", "dd")
                    for  (i, st)  in  enumerate(stuple)
                        if       i == 1  sa = "1****" * st
                        else             sa = "     " * st      end
                    end
                end
                println(stream, "  ", TableStrings.hLine(nx))
            end
        end  # outcomes
        #
        #
        return( nothing )
    end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ b1226c30-fdc8-11ea-1020-c1a5f67bc6aa
begin
	using Pkg
	Pkg.add("DataFramesMeta")
end

# ╔═╡ 3c45bc76-fdc5-11ea-21e4-d7973210099b
begin
	using DataFrames
	using CSV
	using HTTP
	using DataFramesMeta
	url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
	df = CSV.read(HTTP.get(url).body)
	df = @transform(df, bus_id = 1:size(df,1)) # Create bus id variable
end

# ╔═╡ 1ec897ae-fdc5-11ea-27b5-554b053bc087
md"""
# Introduction

In this problem set, we will explore a simplified version of the Rust (1987, \textit{Econometrica}) bus engine replacement model. Let's start by reading in the data.
"""

# ╔═╡ b2919210-fdc5-11ea-24d2-bba0c930d349
md"""
# Static estimation

Reshape the data into "long" panel format, calling your long dataset *df\_long*. I have included code on how to do this in the *PS5starter.jl*

1. Create bus id variable
"""

# ╔═╡ Cell order:
# ╟─1ec897ae-fdc5-11ea-27b5-554b053bc087
# ╠═b1226c30-fdc8-11ea-1020-c1a5f67bc6aa
# ╠═3c45bc76-fdc5-11ea-21e4-d7973210099b
# ╟─b2919210-fdc5-11ea-24d2-bba0c930d349


text:
type PetscNonlinearImplicitSolver <: DenseImplicitSolver
    " The System the solve will be done for "
    system::System

    " The 'stiffness matrix' "
    mat::PetscMat

    " The RHS (force vector) "
    rhs::PetscVec

    " The solution vector "
    solution::PetscVec

    " Whether or not this Solver has been initialized "
    initialized::Bool

    " The ghosted solution vector.  This is last because we're going to leave it unitialized in the beginning. "
    ghosted_solution::GhostedPetscVec

    PetscNonlinearImplicitSolver(system::System) = new(system, PetscMat(), PetscVec(), PetscVec(), false)
end

" Initializes the matrix and rhs to the correct size "
function initialize!(solver::PetscNonlinearImplicitSolver)
    @assert solver.system.initialized

    n_dofs = solver.system.n_dofs

    local_n_dofs = solver.system.local_n_dofs

    setSize!(solver.mat, m_local=(Int32)(local_n_dofs), n_local=(Int32)(local_n_dofs))
    generateSparsity!(solver.system)
    setPreallocation!(solver.mat, solver.system.local_dof_sparsity, solver.system.off_proc_dof_sparsity)

    setSize!(solver.rhs, n_local=(Int32)(local_n_dofs))
    setSize!(solver.solution, n_local=(Int32)(local_n_dofs))

    findGhostedDofs!(solver.system)
    solver.ghosted_solution = GhostedPetscVec(Int32[(Int32)(dof) for dof in solver.system.ghosted_dofs], n_local=(Int32)(local_n_dofs))

    solver.initialized = true
end

" Update the ghosted_solution vector from the actual solution "
function updateGhostedSolution!(solver::PetscNonlinearImplicitSolver)
    copy!(solver.ghosted_solution, solver.solution)
end

"""
    Solve using the built-in dense matrix/vector types

    `nl_max_its`: Maximum number of nonlinear iterations to attempt before giving up
    `nl_rel_tol`: Relative nonlinear residual L2-Norm drop tolerance
    `nl_abs_tol`: Absolute nonlinear residual L2-Norm tolerance
"""
function solve!(solver::PetscNonlinearImplicitSolver; nl_max_its=10, nl_rel_tol=1e-8, nl_abs_tol=1e-20)
    if !solver.initialized
        initialize!(solver)
    end

    initial_residual = 0.0

    # The thing being solved for
    delta = similar(solver.solution)

    # Newton loops
    for i in 1:nl_max_its
        assembleResidualAndJacobian(solver, solver.system)

        current_residual_norm = norm(solver.rhs)

        # First iteration save off the residual
        if i == 1
            initial_residual = current_residual_norm
        end

        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            println(i, " NL |R|: ", current_residual_norm, ", |R|/|R1|: ", current_residual_norm / initial_residual)
        end

        if current_residual_norm / initial_residual < nl_rel_tol
            if MPI.Comm_rank(MPI.COMM_WORLD) == 0
                println("Relative Tolerance Reached!")
            end
            return
        end

        if current_residual_norm < nl_abs_tol
            if MPI.Comm_rank(MPI.COMM_WORLD) == 0
                println("Absolute Tolerance Reached!")
            end
            return
        end

        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            println("  Starting Solve... ")
        end

        # Solve the system
        ksp = PetscKSP()
        setOperators(ksp, solver.mat)
        scale!(solver.rhs, -1.0)

        startLog(main_perf_log, "linear_solve")
        solve!(ksp, solver.rhs, delta)
        stopLog(main_perf_log, "linear_solve")

        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            println("Done.")
        end

        plusEquals!(solver.solution, delta)
    end

    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        println("Warning!! Solve did NOT converge to within the set tolerances!")
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0b0902da-40f0-11eb-3e10-a1fe1bd1f364
begin
	using PyPlot
	using SpecialFunctions
	using HypergeometricFunctions
end

# ╔═╡ 0b507856-40f1-11eb-244f-f15c1038b96f
begin
	rc("text", usetex = true)
	ion()
end

# ╔═╡ eaadeff2-546b-11eb-1289-9bade75bad45
function J(α, β, n, z)
    T = eltype(z/2)
    tot = zero(T)
    for k in 0:n
        tot += binomial(α + n, k) * binomial(β + n, n - k) * ((z + 1)/2)^k * ((z - 1)/2)^(n - k)
    end
    tot
end;

# ╔═╡ Cell order:
# ╠═0b0902da-40f0-11eb-3e10-a1fe1bd1f364
# ╠═0b507856-40f1-11eb-244f-f15c1038b96f
# ╠═eaadeff2-546b-11eb-1289-9bade75bad45


text:
function retrieve(
    c::Client,
    ::Type{Reaction},
    m::Message,
    emoji::Union{AbstractString, AbstractChar},
)
    return get_reactions(c, m.channel_id, m.id, emoji)
end
retrieve(c::Client, ::Type{Reaction}, m::Message, e::Emoji) = get(Reaction, c, m, e.name)

function create(
    c::Client,
    ::Type{Reaction},
    m::Message,
    emoji::Union{AbstractString, AbstractChar},
)
    return create_reaction(c, m.channel_id, m.id, emoji)
end
create(c::Client, ::Type{Reaction}, m::Message, e::Emoji) = create(Reaction, c, m, e.name)

function delete(c::Client, ::Type{Reaction}, m::Message)
    return delete_all_reactions(c, m.channel_id, m.id)
end
function delete(
    c::Client,
    emoji::Union{AbstractString, AbstractChar},
    m::Message,
    u::User,
)
    return if c.state.user !== nothing && u.id == c.state.user.id
        delete_own_reaction(c, m.channel_id, m.id, emoji)
    else
        delete_user_reaction(c, m.channel_id, m.id, emoji, u.id)
    end
end
delete(c::Client, e::Emoji, m::Message, u::User) = delete(c, m, e.name, u)


text:
using Random, LinearAlgebra, TensorRefinement, Test

Random.seed!(1365)

@testset "Factor orthogonalization" begin
	n = [4,5,6,7]; p,q = 7,19
	d = length(n)
	@testset "Type: $S" for T ∈ (Float64, Float32), S ∈ (T, Complex{T})
		tol = 100*eps(T)
		U = rand(S, p, n..., q); nrm = norm(U)
		@testset "pivot=$pivot, rev=$rev" for pivot ∈ (false, true), rev ∈ (false, true)
			Q,R = factorqr!(copy(U), Val(pivot); rev=rev)
			r,s = factorranks(Q)
			rev || (r = s)
			if rev
				@test factorranks(R) == (p,r)
				@test factorranks(Q) == (r,q)
				@test factorsize(R) == ones(Int, d)
				@test factorsize(Q) == n
			else
				@test factorranks(Q) == (p,r)
				@test factorranks(R) == (r,q)
				@test factorsize(R) == ones(Int, d)
				@test factorsize(Q) == n
			end
			V = factorcontract(Q, R; rev=rev)
			@test norm(V-U)/nrm ≈ 0 atol=tol
			if rev
				Q = reshape(Q, r, prod(n)*q)
				E = Q*Q'-I
				@test norm(E)/r ≈ 0 atol=tol
			else
				Q = reshape(Q, p*prod(n), r)
				E = Q'*Q-I
				@test norm(E)/r^2 ≈ 0 atol=tol
			end
		end
    end
end


text:
# ---------------------------------------------------------------------------
# Box Plot

notch_width(q2, q4, N) = 1.58 * (q4-q2)/sqrt(N)

@recipe function f(
    ::Type{Val{:boxplot}},
    x,
    y,
    z;
    notch = false,
    range = 1.5,
    outliers = true,
    whisker_width = :half,
    percentiles = [0.25, 0.5, 0.75]
)
    # if only y is provided, then x will be UnitRange 1:size(y,2)
    if typeof(x) <: AbstractRange
        if step(x) == first(x) == 1
            x = plotattributes[:series_plotindex]
        else
            x = [getindex(x, plotattributes[:series_plotindex])]
        end
    end
    xsegs, ysegs = Segments(), Segments()
    texts = String[]
    glabels = sort(collect(unique(x)))
    warning = false
    emptywarning = false
    outliers_x, outliers_y = zeros(0), zeros(0)
    bw = plotattributes[:bar_width]
    isnothing(bw) && (bw = 0.8)
    @assert whisker_width == :match || whisker_width == :half || whisker_width >= 0 "whisker_width must be :match, :half, or a positive number"
    ww = whisker_width == :match ? bw :
         whisker_width == :half ? bw / 2 :
         whisker_width
    for (i, glabel) in enumerate(glabels)
        # filter y
        values = y[filter(i -> _cycle(x, i) == glabel, 1:length(y))]

        # compute quantiles
        q1, q2, q3, q4, q5 = quantile(values, Base.range(0, stop = 1, length = 5))

        # ensure that only valid percentiles are requested (also removing the min + max to ensure lower + upper T gets rendered)
        filter!(x -> x > 0 && x < 100, percentiles)

        # if our validated percentiles result in an empty list, use default quartiles
        if isempty(percentiles)
            if !emptywarning
                @warn("Zero valid percentiles were requested. Plotting using quartiles.")
                emptywarning = true # Show the warning only one time
            end
            percentiles = [0.25, 0.5, 0.75]
        end

        # compute percentiles
        ps = sort(quantile(values, unique(percentiles)))

        # get minimum and maximum values which are not q1 or q5
        pmin, pmax = extrema(ps)

        # notch
        n = notch_width(q2, q4, length(values))

        # warn on inverted notches?
        if notch && !warning && ((q2 > (q3 - n)) || (q4 < (q3 + n)))
            @warn("Boxplot's notch went outside hinges. Set notch to false.")
            warning = true # Show the warning only one time
        end

        # make the shape
        center = Plots.discrete_value!(plotattributes[:subplot][:xaxis], glabel)[1]
        hw = 0.5_cycle(bw, i) # Box width
        HW = 0.5_cycle(ww, i) # Whisker width
        l, m, r = center - hw, center, center + hw
        lw, rw = center - HW, center + HW

        # internal nodes for notches
        L, R = center - 0.5 * hw, center + 0.5 * hw

        # outliers
        if Float64(range) != 0.0  # if the range is 0.0, the whiskers will extend to the data
            limit = range * (pmax - pmin)
            inside = Float64[]
            for value in values
                if (value < (pmin - limit)) || (value > (pmax + limit))
                    if outliers
                        push!(outliers_y, value)
                        push!(outliers_x, center)
                    end
                else
                    push!(inside, value)
                end
            end
            # change q1 and q5 to show outliers
            # using maximum and minimum values inside the limits
            q1, q5 = Plots.ignorenan_extrema(inside)
            q1, q5 = (min(q1, pmin), max(pmax, q5)) # whiskers cannot be inside the box
        end
        # Box
        push!(xsegs, m, lw, rw, m, m)       # lower T
        push!(ysegs, q1, q1, q1, q1, pmin)    # lower T
        push!(
            texts,
            "Lower fence: $q1",
            "Lower fence: $q1",
            "Lower fence: $q1",
            "Lower fence: $q1",
            "Q1: $pmin",
            "",
        )

        if notch
            # draw segments between requested percentile pairs, accounting for case where notch lies between percentiles
            for (ii, (pa, pb)) in enumerate(zip(ps, ps[2:end]))
                # determine if percentile window being processed falls within the notch
                inlower = (q3 - n > pa) && (q3 - n < pb)
                inupper = (q3 + n > pa) && (q3 + n < pb)

                # account for case where the percentile window falls in the upper portion of notch
                if inupper
                    push!(xsegs, r, r, l, l, L, R, r, r) # upper box
                    push!(ysegs, q3 + n, pb, pb, q3 + n, q3, q3, q3 + n, pb) # upper box
                    push!(texts,
                        "Q$ii: $pb",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Q$ii: $pb",
                        "Q$ii: $pb",
                        "Median: $q3 ± $n",
                        "",
                    )
                end

                # account for case where the percentile window falls in the lower portion of notch
                if inlower
                    push!(xsegs, r, r, R, L, l, l, r, r) # lower box
                    push!(ysegs, pa, q3 - n, q3, q3, q3 - n, pa, pa, q3 - n) # lower box
                    push!(texts, "Q$ii: $pa",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Median: $q3 ± $n",
                        "Q$ii: $pa",
                        "Q$ii: $pa",
                        "Median: $q3 ± $n",
                        "",
                    )
                end

                # account for case where the percentile window is completely outside the notch
                if !inupper & !inlower
                    push!(xsegs, r, r, l, l, r, r)
                    push!(ysegs, pa, pb, pb, pa, pa, pb)
                    push!(texts, "Q$ii: $pa", "Median: $pb", "Median: $pb", "Q$ii: $pa", "Q$ii: $pa", "Median: $pb", "")
                end
            end
        else
            # draw segments between requested percentile pairs
            for (ii, (pa, pb)) in enumerate(zip(ps, ps[2:end]))
                push!(xsegs, r, r, l, l, r, r)
                push!(ysegs, pa, pb, pb, pa, pa, pb)
                push!(texts, "Q$ii: $pa", "Median: $pb", "Median: $pb", "Q$ii: $pa", "Q$ii: $pa", "Median: $pb", "")
            end
        end

        push!(xsegs, m, lw, rw, m, m)             # upper T
        push!(ysegs, q5, q5, q5, q5, pmax)          # upper T
        push!(
            texts,
            "Upper fence: $q5",
            "Upper fence: $q5",
            "Upper fence: $q5",
            "Upper fence: $q5",
            "Q3: $pmax",
            "",
        )

    end

    @series begin
        # To prevent linecolor equal to fillcolor (It makes the median visible)
        if plotattributes[:linecolor] == plotattributes[:fillcolor]
            plotattributes[:linecolor] = plotattributes[:markerstrokecolor]
        end
        primary := true
        seriestype := :shape
        x := xsegs.pts
        y := ysegs.pts
        ()
    end

    # Outliers
    if outliers && !isempty(outliers)
        @series begin
            primary := false
            seriestype := :scatter
            if get!(plotattributes, :markershape, :circle) == :none
                plotattributes[:markershape] = :circle
            end

            fillrange := nothing
            x := outliers_x
            y := outliers_y
            ()
        end
    end

    # Hover
    primary := false
    seriestype := :path
    marker := false
    hover := texts
    linewidth := 0
    x := xsegs.pts
    y := ysegs.pts
    ()
end

Plots.@deps boxplot shape scatter

# ------------------------------------------------------------------------------
# Grouped Boxplot

@userplot GroupedBoxplot

recipetype(::Val{:groupedboxplot}, args...) = GroupedBoxplot(args)

@recipe function f(g::GroupedBoxplot; spacing = 0.1)
    x, y = grouped_xy(g.args...)

    # extract xnums and set default bar width.
    # might need to set xticks as well
    ux = unique(x)
    x = if eltype(x) <: Number
        bar_width --> (0.8 * mean(diff(sort(ux))))
        float.(x)
    else
        bar_width --> 0.8
        xnums = [findfirst(isequal(xi), ux) for xi in x] .- 0.5
        xticks --> (eachindex(ux) .- 0.5, ux)
        xnums
    end

    # shift x values for each group
    group = get(plotattributes, :group, nothing)
    if group != nothing
        gb = RecipesPipeline._extract_group_attributes(group)
        labels, idxs = getfield(gb, 1), getfield(gb, 2)
        n = length(labels)
        bws = plotattributes[:bar_width] / n
        bar_width := bws * clamp(1 - spacing, 0, 1)
        for i in 1:n
            groupinds = idxs[i]
            Δx = _cycle(bws, i) * (i - (n + 1) / 2)
            x[groupinds] .+= Δx
        end
    end

    seriestype := :boxplot
    x, y
end

Plots.@deps groupedboxplot boxplot