text:
################################################################################
# MOI wrapper
################################################################################
# This file is largely inspired from MOI wrappers existing in other
# JuliaOpt packages:
# https://www.juliaopt.org/
# The authors are indebted to the developers of JuliaOpt for
# the current MOI wrapper.
#
#--------------------------------------------------
# Specifications
#--------------------------------------------------
# The MOI wrapper works only for KNITRO version >= 11.0.
#
# - linear coefs are handled directly by KNITRO, via
#  `KN_add_con_linear_struct , KN_add_obj_linear_struct
#
# - quadratic coefs are also handled directly by KNITRO, via
#  `KN_add_con_quad_struct , KN_add_obj_quad_struct
#
# - NLP data are handled via KNITRO's callbacks
#
# NB: If `model.nlp_data` is empty, KNITRO would not use any
#     callbacks during solve.
#
#--------------------------------------------------

import MathOptInterface
const MOI  = MathOptInterface
const MOIU = MathOptInterface.Utilities

# TODO
const SF = Union{MOI.ScalarAffineFunction{Float64},
                 MOI.ScalarQuadraticFunction{Float64}}
const VAF = MOI.VectorAffineFunction{Float64}
const VOV = MOI.VectorOfVariables

const SS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64},
                 MOI.Interval{Float64}}
# LinSets
const LS = Union{MOI.EqualTo{Float64},
                 MOI.GreaterThan{Float64},
                 MOI.LessThan{Float64}}
# VecLinSets
const VLS = Union{MOI.Nonnegatives,
                  MOI.Nonpositives,
                  MOI.Zeros}

##################################################
# Define custom error for MOI wrapper.
struct UpdateObjectiveError <: MOI.NotAllowedError end
struct AddVariableError <: MOI.NotAllowedError end
struct AddConstraintError <: MOI.NotAllowedError end

# Import some utils.
include(joinpath("MOI_wrapper", "utils.jl"))

##################################################
mutable struct VariableInfo
    has_lower_bound::Bool # Implies lower_bound == Inf
    has_upper_bound::Bool # Implies upper_bound == Inf
    is_fixed::Bool        # Implies lower_bound == upper_bound and !has_lower_bound and !has_upper_bound.
    name::String
end
VariableInfo() = VariableInfo(false, false, false, "")

##################################################
# EmptyNLPEvaluator for non-NLP problems.
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    return
end
MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x) = nothing
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ)
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

##################################################
# MOI Optimizer
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{Model, Nothing}
    # We only keep in memory some information about variables
    # as we cannot delete variables, we do not have to store an index.
    variable_info::Vector{VariableInfo}
    # Get number of solve for restart.
    number_solved::Int
    # Specify if NLP is loaded inside KNITRO to avoid double definition.
    nlp_loaded::Bool
    nlp_data::MOI.NLPBlockData
    # Store index of nlp constraints.
    nlp_index_cons::Vector{Cint}
    # Store optimization sense.
    sense::MOI.OptimizationSense
    # Store the structure of the objective.
    objective::Union{MOI.SingleVariable,MOI.ScalarAffineFunction{Float64},MOI.ScalarQuadraticFunction{Float64},Nothing}
    # Constraint counters.
    number_zeroone_constraints::Int
    number_integer_constraints::Int
    # Constraint mappings.
    constraint_mapping::Dict{MOI.ConstraintIndex, Union{Cint, Vector{Cint}}}
    license_manager::Union{LMcontext, Nothing}
    options::Dict{String, Any}
end

function set_options(model::Optimizer, options)
    # Set KNITRO option.
    for (name, value) in options
        sname = string(name)
        MOI.set(model, MOI.RawParameter(sname), value)
    end
    return
end

function Optimizer(;license_manager=nothing, options...)
    # Create KNITRO context.
    if isa(license_manager, LMcontext)
        kc = Model(license_manager)
    else
        kc = Model()
    end
    # Convert Symbol to String in options dictionnary.
    options_dict = Dict{String, Any}()
    for (name, value) in options
        options_dict[string(name)] = value
    end
    model = Optimizer(kc, [], 0, false, empty_nlp_data(),
                      Cint[], MOI.FEASIBILITY_SENSE, nothing, 0, 0,
                      Dict{MOI.ConstraintIndex, Int}(), license_manager, Dict())

    set_options(model, options)
    return model
end

# Print Optimizer.
function Base.show(io::IO, model::Optimizer)
    println(io, "A MathOptInterface model with backend:")
    println(io, model.inner)
    return
end

# copy
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = true
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOI.Utilities.automatic_copy_to(model, src; kws...)
end

function MOI.empty!(model::Optimizer)
    # Free KNITRO model properly.
    if model.inner != nothing
        KN_free(model.inner)
    end
    # Handle properly license manager
    if isa(model.license_manager, LMcontext)
        model.inner = Model(model.license_manager)
    else
        model.inner = Model()
    end
    empty!(model.variable_info)
    model.number_solved = 0
    model.nlp_data = empty_nlp_data()
    model.nlp_loaded = false
    model.nlp_index_cons = Cint[]
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    model.number_zeroone_constraints = 0
    model.number_integer_constraints = 0
    model.constraint_mapping = Dict()
    model.license_manager = model.license_manager
    set_options(model, model.options)
    return
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
           model.nlp_data.evaluator isa EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           model.number_solved == 0 &&
           isa(model.objective, Nothing) &&
           model.number_zeroone_constraints == 0 &&
           model.number_integer_constraints == 0 &&
           !model.nlp_loaded
end

# Some utilities.
number_variables(model::Optimizer) = length(model.variable_info)
number_constraints(model::Optimizer) = KN_get_number_cons(model.inner)

# Getter for solver's name.
MOI.get(model::Optimizer, ::MOI.SolverName) = "Knitro"

# MOI.Silent.
MOI.supports(model::Optimizer, ::MOI.Silent) = true
function MOI.get(model::Optimizer, ::MOI.Silent)
    return KN_get_int_param(model.inner, "outlev") == 0
end

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    outlev = value ? 0 : 2
    KN_set_param(model.inner, "outlev", outlev)
    return
end

# MOI.TimeLimitSec.
MOI.supports(model::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return KN_get_double_param(model.inner, "maxtime_cpu")
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    # By default, maxtime is set to 1e8 in Knitro.
    maxtime = isnothing(value) ? 1e8 : value
    KN_set_param(model.inner, "maxtime_cpu", value)
    return
end

# MOI.RawParameters
function MOI.supports(model::Optimizer, param::MOI.RawParameter)
    name = param.name
    if name in KNITRO_OPTIONS || haskey(KN_paramName2Indx, name)
        return true
    end
    return false
end

function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    if !MOI.supports(model, p)
        throw(MOI.UnsupportedAttribute)
    end
    model.options[p.name] = value
    if p.name == "option_file"
        KN_load_param_file(model.inner, value)
    elseif p.name == "tuner_file"
        KN_load_tuner_file(model.inner, value)
    else
        KN_set_param(model.inner, p.name, value)
    end
    return
end

function MOI.get(model::Optimizer, p::MOI.RawParameter)
    if haskey(model.options, p.name)
        return model.options[p.name]
    end
    error("RawParameter with name $(p.name) is not set.")
end

##################################################
# Optimize
##################################################
function MOI.optimize!(model::Optimizer)
    KN_set_param(model.inner, "datacheck", 0)
    KN_set_param(model.inner, "hessian_no_f", 1)

    # Add NLP structure if specified.
    if !isa(model.nlp_data.evaluator, EmptyNLPEvaluator) && !model.nlp_loaded
        # Instantiate NLPEvaluator once and for all.
        features = MOI.features_available(model.nlp_data.evaluator)
        has_hessian = (:Hess in features)
        # Build initial features for solver.
        init_feat = Symbol[]
        model.nlp_data.has_objective && push!(init_feat, :Grad)
        has_hessian && push!(init_feat, :Hess)
        num_nlp_constraints = length(model.nlp_data.constraint_bounds)
        num_nlp_constraints > 0 && push!(init_feat, :Jac)
        # initialize!
        MOI.initialize(model.nlp_data.evaluator, init_feat)

        # Load NLP structure inside Knitro.
        offset = number_constraints(model)
        load_nlp_constraints(model)
        num_cons = KN_get_number_cons(model.inner)

        # The callbacks must match the signature of the callbacks
        # defined in knitro.h.
        # Objective callback (set both objective and constraint evaluation).
        function eval_f_cb(kc, cb, evalRequest, evalResult, userParams)
            # Evaluate objective if specified in nlp_data.
            if model.nlp_data.has_objective
                evalResult.obj[1] = MOI.eval_objective(model.nlp_data.evaluator,
                                                       evalRequest.x)
            end
            # Evaluate nonlinear term in constraint.
            MOI.eval_constraint(model.nlp_data.evaluator,
                                evalResult.c, evalRequest.x)
            return 0
        end

        # Objective gradient callback.
        function eval_grad_cb(kc, cb, evalRequest, evalResult, userParams)
            # Evaluate non-linear term in objective gradient.
            if model.nlp_data.has_objective
                MOI.eval_objective_gradient(model.nlp_data.evaluator,
                                            evalResult.objGrad,
                                            evalRequest.x)
            end
            # Evaluate non linear part of jacobian.
            MOI.eval_constraint_jacobian(model.nlp_data.evaluator,
                                         evalResult.jac,
                                         evalRequest.x)
            return 0
        end

        if has_hessian
            # Hessian callback.
            function eval_h_cb(kc, cb, evalRequest, evalResult, userParams)
                MOI.eval_hessian_lagrangian(model.nlp_data.evaluator,
                                            evalResult.hess,
                                            evalRequest.x,
                                            evalRequest.sigma,
                                            view(evalRequest.lambda, offset+1:num_cons))
                return 0
            end
        else
            eval_h_cb = nothing
        end

        # Be careful that sometimes objective is not evaluated here.
        # In any case, NLP objective has precedence over model.objective.
        if model.nlp_data.has_objective
            if num_nlp_constraints == 0
                # Add only a callback for objective if no NLP constraint
                cb = KN_add_objective_callback(model.inner, eval_f_cb)
            else
                cb = KN_add_eval_callback(model.inner, true, model.nlp_index_cons, eval_f_cb)
            end
        elseif num_nlp_constraints > 0
            cb = KN_add_eval_callback(model.inner, false, model.nlp_index_cons, eval_f_cb)

            # If a objective is specified in model.objective, load it.
            !isa(model.objective, Nothing) && add_objective!(model, model.objective)
        end

        # Gradient structure.
        nV = (model.nlp_data.has_objective) ? KN_DENSE : 0
        # Get jacobian structure.
        jacob_structure = MOI.jacobian_structure(model.nlp_data.evaluator)
        nnzJ = length(jacob_structure)
        if nnzJ == 0
            KN_set_cb_grad(model.inner, cb, eval_grad_cb, nV=nV)
        else
            # Take care to convert 1-indexing to 0-indexing!
            # KNITRO supports only Int32 array for integer.
            jacIndexVars = Int32[j - 1 for (_, j) in jacob_structure]
            # NLP constraints are set after all other constraints
            # inside Knitro.
            jacIndexCons = Int32[i - 1 + offset for (i, _) in jacob_structure]
            KN_set_cb_grad(model.inner, cb, eval_grad_cb, nV=nV,
                           jacIndexCons=jacIndexCons, jacIndexVars=jacIndexVars)
        end

        if has_hessian
            # Get hessian structure.
            hessian_structure = MOI.hessian_lagrangian_structure(model.nlp_data.evaluator)
            nnzH = length(hessian_structure)
            # Take care to convert 1-indexing to 0-indexing!
            # Knitro supports only Int32 array for integer.
            hessIndexVars1 = Int32[i - 1 for (i, _) in hessian_structure]
            hessIndexVars2 = Int32[j - 1 for (_, j) in hessian_structure]

            KN_set_cb_hess(model.inner, cb, nnzH, eval_h_cb,
                           hessIndexVars1=hessIndexVars1,
                           hessIndexVars2=hessIndexVars2)
        end

        model.nlp_loaded = true
    elseif !isa(model.objective, Nothing)
        add_objective!(model, model.objective)
    end

    KN_solve(model.inner)
    model.number_solved += 1
    return
end

include(joinpath("MOI_wrapper", "variables.jl"))
include(joinpath("MOI_wrapper", "constraints.jl"))
include(joinpath("MOI_wrapper", "objective.jl"))
include(joinpath("MOI_wrapper", "results.jl"))
include(joinpath("MOI_wrapper", "nlp.jl"))


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

mutable struct IoK8sApiCoreV1Capabilities <: SwaggerModel
    add::Any # spec type: Union{ Nothing, Vector{String} } # spec name: add
    drop::Any # spec type: Union{ Nothing, Vector{String} } # spec name: drop

    function IoK8sApiCoreV1Capabilities(;add=nothing, drop=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Capabilities, Symbol("add"), add)
        setfield!(o, Symbol("add"), add)
        validate_property(IoK8sApiCoreV1Capabilities, Symbol("drop"), drop)
        setfield!(o, Symbol("drop"), drop)
        o
    end
end # type IoK8sApiCoreV1Capabilities

const _property_map_IoK8sApiCoreV1Capabilities = Dict{Symbol,Symbol}(Symbol("add")=>Symbol("add"), Symbol("drop")=>Symbol("drop"))
const _property_types_IoK8sApiCoreV1Capabilities = Dict{Symbol,String}(Symbol("add")=>"Vector{String}", Symbol("drop")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1Capabilities }) = collect(keys(_property_map_IoK8sApiCoreV1Capabilities))
Swagger.property_type(::Type{ IoK8sApiCoreV1Capabilities }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1Capabilities[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Capabilities }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Capabilities[property_name]

function check_required(o::IoK8sApiCoreV1Capabilities)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Capabilities }, name::Symbol, val)
end

ext: jl


text:
module PlotNested

using Requires,
      FieldMetadata,
      Interact,
      Plots

import FieldMetadata: @plottable, plottable

export @plottable, plottable, plottables, plotnames, plotdata, plotchecks, plot_selected, plot_all, autoplot

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        plottables(x::DataFrames.DataFrame, fname) = (fname => x,)
    end
    @require TypedTables="9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" begin
        plottables(x::TypedTables.Table, fname) = (fname => x,)
    end
end

##### Flattening

plottables_expr(T, index) = quote
    if plottable($T, Val{$(QuoteNode(index))})
        plottables(getfield(t, $(QuoteNode(index))), $(QuoteNode(index)))
    else
        ()
    end
end

plottables_expr(T::Type{X}, index) where X <: Tuple = quote
    plottables(getfield(t, $(QuoteNode(index))), Symbol(fname, $(QuoteNode(index))))
end

wrap(expressions) = Expr(:tuple, Expr.(:..., expressions)...)

plottables(x) = plottables(x, :unnamed)
plottables(x::AbstractArray, fname) = (fname => x,)
plottables(x::Tuple{N,AbstractArray}, fname) where N = (fname => x,)
plottables(x::Tuple{}, fname) where N = ()
@generated plottables(t::T, fname) where T =
    wrap([plottables_expr(T, fn) for fn in fieldnames(T)])

##### Plot generation
plotnames(x) = getfield.(plottables(x), 1)

plotdata(x) = getfield.(plottables(x), 2)
plotdata(x, range) = [getrange(d, range) for d in plotdata(x)]
getrange(d::AbstractVector, range) = begin
    d[first(range):min(last(range), length(d))]
end
getrange(d::AbstractMatrix, range) = d[first(range):min(last(range), length(d)), :]
getrange(d::AbstractArray{T,3}, range) where T = d[:, :, first(range):min(last(range), length(d))]

plotchecks(x) = [checkbox(false, label=string(name)) for name in plotnames(x)]

plot_selected(x, checklist, args...) = begin
    labels = map(string, plotnames(x))
    d = plotdata(x, args...)
    [plot(d[i], label=labels[i]) for (i, check) in enumerate(checklist) if check]
end
plot_all(x) = [plot(data) for (i, (name, data)) in enumerate(plottables(x))]

autoplot(x) = plot(plot_all(x)...)

end # module


text:
using DelimitedFiles
cd(dirname(@__FILE__))
X=readdlm("X.txt",',',Int16)
genelist=vec(readdlm("genelist.txt",String))

using Statistics
glogmean=log10.(mean(X,dims=2));
glogvar=log10.(var(X,dims=2));
dropoutrate=mean(X.==0,dims=2);
using SparseArrays
nnz(sparse(X))
count(!iszero,X)

using Plots
Plots.scatter(glogmean,glogvar)
Plots.scatter(glogmean,1 .-dropoutrate)

using Distributions
# fit_mle(NegativeBinomial,X[3,:].+0.0)

fit_mle(Normal,X[3,:].+0.0)
# fit(Normal,X[3,:].+0.0)
# using GLM

# http://naobioml.blogspot.com/2017/02/how-to-fit-count-data-with-negative.html

data=X[3,:].+0.0;
function f2(x::Vector)
    sum = 0
    for i in data
        sum += log(pdf(NegativeBinomial(x[1], x[2]), i))
    end
    return -1*sum
end

using Optim, Distributions

lower = [0.0, 0.0]
upper = [Inf, 1]
initial_x = [0.5, 0.5]
# x2 = optimize(DifferentiableFunction(f2), initial_x, lower, upper, Fminbox(), optimizer = GradientDescent)
x2 = optimize(f2, [0.0, 0.0], [Inf, 1],[0.5, 0.5])

function f1(x)
    sum = 0
    for i in data
        sum += log(pdf(Poisson(x), i))
    end
    -1*sum
end
x1 = optimize(f1, 0.0, 7)

#=
 ------------

using Distributions, Optim, StatsPlots
julia> data = vcat([0 for i in 1:70],
           [1 for i in 1:38],
           [2 for i in 1:17],
           [3 for i in 1:10],
           [4 for i in 1:9],
           [5 for i in 1:3],
           [6 for i in 1:2],
           [7 for i in 1:1]);
julia> function f2(x)
           sum = 0.0
           for i in data
               sum += log(pdf(NegativeBinomial(x[1], x[2]), i))
           end
           return -sum
       end
f2 (generic function with 1 method)
julia> opt_result = optimize(f2, [0.0, 0.0], [Inf, 1],[0.5, 0.5])
 * Status: success
 * Candidate solution
    Minimizer: [1.02e+00, 4.72e-01]
    Minimum:   2.224372e+02
 * Found with
    Algorithm:     Fminbox with L-BFGS
    Initial Point: [5.00e-01, 5.00e-01]
 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 7.79e-08 ≰ 1.0e-08
 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    4
    f(x) calls:    226
    ∇f(x) calls:   226
julia> histogram(data, normalize = true, label = "Data", alpha = 0.5, linecolor = "white"); plot!(NegativeBinomial(opt_result.minimizer[1], opt_result.minimizer[2]), label = "Negative Binomial fit")

=#


text:
push!(LOAD_PATH,"../src/")

using Documenter, MissingsAsFalse

makedocs(
    sitename = "MissingsAsFalse.jl",
    pages = Any[
        "Introduction" => "index.md",
        "API" => "api/api.md"],
    format = Documenter.HTML(
        canonical = "https://pdeffebach.github.io/MissingsAsFalse.jl/stable/"
    ))

deploydocs(
    repo = "github.com/pdeffebach/MissingsAsFalse.jl.git",
    target = "build",
    deps = nothing,
    make = nothing)


text:
# Basic algebra
#---------------
Base.copy(t::AbstractTensorMap) = Base.copyto!(similar(t), t)

Base.:-(t::AbstractTensorMap) = mul!(similar(t), t, -one(eltype(t)))
function Base.:+(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(one(T), t2, copyto!(similar(t1, T), t1))
end
function Base.:-(t1::AbstractTensorMap, t2::AbstractTensorMap)
    T = promote_type(eltype(t1), eltype(t2))
    return axpy!(-one(T), t2, copyto!(similar(t1, T), t1))
end

Base.:*(t::AbstractTensorMap, α::Number) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), t, α)
Base.:*(α::Number, t::AbstractTensorMap) =
    mul!(similar(t, promote_type(eltype(t), typeof(α))), α, t)
Base.:/(t::AbstractTensorMap, α::Number) = *(t, one(α)/α)
Base.:\(α::Number, t::AbstractTensorMap) = *(t, one(α)/α)

LinearAlgebra.normalize!(t::AbstractTensorMap, p::Real = 2) = rmul!(t, inv(norm(t, p)))
LinearAlgebra.normalize(t::AbstractTensorMap, p::Real = 2) =
    mul!(similar(t), t, inv(norm(t, p)))

Base.:*(t1::AbstractTensorMap, t2::AbstractTensorMap) =
    mul!(similar(t1, promote_type(eltype(t1),eltype(t2)), codomain(t1)←domain(t2)), t1, t2)
Base.exp(t::AbstractTensorMap) = exp!(copy(t))
Base.:^(t::AbstractTensorMap, p::Integer) =
    p < 0 ? Base.power_by_squaring(inv(t), -p) : Base.power_by_squaring(t, p)

# Special purpose constructors
#------------------------------
Base.zero(t::AbstractTensorMap) = fill!(similar(t), 0)
function Base.one(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SectorMismatch("no identity if domain and codomain are different"))
    return one!(similar(t))
end
function one!(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SectorMismatch("no identity if domain and codomain are different"))
    for (c,b) in blocks(t)
        _one!(b)
    end
    return t
end

"""
    id([A::Type{<:DenseMatrix} = Matrix{Float64},] space::VectorSpace) -> TensorMap

Construct the identity endomorphism on space `space`, i.e. return a `t::TensorMap` with `domain(t) == codomain(t) == V`, where `storagetype(t) = A` can be specified.
"""
id(A, V::ElementarySpace) = id(A, ProductSpace(V))
id(V::VectorSpace) = id(Matrix{Float64}, V)
id(::Type{A}, P::ProductSpace) where {A<:DenseMatrix} =
    one!(TensorMap(s->A(undef, s), P, P))

"""
    isomorphism([A::Type{<:DenseMatrix} = Matrix{Float64},]
                    cod::VectorSpace, dom::VectorSpace)
    -> TensorMap

Return a `t::TensorMap` that implements a specific isomorphism between the codomain `cod`
and the domain `dom`, and for which `storagetype(t)` can optionally be chosen to be of type
`A`. If the two spaces do not allow for such an isomorphism, and are thus not isomorphic,
and error will be thrown. When they are isomorphic, there is no canonical choice for a
specific isomorphism, but the current choice is such that
`isomorphism(cod, dom) == inv(isomorphism(dom, cod))`.

See also [`unitary`](@ref) when `spacetype(cod) isa EuclideanSpace`.
"""
isomorphism(cod::TensorSpace, dom::TensorSpace) = isomorphism(Matrix{Float64}, cod, dom)
isomorphism(P::TensorMapSpace) = isomorphism(codomain(P), domain(P))
isomorphism(A::Type{<:DenseMatrix}, P::TensorMapSpace) =
    isomorphism(A, codomain(P), domain(P))
isomorphism(A::Type{<:DenseMatrix}, cod::TensorSpace, dom::TensorSpace) =
    isomorphism(A, convert(ProductSpace, cod), convert(ProductSpace, dom))
function isomorphism(::Type{A}, cod::ProductSpace, dom::ProductSpace) where {A<:DenseMatrix}
    spacetype(cod) == spacetype(dom) ||
        throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic"))
    for c in union(blocksectors(cod), blocksectors(dom))
        blockdim(cod, c) == blockdim(dom, c) ||
            throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic"))
    end
    t = TensorMap(s->A(undef, s), cod, dom)
    for (c,b) in blocks(t)
        _one!(b)
    end
    return t
end

const EuclideanTensorSpace = TensorSpace{<:EuclideanSpace}
const EuclideanTensorMapSpace = TensorMapSpace{<:EuclideanSpace}
const AbstractEuclideanTensorMap = AbstractTensorMap{<:EuclideanTensorSpace}
const EuclideanTensorMap = TensorMap{<:EuclideanTensorSpace}

"""
    unitary([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace)
    -> TensorMap

Return a `t::TensorMap` that implements a specific unitary isomorphism between the codomain
`cod` and the domain `dom`, for which `spacetype(dom)` (`== spacetype(cod)`) must be a
subtype of `EuclideanSpace`. Furthermore, `storagetype(t)` can optionally be chosen to be
of type `A`. If the two spaces do not allow for such an isomorphism, and are thus not
isomorphic, and error will be thrown. When they are isomorphic, there is no canonical choice
for a specific isomorphism, but the current choice is such that
`isomorphism(cod, dom) == inv(isomorphism(dom, cod)) = adjoint(unitary(dom, cod))`.
"""
unitary(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) = isomorphism(cod, dom)
unitary(P::EuclideanTensorMapSpace) = isomorphism(P)
unitary(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) = isomorphism(A, P)
unitary(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isomorphism(A, cod, dom)

"""
    isometry([A::Type{<:DenseMatrix} = Matrix{Float64},] cod::VectorSpace, dom::VectorSpace)
    -> TensorMap

Return a `t::TensorMap` that implements a specific isometry that embeds the domain `dom`
into the codomain `cod`, and which requires that `spacetype(dom)` (`== spacetype(cod)`) is
a subtype of `EuclideanSpace`. An isometry `t` is such that its adjoint `t'` is the left
inverse of `t`, i.e. `t'*t = id(dom)`, while `t*t'` is some idempotent endomorphism of
`cod`, i.e. it squares to itself. When `dom` and `cod` do not allow for such an isometric
inclusion, an error will be thrown.
"""
isometry(cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isometry(Matrix{Float64}, cod, dom)
isometry(P::EuclideanTensorMapSpace) = isometry(codomain(P), domain(P))
isometry(A::Type{<:DenseMatrix}, P::EuclideanTensorMapSpace) =
    isometry(A, codomain(P), domain(P))
isometry(A::Type{<:DenseMatrix}, cod::EuclideanTensorSpace, dom::EuclideanTensorSpace) =
    isometry(A, convert(ProductSpace, cod), convert(ProductSpace, dom))
function isometry(::Type{A},
                    cod::ProductSpace{S},
                    dom::ProductSpace{S}) where {A<:DenseMatrix, S<:EuclideanSpace}
    spacetype(cod) == spacetype(dom) || throw(SpaceMismatch())
    for c in union(blocksectors(cod), blocksectors(dom))
        blockdim(cod, c) >= blockdim(dom, c) ||
            throw(SpaceMismatch("codomain $cod and domain $dom do not allow for an isometric mapping"))
    end
    t = TensorMap(s->A(undef, s), cod, dom)
    for (c,b) in blocks(t)
        _one!(b)
    end
    return t
end

# In-place methods
#------------------
# Wrapping the blocks in a StridedView enables multithreading if JULIA_NUM_THREADS > 1
# Copy, adjoint! and fill:
function Base.copyto!(tdest::AbstractTensorMap, tsource::AbstractTensorMap)
    codomain(tdest) == codomain(tsource) && domain(tdest) == domain(tsource) ||
        throw(SpaceMismatch())
    for c in blocksectors(tdest)
        copyto!(StridedView(block(tdest, c)), StridedView(block(tsource, c)))
    end
    return tdest
end
function Base.fill!(t::AbstractTensorMap, value::Number)
    for (c,b) in blocks(t)
        fill!(b, value)
    end
    return t
end
function LinearAlgebra.adjoint!(tdst::AbstractEuclideanTensorMap,
                                tsrc::AbstractEuclideanTensorMap)
    (codomain(tdst) == domain(tsrc) && domain(tdst) == codomain(tsrc)) ||
        throw(SpaceMismatch())
    for c in blocksectors(tdst)
        adjoint!(StridedView(block(tdst, c)), StridedView(block(tsrc, c)))
    end
    return tdst
end

# Basic vector space methods: addition and scalar multiplication
LinearAlgebra.rmul!(t::AbstractTensorMap, α::Number) = mul!(t, t, α)
LinearAlgebra.lmul!(α::Number, t::AbstractTensorMap) = mul!(t, α, t)

function LinearAlgebra.mul!(t1::AbstractTensorMap, t2::AbstractTensorMap, α::Number)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), StridedView(block(t2, c)), α)
    end
    return t1
end
function LinearAlgebra.mul!(t1::AbstractTensorMap, α::Number, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        mul!(StridedView(block(t1, c)), α, StridedView(block(t2, c)))
    end
    return t1
end
function LinearAlgebra.axpy!(α::Number, t1::AbstractTensorMap, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpy!(α, StridedView(block(t1, c)), StridedView(block(t2, c)))
    end
    return t2
end
function LinearAlgebra.axpby!(α::Number, t1::AbstractTensorMap,
                                β::Number, t2::AbstractTensorMap)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    for c in blocksectors(t1)
        axpby!(α, StridedView(block(t1, c)), β, StridedView(block(t2, c)))
    end
    return t2
end

# inner product and norm only valid for spaces with Euclidean inner product
function LinearAlgebra.dot(t1::AbstractEuclideanTensorMap, t2::AbstractEuclideanTensorMap)
    (codomain(t1) == codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceMismatch())
    iter = blocksectors(t1)
    if isempty(iter)
        return zero(eltype(t1))*zero(eltype(t2))
    else
        return sum(dim(c)*dot(block(t1,c), block(t2,c)) for c in blocksectors(t1))
    end
end

LinearAlgebra.norm(t::AbstractEuclideanTensorMap, p::Real) =
    isempty(blocks(t)) ? zero(real(eltype(t))) : _norm(blocks(t), p)
function _norm(blockiterator, p::Real)
    if p == Inf
        maximum(norm(b, p) for (c,b) in blockiterator)
    elseif p == -Inf
        minimum(norm(b, p) for (c,b) in blockiterator)
    elseif p == 1
        sum(dim(c)*norm(b, p) for (c,b) in blockiterator)
    else
        s = sum(dim(c)*norm(b, p)^p for (c,b) in blockiterator)
        return s^inv(oftype(s, p))
    end
end

# TensorMap trace
function LinearAlgebra.tr(t::AbstractTensorMap)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("Trace of a tensor only exist when domain == codomain"))
    return sum(dim(c)*tr(b) for (c,b) in blocks(t))
end

# TensorMap multiplication:
function LinearAlgebra.mul!(tC::AbstractTensorMap,
                            tA::AbstractTensorMap,
                            tB::AbstractTensorMap, α = true, β = false)
    if !(codomain(tC) == codomain(tA) && domain(tC) == domain(tB) && domain(tA) == codomain(tB))
        throw(SpaceMismatch())
    end
    for c in blocksectors(tC)
        if hasblock(tA, c) # then also tB should have such a block
            A = block(tA, c)
            B = block(tB, c)
            C = block(tC, c)
            if isbitstype(eltype(A)) && isbitstype(eltype(B)) && isbitstype(eltype(C))
                @unsafe_strided A B C mul!(C, A, B, α, β)
            else
                mul!(StridedView(C), StridedView(A), StridedView(B), α, β)
            end
        elseif β != one(β)
            rmul!(block(tC, c), β)
        end
    end
    return tC
end

# TensorMap inverse
function Base.inv(t::AbstractTensorMap)
    cod = codomain(t)
    dom = domain(t)
    for c in union(blocksectors(cod), blocksectors(dom))
        blockdim(cod, c) == blockdim(dom, c) ||
            throw(SpaceMismatch("codomain $cod and domain $dom are not isomorphic: no inverse"))
    end
    if sectortype(t) === Trivial
        return TensorMap(inv(block(t, Trivial())), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c,b) in blocks(t)
            data[c] = inv(b)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end
function LinearAlgebra.pinv(t::AbstractTensorMap; kwargs...)
    if sectortype(t) === Trivial
        return TensorMap(pinv(block(t, Trivial()); kwargs...), domain(t)←codomain(t))
    else
        data = empty(t.data)
        for (c,b) in blocks(t)
            data[c] = pinv(b; kwargs...)
        end
        return TensorMap(data, domain(t)←codomain(t))
    end
end
function Base.:(\)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    codomain(t1) == codomain(t2) ||
        throw(SpaceMismatch("non-matching codomains in t1 \\ t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) \ block(t2, Trivial())
        return TensorMap(data, domain(t1)←domain(t2))
    else
        cod = codomain(t1)
        data = SectorDict(c=>block(t1,c) \ block(t2,c) for c in blocksectors(codomain(t1)))
        return TensorMap(data, domain(t1)←domain(t2))
    end
end
function Base.:(/)(t1::AbstractTensorMap, t2::AbstractTensorMap)
    domain(t1) == domain(t2) ||
        throw(SpaceMismatch("non-matching domains in t1 / t2"))
    if sectortype(t1) === Trivial
        data = block(t1, Trivial()) / block(t2, Trivial())
        return TensorMap(data, codomain(t1)←codomain(t2))
    else
        data = SectorDict(c=>block(t1,c) / block(t2,c) for c in blocksectors(domain(t1)))
        return TensorMap(data, codomain(t1)←codomain(t2))
    end
end

# TensorMap exponentation:
function exp!(t::TensorMap)
    domain(t) == codomain(t) ||
        error("Exponentional of a tensor only exist when domain == codomain.")
    for (c,b) in blocks(t)
        copyto!(b, LinearAlgebra.exp!(b))
    end
    return t
end

# functions that map ℝ to (a subset of) ℝ
for f in (:cos, :sin, :tan, :cot, :cosh, :sinh, :tanh, :coth, :atan, :acot, :asinh)
    sf = string(f)
    @eval function Base.$f(t::AbstractTensorMap)
        domain(t) == codomain(t) ||
            error("$sf of a tensor only exist when domain == codomain.")
        G = sectortype(t)
        T = similarstoragetype(t, float(eltype(t)))
        if sectortype(t) === Trivial
            local data::T
            if eltype(t) <: Real
                data = real($f(block(t, Trivial())))
            else
                data = $f(block(t, Trivial()))
            end
            return TensorMap(data, codomain(t), domain(t))
        else
            if eltype(t) <: Real
                datadict = SectorDict{G,T}(c=>real($f(b)) for (c,b) in blocks(t))
            else
                datadict = SectorDict{G,T}(c=>$f(b) for (c,b) in blocks(t))
            end
            return TensorMap(datadict, codomain(t), domain(t))
        end
    end
end
# functions that don't map ℝ to (a subset of) ℝ
for f in (:sqrt, :log, :asin, :acos, :acosh, :atanh, :acoth)
    sf = string(f)
    @eval function Base.$f(t::AbstractTensorMap)
        domain(t) == codomain(t) ||
            error("$sf of a tensor only exist when domain == codomain.")
        G = sectortype(t)
        T = similarstoragetype(t, complex(float(eltype(t))))
        if sectortype(t) === Trivial
            data::T = $f(block(t, Trivial()))
            return TensorMap(data, codomain(t), domain(t))
        else
            datadict = SectorDict{G,T}(c=>$f(b) for (c,b) in blocks(t))
            return TensorMap(datadict, codomain(t), domain(t))
        end
    end
end

# # concatenate tensors
# function concatenate(t1::AbstractTensorMap{S}, ts::AbstractTensorMap{S}...; direction::Symbol) where {S}
#     if direction = :domain
#         numin(t1) == 1 || throw(SpaceMismatch("concatenation along domain"))
#         for t2 in ts
#             domain(t1) == domain(t2) || throw(SpaceMismatch())
#         end
#         t = t1
#         for t2 in ts
#             cat(t1)

function catdomain(t1::AbstractTensorMap{S,N₁,1}, t2::AbstractTensorMap{S,N₁,1}) where {S,N₁}
    codomain(t1) == codomain(t2) || throw(SpaceMismatch())

    V1, = domain(t1)
    V2, = domain(t2)
    isdual(V1) == isdual(V2) ||
        throw(SpaceMismatch("cannot hcat tensors whose domain has non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1), eltype(t2)), codomain(t1), V)
    for c in sectors(V)
        block(t, c)[:,1:dim(V1, c)] .= block(t1, c)
        block(t, c)[:,dim(V1, c) .+ (1:dim(V2, c))] .= block(t2, c)
    end
    return t
end
function catcodomain(t1::AbstractTensorMap{S,1,N₂}, t2::AbstractTensorMap{S,1,N₂}) where {S,N₂}
    domain(t1) == domain(t2) || throw(SpaceMismatch())

    V1, = codomain(t1)
    V2, = codomain(t2)
    isdual(V1) == isdual(V2) ||
        throw(SpaceMismatch("cannot hcat tensors whose domain has non-matching duality"))

    V = V1 ⊕ V2
    t = TensorMap(undef, promote_type(eltype(t1),eltype(t2)), V, domain(t1))
    for c in sectors(V)
        block(t, c)[1:dim(V1, c),:] .= block(t1, c)
        block(t, c)[dim(V1,c) .+ (1:dim(V2, c)),:] .= block(t2, c)
    end
    return t
end

# tensor product of tensors
"""
    ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}, ...) -> TensorMap{S}

Compute the tensor product between two `AbstractTensorMap` instances, which results in a
new `TensorMap` instance whose codomain is `codomain(t1) ⊗ codomain(t2)` and whose domain
is `domain(t1) ⊗ domain(t2)`.
"""
function ⊗(t1::AbstractTensorMap{S}, t2::AbstractTensorMap{S}) where S
    cod1, cod2 = codomain(t1), codomain(t2)
    dom1, dom2 = domain(t1), domain(t2)
    cod = cod1 ⊗ cod2
    dom = dom1 ⊗ dom2
    t = TensorMap(zeros, promote_type(eltype(t1), eltype(t2)), cod, dom)
    if sectortype(S) === Trivial
        d1 = dim(cod1)
        d2 = dim(cod2)
        d3 = dim(dom1)
        d4 = dim(dom2)
        m1 = reshape(t1[], (d1, 1, d3, 1))
        m2 = reshape(t2[], (1, d2, 1, d4))
        m = reshape(t[], (d1, d2, d3, d4))
        m .= m1 .* m2
    else
        for (f1l, f1r) in fusiontrees(t1)
            for (f2l, f2r) in fusiontrees(t2)
                c1 = f1l.coupled # = f1r.coupled
                c2 = f2l.coupled # = f2r.coupled
                for c in c1 ⊗ c2
                    degeneracyiter = FusionStyle(c) isa DegenerateNonAbelian ?
                                        (1:Nsymbol(c1, c2, c)) : (nothing,)
                    for μ in degeneracyiter
                        for (fl, coeff1) in merge(f1l, f2l, c, μ)
                            for (fr, coeff2) in merge(f1r, f2r, c, μ)
                                d1 = dim(cod1, f1l.uncoupled)
                                d2 = dim(cod2, f2l.uncoupled)
                                d3 = dim(dom1, f1r.uncoupled)
                                d4 = dim(dom2, f2r.uncoupled)
                                m1 = reshape(t1[f1l,f1r], (d1, 1, d3, 1))
                                m2 = reshape(t2[f2l,f2r], (1, d2, 1, d4))
                                m = reshape(t[fl, fr], (d1, d2, d3, d4))
                                m .+= coeff1 .* coeff2 .* m1 .* m2
                            end
                        end
                    end
                end
            end
        end
    end
    return t
end


text:
@doc raw"""
    update!(filter, y, u)

Update the internal states of the filter in accordance with the output y
and input u.

# Methods
    update!(filter::KalmanFilter, y, u)

Update the internal states of the KalmanFilter in according to
measured output $y$ and controlled input $u$.

The internal states are updated as follows

```math
\begin{aligned}
P_{t+1} & = Q + FP_{t}F^\top \quad - FP_{t}H^\top(R + HP_{t}H^\top)^{-1}H P_{t}F^\top \\
\breve{x}_{t+1}   & = F \breve{x}_{t} + K_{t}(y_t - H\breve{x}_{t})\\
K_{t}     & = FP_{t}H^\top(R + H P_{t} H^\top)^{-1} \\
c_{t+1} & = |H_i\breve x_{t} -y_t|^2_{(R + HP_{t}H^\top)^{-1}} + c_{t} \\
p(i|y_t) & = \frac{1}{(2\pi)^{n/2}\sqrt{|R + HP_tH^\top|}}e^{-1/2 (H\hat x_t - y_t)^\top (R + HP_tH^\top)^{-1}(Hx_t-y_t)}p(i|y_{t-1}).
\end{aligned}
```

The cummulative costs $c_t$ are states used by Minimax Adaptive Estimators, and $p(i|y_t)$ are states used by the standard adaptive estimator.

    update!(filter::MinimaxFilter, y, u)

Update each filter in the internal filterbank of the minimaxfilter object.

    update!(filter::BayesianFilter, y, u)

Update each filter in the internal filterbank of the standard multiple model adaptive filter
and normalize the probability of a given state being active conditioned on past measurements,
$p(i, y_t)$.
"""
function update!(filter::AbstractFilter, y, u)
end

function update!(
    filter::KalmanFilter{T},
    y::AbstractVector{T},
    u::AbstractVector{T},
    ) where {T<:Number}

    x = filter.x
    F = filter.F
    H = filter.H
    B = filter.B
    P = filter.P
    n = length(x)
    offset = filter.offset
    Re = (filter.R + H * P * H')
    K = F * P * H' / Re

    filter.c[] += ((H * x - y)' / Re) * (H * x - y)
    filter.pdy[] *= 1 / (2*π)^(n/2) / sqrt(det(Re)) * exp(-1/2 * ((H*x -y)' / Re)*(H*x - y))
    filter.x[:] = F*x + B*u + K*(y-H*x) + offset
    filter.P[:,:] = filter.Q + F*P*F' - F*P*(H' / Re) * H *P*F'
    return nothing
end

function update!(
    filter::MinimaxMMAE{T},
    y::AbstractVector{T},
    u::AbstractVector{T},
) where {T<:Number}
    for f in filter.filterbank
        update!(f, y, u)
    end

    return nothing
end

function update!(
    filter::BayesianMMAE{T},
    y::AbstractVector{T},
    u::AbstractVector{T},
) where {T<:Number}
    sum_pdys = 0
    for f in filter.filterbank
        update!(f, y, u)
        sum_pdys += f.pdy[]
    end

    for f in filter.filterbank
        f.pdy[] = f.pdy[]/sum_pdys
    end

    return nothing
end

@doc raw"""
    predict(filter)

Predict the states at the next time instance.

# Methods

    xhat = predict(filter::KalmanFilter)

Get the kalman filter prediction of the state at the next time-step, $\hat x_{t+1}$.

    yhat, val = predict(minimaxfilter)

Predict the output of the next time-step as the minimizing argument of
 the quadratically constrainde convex program
```math
J_N^\star(y)= \min_{\hat y} \max_i |\hat y - H_i\breve x_{N,i}|^2_{(I-\gamma^{-2}H_iP_{N,i}H_i^\top)^{-1}} - \gamma^2 c_{N,i}.
```

    yhat = predict(bayesianfilter)

Predict the output at the next time-step as the expected value
`` E[\hat y_{t+1}] = \sum_{i=1}^K \breve y_i p(i|y_t) ``, where $\breve y_i$ are the
corresponding KalmanFilter estimates.
"""
function predict(filter::AbstractFilter)
end

function predict(filter::KalmanFilter)
    return copy(filter.x)
end

function predict(filter::MinimaxMMAE, γ::Real; silent = false)
    γ <= 0 ? DomainError("γ must be positive") : nothing
    m = size(filter.filterbank[1].H)[1]
    K = length(filter.filterbank)

    # Setting up optimization model
    model = Model(filter.optimizer)
    silent ? set_silent(model) : unset_silent(model)
    @variable(model, yhat[1:m])
    @variable(model, t)
    for f in filter.filterbank
        W = (I - γ^(-2)*f.H*f.P*f.H')
        @constraint(
            model,
            (yhat - f.H*f.x)'*pinv(W)*(yhat - f.H*f.x) - γ^2*f.c[] <= t
        )
    end

    @objective(model, Min, t)
    optimize!(model)
    yhat = value.(yhat)
    val = value(t)
    return yhat, val
end

function predict(filter::MinimaxMMAE; γmax = 100, silent = true)
    γmin = sqrt(maximum([maximum(eigvals(kf.H*kf.P*kf.H')) for kf in filter.filterbank]))
    γ = (γmin + γmax)/2
    yhat = []
    val = 0
    for k = 1:20
        yhat, val = predict(filter, γ, silent = silent)
        val < 0 ? γmax = γ : γmin = γ
        γ = (γmin + γmax)/2
    end
    return yhat, val, γ
end

function predict(filter::BayesianMMAE)
    yhat = [f.pdy[]*f.H*f.x for f in filter.filterbank] |> sum
    return yhat
end

# Lazy wrappers

function update!(filter::AbstractFilter{T}, y::Number) where {T<:Number}
    return update!(filter, [T(y)], 0)
end

function update!(filter::AbstractFilter{T}, y::AbstractVector) where {T<:Number}
    return update!(filter, T.(y), 0)
end

function update!(filter::AbstractFilter{T}, y::Number, u::Number) where {T<:Number}
    return update!(filter, [T(y)], [T(u)])
end

function update!(filter::AbstractFilter{T}, y::Number, u::AbstractVector) where {T<:Number}
    return update!(filter, [T(y)], T.(u))
end

function update!(filter::AbstractFilter{T}, y::AbstractVector, u::Number) where {T<:Number}
    return update!(filter, T.(y), [T.(u)])
end

function update!(
    filter::AbstractFilter{T},
    y::AbstractVector,
    u::AbstractVector,
) where {T<:Number}

    return update!(filter, T.(y), T.(u))
end


text:
function setup_server(env = dirname(SymbolServer.Pkg.Types.Context().env.project_file), depot = first(SymbolServer.Pkg.depots()), cache = joinpath(dirname(pathof(SymbolServer)), "..", "store"))
    server = StaticLint.FileServer()
    ssi = SymbolServerInstance(depot, cache)
    _, server.symbolserver = SymbolServer.getstore(ssi, env)
    server.symbol_extends  = SymbolServer.collect_extended_methods(server.symbolserver)
    server
end

"""
    lint_string(s, server; gethints = false)

Parse a string and run a semantic pass over it. This will mark scopes, bindings,
references, and lint hints. An annotated `EXPR` is returned or, if `gethints = true`,
it is paired with a collected list of errors/hints.
"""
function lint_string(s::String, server = setup_server(); gethints = false)
    empty!(server.files)
    f = File("", s, CSTParser.parse(s, true), nothing, server)
    env = getenv(f, server)
    setroot(f, f)
    setfile(server, "", f)
    semantic_pass(f)
    check_all(f.cst, LintOptions(), env)
    if gethints
        return f.cst, [(x, string(haserror(x) ? LintCodeDescriptions[x.meta.error] : "Missing reference", " at offset ", offset)) for (offset, x) in collect_hints(f.cst, env)]
    else
        return f.cst
    end
end

"""
    lint_file(rootpath, server)

Read a file from disc, parse and run a semantic pass over it. The file should be the
root of a project, e.g. for this package that file is `src/StaticLint.jl`. Other files
in the project will be loaded automatically (calls to `include` with complicated arguments
are not handled, see `followinclude` for details). A `FileServer` will be returned
containing the `File`s of the package.
"""
function lint_file(rootpath, server = setup_server(); gethints = false)
    empty!(server.files)
    root = loadfile(server, rootpath)
    semantic_pass(root)
    for f in values(server.files)
        check_all(f.cst, LintOptions(), getenv(f, server))
    end
    if gethints
        hints = []
        for (p,f) in server.files
            append!(hints, [(x, string(haserror(x) ? LintCodeDescriptions[x.meta.error] : "Missing reference", " at offset ", offset, " of ", p)) for (offset, x) in collect_hints(f.cst, getenv(f, server))])
        end
        return root, hints
    else
        return root
    end
end


text:
import MultivariatePolynomials
const MP = MultivariatePolynomials

import MultivariateBases
const MB = MultivariateBases

@testset "_merge_sorted" begin
    @test SumOfSquares.Certificate._merge_sorted([4, 1], [3, 0]) == [4, 3, 1, 0]
    @test SumOfSquares.Certificate._merge_sorted((4, 1), (3, 0)) == (4, 3, 1, 0)
    @test SumOfSquares.Certificate._merge_sorted([4, 1], [3, 2]) == [4, 3, 2, 1]
    @test SumOfSquares.Certificate._merge_sorted((4, 1), (3, 2)) == (4, 3, 2, 1)
end

@testset "with_variables" begin
    @polyvar x y z
    p = x + z
    v = SumOfSquares.Certificate.with_variables(p, y)
    @test v.inner === p
    @test MP.variables(v.inner) == [x, z]
    @test MP.variables(v) == [x, y, z]

    @ncpolyvar a b
    q = a * b + b * a
    v = SumOfSquares.Certificate.with_variables(q, FullSpace())
    @test v.inner === q
    @test MP.variables(v.inner) == [a, b, a]
    @test MP.variables(v) == [a, b, a]
end

@testset "Monomial selection for certificate" begin
    @polyvar x y z
    @ncpolyvar a b
    @testset "Multipartite error not commutative" for parts in [([a],), ([a], [b])]
        err = ArgumentError("Multipartite Newton polytope not supported with noncommutative variables.")
        @test_throws err SumOfSquares.Certificate.monomials_half_newton_polytope([a*b, b^2], parts)
    end
    @testset "Multipartite error not disjoint: $parts" for parts in [([x], [x]),
                                                                     ([x], [x, y]),
                                                                     ([x], [y], [x, y])]
        err = ArgumentError("Parts are not disjoint in multipartite Newton polytope estimation: $parts.")
        @test_throws err SumOfSquares.Certificate.monomials_half_newton_polytope([x*y, y^2], parts)
    end
    @testset "Unipartite" begin
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x*y, y^2], tuple()) == [y]
        @test isempty(SumOfSquares.Certificate.monomials_half_newton_polytope([x, y], tuple()))
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], tuple()) == [x, y]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], ([x, y],)) == [x, y]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], ([y, x],)) == [x, y]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, x^3*y^2, x^4*y^4], tuple()) == [x^2*y^2, x]
    end
    @testset "Non-commutative" begin
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([a^4, a^3*b, a*b*a^2, a*b*a*b], tuple()) == [a^2, a*b]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([a^2, a^10*b^20*a^11, a^11*b^20*a^10, a^10*b^20*a^20*b^20*a^10], tuple()) == [a^10*b^20*a^10, a]
    end
    @testset "Multipartite" begin
        # In the part [y, z], the degree is between 0 and 2
        X = [x^4, x^2*y^2, x^2*z^2, x^2*y*z, y*z]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, tuple(), apply_post_filter = false) == [x^2, x*y, x*z, y*z, x, y, z]
        function full_test(X, Y, part1, part2)
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part1,), apply_post_filter = false) == Y
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2,), apply_post_filter = false) == Y
            a = SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2,), apply_post_filter = false)
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part1, part2), apply_post_filter = false) == Y
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2, part1), apply_post_filter = false) == Y
        end
        full_test(X, monovec([x^2, x*y, x*z,      x, y, z]), [x], [y, z])
        full_test(X, monovec([x^2, x*y, x*z, y*z, x,    z]), [y], [x, z])
        full_test(X, monovec([x^2, x*y, x*z, y*z, x, y   ]), [z], [x, y])
        # FIXME: With recursive merging, it should give [x^2, x*y, x*z, x]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^4, x^2*y^2, x^2*z^2, x^2*y*z, y*z], ([x], [y], [z]), apply_post_filter = false) == [x^2, x*y, x*z, y*z, x, y, z]
    end
end

@testset "Random SOS should be SOS" begin
    @polyvar x y
    x = [1, x, y, x^2, y^2, x*y]
    @test_throws ArgumentError randsos(x, monotype=:Unknown)
    for i in 1:10
        for monotype in [:Classic, :Gram]
            p = randsos(x, monotype=monotype)
            @test p isa GramMatrix
            @test isposdef(Matrix(p.Q))
        end
    end
end

function _certificate_api(certificate::Certificate.AbstractCertificate)
    @test Certificate.cone(certificate) isa SumOfSquares.SOSLikeCone
    @test SumOfSquares.matrix_cone_type(typeof(certificate)) <: MOI.AbstractVectorSet
end
function _basis_check_each(basis::MB.AbstractPolynomialBasis, basis_type)
    @test basis isa basis_type
    if basis isa MB.AbstractMonomialBasis
        # This fails if `basis` is `Vector{Monomial{true}}` instead of `MonomialVector{true}`
        # for DynamicPolynomials. This is important as
        # `polynomial(::AbstractMatrix, ::MonomialVector, ::Type)` is implemented but
        # `polynomial(Q::AbstractMatrix, X::AbstractVector, ::Type)` falls back to
        # `dot(X, Q * X)` which does not work (`promote_operation` calls `zero(eltype(X))`
        # which gives `Polynomial{true, Int}` which then tries to multiply a
        # `ScalarAffineFunction{Float64}` with an `Int`).
        monos = basis.monomials
        @test typeof(monos) == typeof(monovec(monos))
        @test issorted(monos, rev=true)
    end
end
function _basis_check(basis, basis_type)
    @test basis isa MB.AbstractPolynomialBasis || basis isa Vector{<:MB.AbstractPolynomialBasis}
    if basis isa Vector
        # FIXME `basis_type` is `Vector{MB.MonomialBasis}` instead of `Vector{MB.MonomialBasis{...}}`
        # Once this is fixed, we should check
        # @test basis isa basis_type
        for b in basis
            _basis_check_each(b, eltype(basis_type))
        end
    else
        _basis_check_each(basis, basis_type)
    end
end

function certificate_api(certificate::Certificate.AbstractIdealCertificate)
    _certificate_api(certificate)
    @polyvar x
    poly = x + 1
    domain = @set x == 1
    @test Certificate.reduced_polynomial(certificate, poly, domain) isa MP.AbstractPolynomial
    _basis_check(Certificate.gram_basis(certificate, poly),
                 Certificate.gram_basis_type(typeof(certificate)))
    zbasis = Certificate.zero_basis(certificate)
    @test zbasis <: MB.AbstractPolynomialBasis
    @test zbasis == Certificate.zero_basis_type(typeof(certificate))
end

function certificate_api(certificate::Certificate.AbstractPreorderCertificate)
    _certificate_api(certificate)
    @polyvar x
    poly = x + 1
    domain = @set x >= 1
    processed = Certificate.preprocessed_domain(certificate, domain, poly)
    for idx in Certificate.preorder_indices(certificate, processed)
        _basis_check(Certificate.multiplier_basis(certificate, idx, processed),
                     Certificate.multiplier_basis_type(typeof(certificate)))
        @test Certificate.generator(certificate, idx, processed) isa MP.AbstractPolynomial
    end
    icert = Certificate.ideal_certificate(certificate)
    @test icert isa Certificate.AbstractIdealCertificate
    @test typeof(icert) == Certificate.ideal_certificate(typeof(certificate))
end

@testset "API" begin
    @polyvar x
    cone = SumOfSquares.SOSCone()
    BT = MB.MonomialBasis
    maxdegree = 2
    function _test(certificate::Certificate.AbstractIdealCertificate)
        certificate_api(certificate)
        preorder = Certificate.Putinar(certificate, cone, BT, maxdegree)
        certificate_api(preorder)
        sparsities = Sparsity.Pattern[Sparsity.Variable()]
        if certificate isa Certificate.MaxDegree
            push!(sparsities, Sparsity.Monomial(ChordalCompletion(), 1))
        end
        @testset "$(typeof(sparsity))" for sparsity in sparsities
            certificate_api(Certificate.Sparsity.Preorder(sparsity, preorder))
        end
    end
    basis = BT(monovec([x^2, x]))
    @testset "$(typeof(certificate))" for certificate in [
        Certificate.MaxDegree(cone, BT, maxdegree),
        Certificate.FixedBasis(cone, basis),
        Certificate.Newton(cone, BT, tuple())
    ]
        _test(certificate)
        _test(Certificate.Remainder(certificate))
        if certificate isa Certificate.MaxDegree
            _test(Certificate.Sparsity.Ideal(Sparsity.Variable(), certificate))
        end
        @testset "$(typeof(sparsity))" for sparsity in [SignSymmetry(), Sparsity.Monomial(ChordalCompletion(), 1)]
            _test(Certificate.Sparsity.Ideal(sparsity, certificate))
            _test(Certificate.Sparsity.Ideal(sparsity, Certificate.Remainder(certificate)))
        end
    end
end

include("ceg_test.jl")
include("csp_test.jl")
include("sparsity.jl")
include("symmetry.jl")

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# Autogenerated wrapper script for libsqm_jll for x86_64-linux-gnu-libgfortran3
export libsqm

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsqm")
JLLWrappers.@declare_library_product(libsqm, "libsqm.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsqm,
        "lib/libsqm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using FMM3D
using Random
using LinearAlgebra

#  This is a sample code to demonstrate how to use
#  the fmm libraries
#

# sample with one density, sources to sources,
# charge interactions, and potential only
#
n = 2000
nd = 1
sources = rand(3,n)
eps = 10.0^(-5)

charges = rand(n)
vals = lfmm3d(eps,sources,charges=charges,pg=1)

# sample with a vector of densities, sources to
# sources and targets, dipole interactions,
# potential and gradietns

nd = 3
nt = 1870
targ = rand(3,nt)
dipvecs = rand(nd,3,n)
vals = lfmm3d(eps,sources,dipvecs=dipvecs,
             targets=targ,nd=nd,pg=2,pgt=2)


text:
using Gadfly

n = 10
plot(x=rand(n), y=rand(n), color=sqrt(rand(n)), Scale.color_sqrt)


text:
module RefFEsTestsAll

using Test

@testset "PolynomialBases" begin include("PolynomialBasesTests.jl") end

@testset "Polytopes" begin include("PolytopesTests.jl") end

@testset "DOFBases" begin include("DOFBasesTests.jl") end

@testset "RefFEs" begin include("RefFEsTests.jl") end

@testset "PDiscRefFEs" begin include("PDiscRefFEsTests.jl") end

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

end # module


text:
struct DiffEqSystem <: AbstractSystem
    eqs::Vector{Operation}
    ivs::Vector{Variable}
    dvs::Vector{Variable}
    vs::Vector{Variable}
    ps::Vector{Variable}
    iv_name::Symbol
    dv_name::Symbol
    p_name::Symbol
end

function DiffEqSystem(eqs, ivs, dvs, vs, ps)
    iv_name = ivs[1].subtype
    dv_name = dvs[1].subtype
    p_name = isempty(ps) ? :Parameter : ps[1].subtype
    DiffEqSystem(eqs, ivs, dvs, vs, ps, iv_name, dv_name, p_name)
end

function DiffEqSystem(eqs; iv_name = :IndependentVariable,
                           dv_name = :DependentVariable,
                           v_name = :Variable,
                           p_name = :Parameter)
    targetmap =  Dict(iv_name => iv_name, dv_name => dv_name, v_name => v_name,
                       p_name => p_name)
    ivs, dvs, vs, ps = extract_elements(eqs, targetmap)
    DiffEqSystem(eqs, ivs, dvs, vs, ps, iv_name, dv_name, p_name)
end

function DiffEqSystem(eqs, ivs;
                      dv_name = :DependentVariable,
                      v_name = :Variable,
                      p_name = :Parameter)
    targetmap =  Dict(dv_name => dv_name, v_name => v_name, p_name => p_name)
    dvs, vs, ps = extract_elements(eqs, targetmap)
    DiffEqSystem(eqs, ivs, dvs, vs, ps, ivs[1].subtype, dv_name, p_name)
end

function generate_ode_function(sys::DiffEqSystem)
    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in 1:length(sys.dvs)]
    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
    sys_exprs = build_equals_expr.(sys.eqs)
    dvar_exprs = [:(du[$i] = $(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in 1:length(sys.dvs)]
    exprs = vcat(var_exprs,param_exprs,sys_exprs,dvar_exprs)
    block = expr_arr_to_block(exprs)
    :((du,u,p,t)->$(block))
end

isintermediate(eq) = eq.args[1].diff == nothing

function build_equals_expr(eq)
    @assert typeof(eq.args[1]) <: Variable
    if !(isintermediate(eq))
        # Differential statement
        :($(Symbol("$(eq.args[1].name)_$(eq.args[1].diff.x.name)")) = $(eq.args[2]))
    else
        # Intermediate calculation
        :($(Symbol("$(eq.args[1].name)")) = $(eq.args[2]))
    end
end

function calculate_jacobian(sys::DiffEqSystem,simplify=true)
    diff_idxs = map(eq->eq.args[1].diff !=nothing,sys.eqs)
    diff_exprs = sys.eqs[diff_idxs]
    rhs = [eq.args[2] for eq in diff_exprs]
    # Handle intermediate calculations by substitution
    calcs = sys.eqs[.!(diff_idxs)]
    for i in 1:length(calcs)
        find_replace!.(rhs,calcs[i].args[1],calcs[i].args[2])
    end
    sys_exprs = calculate_jacobian(rhs,sys.dvs)
    sys_exprs = Expression[expand_derivatives(expr) for expr in sys_exprs]
    if simplify
        sys_exprs = Expression[simplify_constants(expr) for expr in sys_exprs]
    end
    sys_exprs
end

function generate_ode_jacobian(sys::DiffEqSystem,simplify=true)
    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in 1:length(sys.dvs)]
    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
    diff_idxs = map(eq->eq.args[1].diff !=nothing,sys.eqs)
    diff_exprs = sys.eqs[diff_idxs]
    jac = calculate_jacobian(sys,simplify)
    jac_exprs = [:(J[$i,$j] = $(Expr(jac[i,j]))) for i in 1:size(jac,1), j in 1:size(jac,2)]
    exprs = vcat(var_exprs,param_exprs,vec(jac_exprs))
    block = expr_arr_to_block(exprs)
    :((J,u,p,t)->$(block))
end

function DiffEqBase.DiffEqFunction(sys::DiffEqSystem)
    expr = generate_ode_function(sys)
    DiffEqFunction{true}(eval(expr))
end

export DiffEqSystem, DiffEqFunction
export generate_ode_function


text:
const _SHOW_LIMIT = 20

_values(od::OrderedDict) = od.vals # This might be unstable

function format_idkey(key)
    reg = Regex("[^A-Za-z0-9]")
    key = replace(key, reg => "_")
    return key
end

function _preview(io::IO, str, hlim = displaysize(io)[2])
    hlim = max(40, floor(Int, hlim * 0.8))
    (length(str) <= hlim) ? str : string(SubString(str, 1, hlim), "...")
end

function _show_data_preview(f::Function, io::IO, col)
    vlim = max(20, displaysize(io)[1])
    for (i, dat) in enumerate(col)
        val = string(f(dat))
        str = _preview(io, val)
        print(io, "\n[", i, "] ", str)
        if i == vlim
            print(io, "\n[...]")
            break
        end
    end
end

function _push_csv!(col, str::String)
    vals = split(str, ",", keepempty = false)
    for vali in vals
        push!(col, strip(vali))
    end
    return col
end

function _push_csv!(col, str::String, strs::String...)
    _push_csv!(col, str)
    for stri in strs
        _push_csv!(col, stri)
    end
    return col
end


text:
function plot_comp_static_firm(x, y; xlabel = "", ylabel = "", Filename::String, Folder::String)
    y_LU = y[1]
    y_BU = y[2]
    y_LS = y[3]
    y_BS = y[4]

    plot(x, y_LU ./ y_LU[1], xlabel = xlabel, ylabel = ylabel, label = "Small, unsophisticated", color = :black, ls = :dash)
    plot!(x, y_BU ./ y_BU[1], color = :black, label = "Big, unsophisticated")
    plot!(x, y_LS ./ y_LS[1], color = :red, ls = :dash, label = "Small, sophisticated")
    plot!(x, y_BS ./ y_BS[1], color = :red, label = "Big, sophisticated")

    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_normalized.pdf")
    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_normalized.png")

    plot(x, y_LU, color = :black, ls = :dash, xlabel = xlabel, ylabel = ylabel, label = "Small, unsophisticated")
    plot!(x, y_BU, color = :black, label = "Big, unsophisticated")
    plot!(x, y_LS, color = :red, ls = :dash, label = "Small, sophisticated")
    plot!(x, y_BS, color = :red, label = "Big, sophisticated")

    savefig("../Graphs/v3/" * Folder * "/" * Filename * ".pdf")
    savefig("../Graphs/v3/" * Folder * "/" * Filename * ".png")
end
function plot_comp_static_firm_export(x, y; xlabel = "", ylabel = "", Filename::String, Folder::String, legend = :no, norm = :no)
    if norm == :yes
        y_LU = y[1]./y[1][1]
        y_BU = y[2]./y[2][1]
        y_LS = y[3]./y[3][1]
        y_BS = y[4]./y[4][1]
    elseif norm == :no
        y_LU = y[1]
        y_BU = y[2]
        y_LS = y[3]
        y_BS = y[4]
    else
        error("Keyword norm needs to be either :yes or :no")
    end

    if legend == :no
        p = plot(x, y_LU, color = :black, ls = :dash, xlabel = xlabel, ylabel = ylabel, label = "")
        plot!(x, y_BU, color = :black, label = "")
        plot!(x, y_LS, color = :red, ls = :dash, label = "")
        plot!(x, y_BS, color = :red, label = "")
    elseif legend == :yes
        p = plot(x, y_LU, color = :black, ls = :dash, xlabel = xlabel, ylabel = ylabel, label = "Small, unsophisticated")
        plot!(x, y_BU, color = :black, label = "Big, unsophisticated")
        plot!(x, y_LS, color = :red, ls = :dash, label = "Small, sophisticated")
        plot!(x, y_BS, color = :red, label = "Big, sophisticated")
    end
    return p
end
function plot_comp_static_firm1d(x, y; xlabel="", ylabel="", Filename::String, Folder::String)
    y_L = y[1]
    y_B = y[2]

    plot(x, y_L ./ y_L[1], xlabel=xlabel, ylabel=ylabel, label="low data collection ability", color=:black)
    plot!(x, y_B ./ y_B[1], color=:red, label="high data collection ability")

    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_normalized_1d.pdf")
    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_normalized_1d.png")

    plot(x, y_L, xlabel=xlabel, ylabel=ylabel, label="low data collection ability", color=:black)
    plot!(x, y_B, color=:red, label="high data collection ability")

    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_1d.pdf")
    savefig("../Graphs/v3/" * Folder * "/" * Filename * "_1d.png")
end
function plot_comp_static_firm_export1d(x, y; xlabel = "", ylabel = "", Filename::String, Folder::String, legend = :no, norm = :no)
    if norm == :yes
        y_L = y[1]./y[1][1]
        y_B = y[2]./y[2][1]
    elseif norm == :no
        y_L = y[1]
        y_B = y[2]
    else
        error("Keyword norm needs to be either :yes or :no")
    end

    if legend == :no
        p = plot(x, y_L, color=:black, xlabel=xlabel, ylabel=ylabel, label="")
        plot!(x, y_B, color=:red, label="")
    elseif legend == :yes
        p = plot(x, y_L, color=:black, xlabel=xlabel, ylabel=ylabel, label="low data collection ability")
        plot!(x, y_B, color=:red, label="high data collection ability")
    end
    return p
end


text:
"Encapsulates the state of a given player"
type PlayerState
    loc::Loc
    walls::Int
end

PlayerState(loc::Loc) = PlayerState(loc, 10)

loc(ps::PlayerState) = ps.loc
walls(ps::PlayerState) = ps.walls

function potential_moves(ps::PlayerState)
    # TODO rewrite in terms of adjacent(loc::Loc). might be slower?
    our_loc = loc(ps)
    Vector4{PlayerMove}(PlayerMove(our_loc, above(our_loc)),
                        PlayerMove(our_loc, below(our_loc)),
                        PlayerMove(our_loc, right(our_loc)),
                        PlayerMove(our_loc, left(our_loc)))
end

function moveto!(p::PlayerState, loc::Loc)
    p.loc = loc
end


text:
using CLIMAParameters.Planet: R_d, cv_d, T_0, e_int_v0, e_int_i0

"""
    linearized_air_pressure(ρ, ρe_tot, ρe_pot, ρq_tot=0, ρq_liq=0, ρq_ice=0)

The air pressure, linearized around a dry rest state, from the equation of state
(ideal gas law) where:

 - `ρ` (moist-)air density
 - `ρe_tot` total energy density
 - `ρe_pot` potential energy density
and, optionally,
 - `ρq_tot` total water density
 - `ρq_liq` liquid water density
 - `ρq_ice` ice density
"""
function linearized_air_pressure(
    param_set::AbstractParameterSet,
    ρ::FT,
    ρe_tot::FT,
    ρe_pot::FT,
    ρq_tot::FT = FT(0),
    ρq_liq::FT = FT(0),
    ρq_ice::FT = FT(0),
) where {FT <: Real, PS}
    _R_d::FT = R_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _T_0::FT = T_0(param_set)
    _e_int_v0::FT = e_int_v0(param_set)
    _e_int_i0::FT = e_int_i0(param_set)
    return ρ * _R_d * _T_0 +
           _R_d / _cv_d * (
        ρe_tot - ρe_pot - (ρq_tot - ρq_liq) * _e_int_v0 +
        ρq_ice * (_e_int_i0 + _e_int_v0)
    )
end

@inline function linearized_pressure(
    ::DryModel,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    return linearized_air_pressure(param_set, state.ρ, state.ρe, ρe_pot)
end
@inline function linearized_pressure(
    ::EquilMoist,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(
        param_set,
        state.ρ,
        state.ρe,
        ρe_pot,
        state.moisture.ρq_tot,
    )
end
@inline function linearized_pressure(
    ::NonEquilMoist,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(
        param_set,
        state.ρ,
        state.ρe,
        ρe_pot,
        state.moisture.ρq_tot,
        state.moisture.ρq_liq,
        state.moisture.ρq_ice,
    )
end

abstract type AtmosLinearModel <: BalanceLaw end

function vars_state(lm::AtmosLinearModel, st::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        turbulence::vars_state(lm.atmos.turbulence, st, FT)
        hyperdiffusion::vars_state(lm.atmos.hyperdiffusion, st, FT)
        moisture::vars_state(lm.atmos.moisture, st, FT)
    end
end
vars_state(lm::AtmosLinearModel, ::AbstractStateType, FT) = @vars()
vars_state(lm::AtmosLinearModel, st::Auxiliary, FT) =
    vars_state(lm.atmos, st, FT)

function update_auxiliary_state!(
    dg::DGModel,
    lm::AtmosLinearModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end
function flux_second_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
integral_set_auxiliary_state!(lm::AtmosLinearModel, aux::Vars, integ::Vars) =
    nothing
reverse_integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
reverse_integral_set_auxiliary_state!(
    lm::AtmosLinearModel,
    aux::Vars,
    integ::Vars,
) = nothing
flux_second_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
) = nothing
function wavespeed(
    lm::AtmosLinearModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ref = aux.ref_state
    return soundspeed_air(lm.atmos.param_set, ref.T)
end

function boundary_state!(
    nf::NumericalFluxFirstOrder,
    atmoslm::AtmosLinearModel,
    args...,
)
    atmos_boundary_state!(nf, AtmosBC(), atmoslm, args...)
end
function boundary_state!(
    nf::NumericalFluxSecondOrder,
    atmoslm::AtmosLinearModel,
    args...,
)
    nothing
end
init_state_auxiliary!(
    lm::AtmosLinearModel,
    aux::MPIStateArray,
    grid,
    direction,
) = nothing
init_state_prognostic!(
    lm::AtmosLinearModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end

function flux_first_order!(
    lm::AtmosAcousticLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL = linearized_pressure(
        lm.atmos.moisture,
        lm.atmos.param_set,
        lm.atmos.orientation,
        state,
        aux,
    )
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ - e_pot) * state.ρu
    nothing
end
source!(::AtmosAcousticLinearModel, _...) = nothing

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticGravityLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("AtmosAcousticGravityLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end
function flux_first_order!(
    lm::AtmosAcousticGravityLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)

    flux.ρ = state.ρu
    pL = linearized_pressure(
        lm.atmos.moisture,
        lm.atmos.param_set,
        lm.atmos.orientation,
        state,
        aux,
    )
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ) * state.ρu
    nothing
end
function source!(
    lm::AtmosAcousticGravityLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    ::NTuple{1, Dir},
) where {Dir <: Direction}
    if Dir === VerticalDirection || Dir === EveryDirection
        ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
        source.ρu -= state.ρ * ∇Φ
    end
    nothing
end

function numerical_flux_first_order!(
    numerical_flux::RoeNumericalFlux,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}
    @assert balance_law.atmos.moisture isa DryModel

    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    atmos = balance_law.atmos
    param_set = atmos.param_set

    ρu⁻ = state_prognostic⁻.ρu

    ref_ρ⁻ = state_auxiliary⁻.ref_state.ρ
    ref_ρe⁻ = state_auxiliary⁻.ref_state.ρe
    ref_T⁻ = state_auxiliary⁻.ref_state.T
    ref_p⁻ = state_auxiliary⁻.ref_state.p
    ref_h⁻ = (ref_ρe⁻ + ref_p⁻) / ref_ρ⁻
    ref_c⁻ = soundspeed_air(param_set, ref_T⁻)

    pL⁻ = linearized_pressure(
        atmos.moisture,
        param_set,
        atmos.orientation,
        state_prognostic⁻,
        state_auxiliary⁻,
    )

    ρu⁺ = state_prognostic⁺.ρu

    ref_ρ⁺ = state_auxiliary⁺.ref_state.ρ
    ref_ρe⁺ = state_auxiliary⁺.ref_state.ρe
    ref_T⁺ = state_auxiliary⁺.ref_state.T
    ref_p⁺ = state_auxiliary⁺.ref_state.p
    ref_h⁺ = (ref_ρe⁺ + ref_p⁺) / ref_ρ⁺
    ref_c⁺ = soundspeed_air(param_set, ref_T⁺)

    pL⁺ = linearized_pressure(
        atmos.moisture,
        param_set,
        atmos.orientation,
        state_prognostic⁺,
        state_auxiliary⁺,
    )

    # not sure if arithmetic averages are a good idea here
    h̃ = (ref_h⁻ + ref_h⁺) / 2
    c̃ = (ref_c⁻ + ref_c⁺) / 2

    ΔpL = pL⁺ - pL⁻
    Δρuᵀn = (ρu⁺ - ρu⁻)' * normal_vector

    fluxᵀn.ρ -= ΔpL / 2c̃
    fluxᵀn.ρu -= c̃ * Δρuᵀn * normal_vector / 2
    fluxᵀn.ρe -= h̃ * ΔpL / 2c̃
end

function numerical_flux_first_order!(
    ::HLLCNumericalFlux,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    # There is no intermediate speed for the AtmosLinearModel.
    # As a result, HLLC simplifies to Rusanov.
    numerical_flux_first_order!(
        RusanovNumericalFlux(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
##
## --------------------------------------------------
##
# Macros to register POLY{T},  POLY{αs..., T}
#
# We use `Vector{T}` with N =  d+1 to store  a degree d  polynomial - no trailing zeros permitted.
##
macro register0(name, parent)
    poly = esc(name)
    parent_type = esc(parent)
    quote
        struct $poly{T,N} <: $parent_type{T,N}
            coeffs::Vector{T}
            var::Symbol
            function $poly{T,N}(coeffs::Vector{T}, var::Polynomials.SymbolLike=:x) where {T, N}
                M = length(coeffs)
                (M  != N || (N > 0 &&iszero(coeffs[end]))) &&  throw(ArgumentError("wrong  size"))
                new{T,N}(coeffs, Symbol(var))

            end
            function $poly{T,N}(coeffs::NTuple{N,T}, var::Polynomials.SymbolLike=:x) where {T, N}
                M = length(coeffs)
                (M  != N || (N > 0 &&iszero(coeffs[end]))) &&  throw(ArgumentError("wrong  size"))
                new{T,N}(collect(coeffs), Symbol(var))
            end

            function $poly{T}(coeffs::Vector{S},  var::Polynomials.SymbolLike=:x) where {T,S}
                N = findlast(!iszero, coeffs)
                if N ==  nothing
                    new{T,0}(T[], Symbol(var))
                else
                    cs = T.(coeffs[1:N])
                    new{T,N}(cs,  Symbol(var))
                end
            end

            function $poly(coeffs::Vector{S},  var::Polynomials.SymbolLike=:x) where {S}
                $poly{S}(coeffs, var)
            end

        end

        Base.length(p::$poly{T,N}) where {T,N} = N
        (p::$poly)(x::S) where  {S} = eval_cop(typeof(p), p.coeffs, x)

        Polynomials.@register $poly

        # work around N parameter in promote(p,q) usage in defaults
        Base.:*(p::$poly{T}, q::$poly{T}) where {T}  = ⊗(p,q)
        Base.:+(p::$poly{T}, q::$poly{T}) where {T}  = ⊕(p,q)
        Base.divrem(p::$poly{T}, q::$poly{T}) where {T}  = _divrem(p,q)
    end
end

# register macro for polynomial families with parameters
# would be nice to consolidate, but ...
macro registerN(name,  parent, params...)
    poly = esc(name)
    parent_type = esc(parent)
    αs = tuple(esc.(params)...)
    quote
        struct $poly{$(αs...), T,N} <: $parent_type{$(αs...), T,N}
            coeffs::Vector{T}
            var::Symbol
            function $poly{$(αs...), T,N}(coeffs::Vector{T}, var::Polynomials.SymbolLike=:x) where {$(αs...), T, N}
                M = length(coeffs)
                (M  != N || (N > 0 &&iszero(coeffs[end]))) &&  throw(ArgumentError("wrong  size"))
                new{$(αs...), T,N}(coeffs, Symbol(var))
            end

            function $poly{$(αs...),T,N}(coeffs::NTuple{N,T}, var::Polynomials.SymbolLike=:x) where {$(αs...), T, N}
                M = length(coeffs)
                (M  != N || (N > 0 &&iszero(coeffs[end]))) &&  throw(ArgumentError("wrong  size"))
                new{$(αs...),T,N}(collect(coeffs), Symbol(var))
            end

            function $poly{$(αs...),T}(coeffs::Vector{S},  var::Polynomials.SymbolLike=:x) where {$(αs...),T,S}
                N = findlast(!iszero, coeffs)
                if N ==  nothing
                    new{$(αs...),T,0}(T[], Symbol(var))
                else
                    new{$(αs...),T,N}(T.(coeffs[1:N]), Symbol(var))
                end
            end

            function $poly{$(αs...)}(coeffs::Vector{T},  var::Polynomials.SymbolLike=:x) where {$(αs...), T}
                $poly{$(αs...),T}(coeffs, var)
            end

        end

        Base.length(ch::$poly{$(αs...),T,N}) where {$(αs...),T,N} = N
        (ch::$poly)(x::S) where  {S} = eval_cop(typeof(ch),  ch.coeffs, x)

        Base.convert(::Type{P}, q::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = q
        Base.convert(::Type{$poly{$(αs...)}}, q::Q) where {$(αs...),T, Q <: $poly{$(αs...),T}} = q
        Base.promote(p1::P, p2::Q) where {$(αs...),T, P<:$poly{$(αs...),T}, Q <: $poly{$(αs...),T}} = p1,p2
        Base.promote_rule(::Type{<:$poly{$(αs...),T}}, ::Type{<:$poly{$(αs...),S}}) where {$(αs...),T,S} =
            $poly{$(αs...),promote_type(T, S)}
        Base.promote_rule(::Type{<:$poly{$(αs...),T}}, ::Type{S}) where {$(αs...),T,S<:Number} =
            $poly{$(αs...),promote_type(T,S)}

        $poly{$(αs...),T}(n::Number, var::Polynomials.SymbolLike = :x) where {$(αs...),T} =
            n * one($poly{$(αs...),T}, var)
        $poly{$(αs...)}(n::S, var::Polynomials.SymbolLike = :x)        where {$(αs...), S<:Number} =
            $poly{$(αs...),S}(n,var)
        $poly{$(αs...),T}(var::Polynomials.SymbolLike=:x)              where {$(αs...), T} =
            variable($poly{$(αs...),T}, var)
        $poly{$(αs...)}(var::Polynomials.SymbolLike=:x)                where {$(αs...)} =
            variable($poly{$(αs...)}, var)

        # work around N parameter in promote
        Base.:*(p1::$poly{$(αs...),T}, p2::$poly{$(αs...),T}) where {$(αs...),T}  = ⊗(p1,p2)
        Base.:+(p1::$poly{$(αs...),T}, p2::$poly{$(αs...),T}) where {$(αs...),T}  = ⊕(p1,p2)
        Base.divrem(p1::$poly{$(αs...),T}, p2::$poly{$(αs...),T}) where {$(αs...),T}  = _divrem(p1,p2)

    end
end

## Delegation macros
"""
    ϟ(P) = Q; ϟ(P{T}) where {T} =   Q{T}

Used to delegate methods of P to Q.  ([slash]upkoppa[tab])
"""
ϟ(P::Type{<:Polynomials.AbstractPolynomial}) = throw(ArgumentError("No default for  delegation"))

## Make monic version
macro register_monic(name)

    monic=esc(name)

    quote
        SpecialPolynomials.ismonic(::Type{<:$monic}) = true
        SpecialPolynomials.abcde(P::Type{<:$monic})  = abcde(ϟ(P))
        SpecialPolynomials.basis_symbol(P::Type{<:$monic})  = basis_symbol(ϟ(P)) * "̃"
        SpecialPolynomials.weight_function(P::Type{<:$monic}) = weight_function(ϟ(P))
        Polynomials.domain(P::Type{<:$monic}) = domain(ϟ(P))
        SpecialPolynomials.B̃n(P::Type{<:$monic}, n::Int) =  B̃n(ϟ(P),n)
        SpecialPolynomials.B̃n(P::Type{<:$monic}, v::Val{N}) where {N} =  B̃n(ϟ(P),v)
        SpecialPolynomials.C̃n(P::Type{<:$monic}, n::Int) =  C̃n(ϟ(P),n)
        SpecialPolynomials.C̃n(P::Type{<:$monic}, v::Val{N}) where {N} =  C̃n(ϟ(P),v)
        SpecialPolynomials.ẫn(P::Type{<:$monic}, v::Val{N}) where {N} =  ẫn(ϟ(P),v)
        SpecialPolynomials.b̂̃n(P::Type{<:$monic}, v::Val{N}) where {N} =  b̂̃n(ϟ(P),v)
        SpecialPolynomials.ĉ̃n(P::Type{<:$monic}, v::Val{N}) where {N} =  ĉ̃n(ϟ(P),v)
    end
end

## Make ortho
macro register_orthonormal(name)

    orthonormal=esc(name)

    quote
        SpecialPolynomials.isorthonormal(::Type{<:$orthonormal}) = true
        SpecialPolynomials.abcde(P::Type{<:$orthonormal})  = abcde(ϟ(P))
        SpecialPolynomials.basis_symbol(P::Type{<:$orthonormal})  = basis_symbol(ϟ(P)) * "̃"
        SpecialPolynomials.weight_function(P::Type{<:$orthonormal}) = weight_function(ϟ(P))
        Polynomials.domain(P::Type{<:$orthonormal}) = domain(ϟ(P))
        SpecialPolynomials.k0(::Type{P}) where{P <: $orthonormal} = k0(ϟ(P)) / sqrt(norm2(ϟ(P),0))
        SpecialPolynomials.k1k0(::Type{P},  n::Int) where{P <: $orthonormal} = k1k0(ϟ(P),n) /  ω₁₀(ϟ(P), n)
        SpecialPolynomials.ω₁₀(P::Type{<:$orthonormal}, n::Int) = one(eltype(P))

        SpecialPolynomials.B̃n(P::Type{<:$orthonormal}, n::Int) =  B̃n(ϟ(P),n)
        SpecialPolynomials.B̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  B̃n(ϟ(P),v)
        SpecialPolynomials.C̃n(P::Type{<:$orthonormal}, n::Int) =  C̃n(ϟ(P),n)
        SpecialPolynomials.C̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  C̃n(ϟ(P),v)
        SpecialPolynomials.ẫn(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  ẫn(ϟ(P),v)
        SpecialPolynomials.b̂̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  b̂̃n(ϟ(P),v)
        SpecialPolynomials.ĉ̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  ĉ̃n(ϟ(P),v)
    end
end

# P̃ = P(αx+β), so P̃'' = α^2P''(α x+ β); P̃' = αP'(αx+β)
macro register_shifted(name, alpha, beta)
    shifted = esc(name)
    quote
        function SpecialPolynomials.abcde(P::Type{<:$shifted})
            α,β = $alpha, $beta
            oP = one(eltype(P))
            a,b,c,d,e = abcde(ϟ(P))
            ã = oP*a
            b̃ = (oP * 2a*β)/α + (oP*b)/α
            c̃ = (oP * (a*β^2 +b*β + c))/α^2
            d̃ = oP *  d
            ẽ = (oP * (d*β + e))/α
            NamedTuple{(:a,:b,:c,:d,:e)}((ã, b̃, c̃, d̃, ẽ))
        end

        SpecialPolynomials.basis_symbol(P::Type{<:$shifted})  = basis_symbol(ϟ(P)) * "ˢ"
        function  Polynomials.domain(P::Type{<:$shifted})
            α,β = $alpha, $beta
            a, b = extrema(ϟ(P))
            Polynomials.Interval((a- β)/α, (b - β)/α)
        end
        function SpecialPolynomials.k1k0(P::Type{<:$shifted}, n::Int)
            α,β = $alpha, $beta
            α * k1k0(ϟ(P), n)
        end

    end
end

# Weight Functions
macro register_weight_function(W, P, ω)
    W′ = esc(W)
    P′ = esc(P)
    ω′ = esc(ω)

    quote
        SpecialPolynomials.:ϟ(::Type{<:$W′}) =  $P′
        SpecialPolynomials.weight_function(::Type{<:$W′}) =  $ω′
        (ch::$W′)(x)  = clenshaw_eval(ch, x)
    end
end

macro register_discrete_weight_function(W, xs, ws)
    W′ = esc(W)
    xs′ = esc(xs)
    ws′ = esc(ws)

    quote
        SpecialPolynomials.xs_ws(::Type{<:$W′}) = ($xs′, $ws′)
        (p::$W′)(x)  = clenshaw_eval(p, x)
    end
end


text:
using Test, TrialDocs

@test hello("Julia") == "Hello, Julia"
@test domath(2.0) ≈ 7.0

max_stars_repo_name: sandreza/TrialDocs

max_issues_repo_name: sandreza/TrialDocs

max_forks_repo_name: sandreza/TrialDocs


text:
# chebyshev points and chebyshev differentiation matrix
# Ref: Spectral methods in MATLAB by Trefethen, chapters 1, 5 and 6.
function cheb_pts(n)
	return cos.((0:n)*pi/n)
end
function cheb_diff_matrix(n)
	x = cheb_pts(n)
	np = n+1
	D = zeros(np,np)
	E = view(D, reverse(1:np*np))
	D[1] = (2*n*n + 1)/6.
	E[1] = -D[1]

	nhalf = fld(np, 2)
	xint = view(x,2:n)
	D[2:n] .= (-1).^(1:n-1)./(xint .- 1)./2.0
	E[2:n] .= -D[2:n]
	D[np] = (-1)^np/2.0
	E[np] = -D[np]

	if mod(np,2) != 0
		j = nhalf + 1
		A = view(D, (j-1)*np+1:j*np)
		A .=  1.0./(x .- x[j]).*
		((-1).^((1 + j):(np + j)))
		A[j] = -x[j]/2/(1. - x[j]*x[j])
		A[1] *= 2.0
		A[np] *= 2.0
	end

	for j=2:nhalf
		A = view(D, (j-1)*np+1:j*np)
		B = view(E, (j-1)*np+1:j*np)
		A .= 1.0./(x .- x[j]).*
		((-1).^((1 + j):(np + j)))
		B .= -A
		A[1] *= 2.0
		B[1] *= 2.0
		A[j] = -x[j]/2/(1. - x[j]*x[j])
		B[j] = -A[j]
		A[np] *= 2.0
		B[np] *= 2.0
	end

	return D
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
# Admission prob is a function of a single indicator

"""
	$(SIGNATURES)

Switches governing admissions that are function of a single indicator.
"""
mutable struct AdmissionsOneVarSwitches{I1, F1 <: Real} <: AbstractAdmissionsSwitches{I1, F1}
    nColleges :: I1
    # Name of the variable that determines admissions probability
    # pctVar :: Symbol
    # Minimum probability of all college sets
    minCollSetProb :: F1
end

"""
	$(SIGNATURES)

Admissions are governed by a single indicator, such as a test score.
For each college, admissions probability is an increasing function of that indicator.
The function is stored elsewhere.
Students may be allowed to attend other colleges with a fixed probability.

College sets are of the form `1 : n`.
Probability of being admitted to all colleges = admissions probability at top college.
Probability of being admitted to `1 : n` = admissions probablity of college `n` times 1 minus admission probability at `n+1` or higher.
"""
mutable struct AdmissionsOneVar{I1, F1 <: Real} <: AbstractAdmissionsRule{I1, F1}
    switches :: AdmissionsOneVarSwitches{I1, F1}
    # admissionProbFctV :: Vector
end

function Base.show(io :: IO, a :: AdmissionsOneVar)
    nc = n_colleges(a);
    # rankVar = percentile_var(a);
    print(io, typeof(a));  # , " with ranking variable $rankVar");
end

function StructLH.describe(a :: AdmissionsOneVarSwitches)
    # rankVar = percentile_var(a);
    return [
        "Admission rule"  "based on one indicator";
        # "Admission probability determined by"  "$rankVar";
    ];
end

function validate_admissions(a :: AdmissionsOneVar{I1, F1}) where {I1, F1}
    isValid = true;
    isValid  ||  println(a)
    return isValid
end

make_admissions(switches :: AdmissionsOneVarSwitches{I1, F1}) where {I1, F1} =
    AdmissionsOneVar(switches);

make_test_adm_onevar_switches(nc) =
    AdmissionsOneVarSwitches(nc, 0.05);

# By default: With the admission prob functions inside
function make_test_admissions_onevar(nc)
    a = make_admissions(make_test_adm_onevar_switches(nc));
    return a
end

function test_admission_prob_fct(x :: Real, ic :: Integer)
    if ic == 1
        prob = 1.0;
    else
        prob = 0.7 / (1.0 + 0.2 * exp(-x));
    end
    return prob
end

## ----------  Probability of being admitted to each college set

function prob_coll_set(a :: AdmissionsOneVar{I1, F1},
    admProbFct :: AF1,
    iSet :: Integer, hsGpaPct :: F2) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    return prob_coll_sets(a, admProbFct, hsGpaPct)[iSet];
end

function prob_coll_sets(a :: AdmissionsOneVar{I1, F1},
    admProbFct :: AF1,
    hsGpaPct :: F2) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1, F2 <: Real}

    probAdmitV = [prob_admit(admProbFct, iCollege, hsGpaPct)
        for iCollege = 1 : n_colleges(a)];
    @assert all(0.0 .<= probAdmitV .<= 1.0)  "Out of bounds: $probAdmitV"
    return prob_coll_sets_from_probs(a, probAdmitV)
end

# Compute prob of being admitted to each college set from probs of being admitted to each college.
function prob_coll_sets_from_probs(a :: AdmissionsOneVar{I1, F1},
    probAdmitV :: AbstractVector{F2}) where {I1, F1, F2 <: Real}

    nSets = n_college_sets(a);
    @assert nSets == length(probAdmitV)  "Assuming college sets are 1:n";
    @assert all(0.0 .<= probAdmitV .<= 1.0)  "Out of bounds";
    @assert probAdmitV[1] == 1.0  "Expecting everyone gets into lowest college";

    probV = zeros(nSets);
    # For top college: simply its admissions probability.
    probV[nSets] = probAdmitV[nSets];
    for iSet = (nSets-1) : -1 : 1
        # For lower colleges: (prob being admitted) * (prob not admitted higher).
        probAdmitHigher = sum(probV[(iSet+1) : nSets]);
        probV[iSet] = probAdmitV[iSet] * (1.0 - probAdmitHigher);
    end
    @check sum(probV) ≈ 1.0
    return probV
end

# ----------------


text:
module CODATA2018

using PhysicalConstant, Unitful, Measurements

import PhysicalConstant: @constant, @derived_constant

@constant(α, "Fine-structure constant", 7.297_352_5664e-3,
          BigFloat(72_973_525_664)/BigFloat(10_000_000_000_000), Unitful.NoUnits,
          1.7e-12, BigFloat(17)/BigFloat(10_000_000_000_000), "CODATA 2018")
@constant(a_0, "Bohr radius", 0.529_177_210_67e-10,
          BigFloat(52_917_721_067)/BigFloat(1_000_000_000_000_000_000_000), u"m",
          1.2e-20, BigFloat(12)/BigFloat(1_000_000_000_000_000_000_000), "CODATA 2018")
@constant(atm, "Standard atmosphere", 101_325.0, BigFloat(101_325), u"Pa",
          0.0, BigFloat(0), "CODATA 2018")
@constant(b, "Wien wavelength displacement law constant", 2.897_7729e-3,
          BigFloat(28_977_729)/BigFloat(10_000_000_000), u"m * K",
          1.7e-9, BigFloat(17)/BigFloat(10_000_000_000), "CODATA 2018")
@constant(C, "Gravitational velocity in vacuum", 3.699_3e44, BigFloat(369_930_000_000_000_000_000_000_000_000_000_000_000_000_000), u"m/s",
          0.0, BigFloat(0), "CODATA 2018")
@constant(c, "Speed of light in vacuum", 299_792_458.0, BigFloat(299_792_458.0), u"m/s",
          0.0, BigFloat(0), "CODATA 2018")
#@constant(t, "Cosmic Temperature", 2.725820831,
#          BigFloat(272_582_083.1), u"K",
#          0.0, BigFloat(0), "CODATA 2018")
#@constant(t_cc, "Kotov Cosmic Periodicity", 9600.061,
#          BigFloat(960_006_1), u"s",
#          0.0, BigFloat(0), "CODATA 2018")
@constant(µ_0, "Magnetic constant", 1.2566370614359173e-6, 4*big(pi)/BigFloat(10_000_000),
          u"N * A^-2", 0.0, BigFloat(0.0), "CODATA 2018")
@constant(ε_0, "Electric constant", 8.854187817620389e-12,
          inv(ustrip(big(µ_0)) * ustrip(big(c))^2), u"F * m^-1",
          0.0, BigFloat(0.0), "CODATA 2018")
@constant(e, "Elementary charge", 1.602_176_6208e-19,
          big(16_021_766_208)/100_000_000_000_000_000_000_000_000_000,
          u"C", 9.8e-28, big(98)/100_000_000_000_000_000_000_000_000_000, "CODATA 2018")
@constant(Gg, "Newtonian constant of gravitation", 6.674_08e-11,
          big(667_408)/big(10_000_000_000_000_000), u"m^3 * kg^-1 * s^-2",
          3.1e-15, big(31)/big(10_000_000_000_000_000), "CODATA 2018")
@constant(G, "Sanchez constant of gravitation", 6.675_45e-11,
          big(667_545)/big(10_000_000_000_000_000), u"m^3 * kg^-1 * s^-2",
          3.1e-15, big(31)/big(10_000_000_000_000_000), "CODATA 2018")
@constant(g_n, "Standard acceleration of gravitation", 9.806_65, big(980_665)/big(100_000),
          u"m * s^-2", 0, 0, "CODATA 2018")
@constant(h, "Planck constant", 6.626_070_040e-34,
          6_626_070_040/10_000_000_000_000_000_000_000_000_000_000_000_000_000_000,
          u"J*s", 8.1e-42, 81/10_000_000_000_000_000_000_000_000_000_000_000_000_000_000,
          "CODATA 2018")
@derived_constant(ħ, "Planck constant over 2pi", 1.0545718001391127e-34,
                  ustrip(big(h))/(2 * big(pi)), u"J*s", measurement(h)/2pi,
                  measurement(BigFloat, h)/(2 * big(pi)), "CODATA 2018")
@constant(k_B, "Boltzmann constant", 1.380_648_52e-23,
          BigFloat(138_064_852)/BigFloat(10_000_000_000_000_000_000_000_000_000_000), u"J * K^-1",
          7.9e-30, BigFloat(79)/BigFloat(10_000_000_000_000_000_000_000_000_000_000), "CODATA 2018")
@constant(µ_B, "Bohr magneton", 927.400_9994e-26,
          BigFloat(9274_009_994)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000),
          u"J * T^-1", 5.7e-32,
          BigFloat(57)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000), "CODATA 2018")
@constant(m_e, "Electron mass", 9.109_383_56e-31,
          BigFloat(910_938_356)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 1.1e-38,
          BigFloat(11)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_h, "Hydrogen mass", 1.673_723_600e-27,
          BigFloat(1673_723_600)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_n, "Neutron mass", 1.674_927_471e-27,
          BigFloat(1674_927_471)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_p, "Proton mass", 1.672_621_898e-27,
          BigFloat(1672_621_898)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_m, "Muon mass", 1.836_152_67e-28,
          BigFloat(1836_152_67)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_t, "Tau mass", 3.167_735_502e-27,
          BigFloat(3167_735_502)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_u, "Atomic mass constant", 1.660_539_040e-27,
          BigFloat(1660_539_040)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.0e-35,
          BigFloat(20)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(N_A, "Avogadro constant", 6.022_140_857e23,
          BigFloat(602_214_085_700_000_000_000_000), u"mol^-1",
          7.4e15, BigFloat(7_400_000_000_000_000), "CODATA 2018")
@constant(r_0, "Bare Hydrogen radius", 0.529_177_210_30e-10,
          BigFloat(52_917_721_03)/BigFloat(1_000_000_000_000_000_000_000), u"m",
          1.2e-20, BigFloat(12)/BigFloat(1_000_000_000_000_000_000_000), "CODATA 2018")
@constant(R, "Molar gas constant", 8.314_4598, BigFloat(8_314_4598)/BigFloat(10_000_000),
          u"J * mol^-1 * K^-1", 4.8e-6, BigFloat(48)/BigFloat(10_000_000), "CODATA 2018")
@constant(R_∞, "Rydberg constant", 10_973_731.568_508,
          BigFloat(10_973_731_568_508)/BigFloat(1_000_000), u"m^-1",
          6.5e-5, BigFloat(65)/BigFloat(1_000_000), "CODATA 2018")
@constant(σ, "Stefan-Boltzmann constant", 5.670_367e-8,
          BigFloat(5670_367)/BigFloat(100_000_000_000_000), u"m^2",
          1.3e-13, BigFloat(13)/BigFloat(100_000_000_000_000), "CODATA 2018")
@constant(σ_e, "Thomson cross section", 0.665_245_871_58e-28,
          BigFloat(66_524_587_158)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          u"m^2", 9.1e-38,
          BigFloat(91)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(Z_0, "Characteristic impedance of vacuum", 376.73031346177066,
          BigFloat(1199_169_832)/BigFloat(10_000_000) * big(pi), u"Ω",
          0, BigFloat(0.0), "CODATA 2018")

end