text:
module Terms

using Printf: @printf, @sprintf
using StaticArrays: SVector
using ..Spatials: AbstractBond, pidtype, Bonds, AbstractLattice, acrossbonds
using ..DegreesOfFreedom: Index, IDFConfig, Table, OID, Operators, oidtype, otype, Boundary, coordpresent, coordabsent
using ...Interfaces: add!
using ...Prerequisites: Float, atol, rtol, decimaltostr, rawtype
using ...Prerequisites.TypeTraits: efficientoperations, indtosub, corder
using ...Prerequisites.CompositeStructures: CompositeTuple, NamedContainer
using ...Mathematics.AlgebraOverFields: SimpleID, ID, Element, Elements

import ..DegreesOfFreedom: isHermitian
import ...Interfaces: id, rank, kind, expand, expand!, dimension, update!, reset!
import ...Mathematics.AlgebraOverFields: idtype

export Subscript, Subscripts, @subscript
export Coupling, Couplings, @couplings
export TermFunction, TermAmplitude, TermCouplings, TermModulate, Term, statistics, abbr
export Parameters, GenOperators, AbstractGenerator, Generator

const wildcard = '*'
const constant = ':'
"""
    Subscript(  ipattern::NTuple{N1, Any},
                opattern::NTuple{N2, Any},
                mapping::Union{Function, Nothing}=nothing,
                constrain::Union{Function, Nothing}=nothing,
                identifier::Union{Symbol, Char}=wildcard
                ) where {N1, N2}
    Subscript{N}() where N
    Subscript(opattern::NTuple{N, Int}) where N

The subscripts of some orbital/spin degrees of freedom.
"""
struct Subscript{N1, N2, I<:Tuple, O<:Tuple, M<:Union{Function, Nothing}, C<:Union{Function, Nothing}, D<:Union{Symbol, Char}}
    ipattern::I
    opattern::O
    mapping::M
    constrain::C
    identifier::D
    function Subscript( ipattern::NTuple{N1, Any},
                        opattern::NTuple{N2, Any},
                        mapping::Union{Function, Nothing}=nothing,
                        constrain::Union{Function, Nothing}=nothing,
                        identifier::Union{Symbol, Char}=wildcard
                        ) where {N1, N2}
        new{N1, N2, typeof(ipattern), typeof(opattern), typeof(mapping), typeof(constrain), typeof(identifier)}(ipattern, opattern, mapping, constrain, identifier)
    end
end
Subscript{N}() where {N} = Subscript((wildcard,), NTuple{N, Char}(wildcard for i = 1:N), nothing, nothing, wildcard)
Subscript(opattern::NTuple{N, Int}) where N = Subscript((), opattern, nothing, nothing, constant)

"""
    ==(sub1::Subscript, sub2::Subscript) -> Bool

Judge whether two subscripts are equivalent to each other.
"""
Base.:(==)(sub1::Subscript, sub2::Subscript) = (sub1.ipattern == sub2.ipattern) && (sub1.opattern == sub2.opattern) && (sub1.identifier == sub2.identifier)

"""
    isequal(sub1::Subscript, sub2::Subscript) -> Bool

Judge whether two subscripts are equivalent to each other.
"""
function Base.isequal(sub1::Subscript, sub2::Subscript)
    isequal(sub1.ipattern, sub2.ipattern) && isequal(sub1.opattern, sub2.opattern) && isequal(sub1.identifier, sub2.identifier)
end

"""
    show(io::IO, subscript::Subscript)

Show a subscript.
"""
function Base.show(io::IO, subscript::Subscript)
    if (subscript.identifier == constant) || (subscript.identifier == wildcard)
        @printf io "(%s)" join(subscript.opattern, ", ")
    else
        @printf io "(%s) with %s" join(subscript.opattern, ", ") subscript.identifier
    end
end

"""
    rank(subscript::Subscript) -> Int
    rank(::Type{<:Subscript{N}}) where N -> Int

Get the number of the independent variables that are used to describe the subscripts of some orbital/spin degrees of freedom.
"""
rank(subscript::Subscript) = subscript |> typeof |> rank
rank(::Type{<:Subscript{N}}) where {N} = N

"""
    dimension(subscript::Subscript) -> Int
    dimension(::Type{<:(Subscript{N1, N2} where N1)}) where N2 -> Int

Get the number of the whole variables that are used to describe the subscripts of some orbital/spin degrees of freedom.
"""
dimension(subscript::Subscript) = subscript |> typeof |> dimension
dimension(::Type{<:(Subscript{N1, N2} where N1)}) where {N2} = N2

"""
    (subscript::Subscript{N})(::Val{'M'}, values::Vararg{Int, N}) where N -> NTuple{dimension(subscript), Int}
    (subscript::Subscript{N})(::Val{'C'}, values::Vararg{Int, N}) where N -> Bool

* Construct the subscripts from a set of independent variables.
* Judge whether a set of independent variables are valid to construct the subscripts.
"""
function (subscript::Subscript)() end
(subscript::Subscript{N})(::Val{'M'}, values::Vararg{Int, N}) where {N} = subscript.mapping(values...)
(subscript::Subscript{0, N, <:Tuple, <:Tuple, Nothing, Nothing})(::Val{'M'}) where {N} = subscript.opattern
(subscript::Subscript{1, N, <:Tuple, <:Tuple, Nothing, Nothing})(::Val{'M'}, value::Int) where {N} = NTuple{N, Int}(value for i = 1:N)
(subscript::Subscript{N})(::Val{'C'}, values::Vararg{Int, N}) where N = subscript.constrain(values...)
(subscript::Subscript{N1, N2, <:Tuple, <:Tuple, <:Union{Function, Nothing}, Nothing})(::Val{'C'}, values::Vararg{Int, N1}) where {N1, N2} = true

"""
    @subscript expr::Expr with constrain::Expr -> Subscript

Construct a subscript from a map and optionally with a constrain.
"""
macro subscript(expr::Expr, with::Symbol=:with, constrain::Union{Expr, Symbol}=:nothing)
    @assert with == :with "@subscript error: pattern and constrain must be separated by :with."
    subscriptexpr(expr, constrain)
end
function subscriptexpr(pattern::Expr, ::Symbol)
    @assert pattern.head == :tuple "subscriptexpr error: wrong input pattern (not a tuple)."
    ip, op = Tuple(subscriptipattern(pattern.args)), Tuple(pattern.args)
    mapname, identifier = gensym(), QuoteNode(gensym())
    return quote
        $mapname($(ip...)) = $pattern
        Subscript($ip, $op, $mapname, nothing, $identifier)
    end
end
function subscriptexpr(pattern::Expr, constrain::Expr)
    @assert pattern.head == :tuple "subscriptexpr error: wrong input pattern (not a tuple)."
    ip, op = Tuple(subscriptipattern(pattern.args)), Tuple(pattern.args)
    mapname, constrainname, identifier = gensym(), gensym(), QuoteNode(gensym())
    return quote
        $mapname($(ip...)) = $pattern
        $constrainname($(ip...)) = $(constrain)
        Subscript($ip, $op, $mapname, $constrainname, $identifier)
    end
end
function subscriptipattern(opattern)
    result = []
    for op in opattern
        @assert (op ≠ Symbol(wildcard)) && (op ≠ Symbol(constant)) "subscriptipattern error: wrong symbols."
        isa(op, Symbol) && (op ∉ result) && push!(result, op)
    end
    return result
end

"""
    Subscripts(contents::Subscript...)

A complete set of all the independent subscripts of the orbital/spin degrees of freedom.
"""
struct Subscripts{T<:Tuple, R, D} <: CompositeTuple{T}
    contents::T
    function Subscripts(contents::T) where {T<:Tuple}
        R = sum(rank(fieldtype(T, i)) for i = 1:fieldcount(T))
        D = sum(dimension(fieldtype(T, i)) for i = 1:fieldcount(T))
        new{T, R, D}(contents)
    end
end
Subscripts(contents::Subscript...) = Subscripts(contents)

"""
    rank(subscripts::Subscripts) -> Int
    rank(::Type{<:Subscripts{<:Tuple, R}}) where R -> Int

Get the total number of the independent variables of the complete subscript set.
"""
rank(subscripts::Subscripts) = subscripts |> typeof |> rank
rank(::Type{<:Subscripts{<:Tuple, R}}) where {R} = R

"""
    rank(subscripts::Subscripts, i::Int) -> Int
    rank(::Type{<:Subscripts{T}}, i::Int) where T -> Int

Get the number of the independent variables of a component of the complete subscript set.
"""
rank(subscripts::Subscripts, i::Int) = rank(typeof(subscripts), i)
rank(::Type{<:Subscripts{T}}, i::Int) where {T} = rank(fieldtype(T, i))

"""
    dimension(subscripts::Subscripts) -> Int
    dimension(::Type{<:(Subscripts{<:Tuple, R, D} where R)}) where D -> Int

Get the total number of the whole variables of the complete subscript set.
"""
dimension(subscripts::Subscripts) = subscripts |> typeof |> dimension
dimension(::Type{<:(Subscripts{<:Tuple, R, D} where R)}) where {D} = D

"""
    dimension(subscripts::Subscripts, i::Int) -> Int
    dimension(::Type{<:Subscripts{T}}, i::Int) where T -> Int

Get the total number of the whole variables of a component of the complete subscript set.
"""
dimension(subscripts::Subscripts, i::Int) = dimension(typeof(subscripts), i)
dimension(::Type{<:Subscripts{T}}, i::Int) where {T} = dimension(fieldtype(T, i))

"""
    (subscripts::Subscripts)(::Val{'M'}, values::NTuple{N, Int}) where N -> NTuple{dimension(subscripts), Int}
    (subscripts::Subscripts)(::Val{'C'}, values::NTuple{N, Int}) where N -> Bool

* Construct the complete set of subscripts from a complete set of independent variables.
* Judge whether a complete set of independent variables are valid to construct the complete subscripts.
"""
function (subscripts::Subscripts)() end
@generated function (subscripts::Subscripts)(::Val{'M'}, values::NTuple{N, Int}) where N
    @assert rank(subscripts) == N "subscripts error: dismatched ranks of input subscripts and values."
    exprs, count = [], 1
    for i = 1:length(subscripts)
        vs = Tuple(:(values[$j]) for j = count:(count+rank(subscripts, i)-1))
        push!(exprs, :(subscripts.contents[$i](Val('M'), $(vs...))...))
        count += rank(subscripts, i)
    end
    return Expr(:tuple, exprs...)
end
@generated function (subscripts::Subscripts)(::Val{'C'}, values::NTuple{N, Int}) where N
    @assert rank(subscripts) == N "subscripts error: dismatched ranks of input subscripts and values."
    (length(subscripts) == 0) && return :(true)
    count = rank(subscripts, 1)
    vs = Tuple(:(values[$j]) for j = 1:count)
    expr = :(subscripts.contents[1](Val('C'), $(vs...)))
    for i = 2:length(subscripts)
        vs = Tuple(:(values[$j]) for j = (count+1):(count+rank(subscripts, i)))
        expr = Expr(:(&&), expr, :(subscripts.contents[$i](Val('C'), $(vs...))))
        count += rank(subscripts, i)
    end
    return expr
end

"""
    *(sub1::Subscript, sub2::Subscript) -> Subscripts
    *(subs::Subscripts, sub::Subscript) -> Subscripts
    *(sub::Subscript, subs::Subscripts) -> Subscripts
    *(subs1::Subscripts, subs2::Subscripts) -> Subscripts

Get the multiplication between subscripts or complete sets of subscripts.
"""
Base.:*(sub1::Subscript, sub2::Subscript) = Subscripts(sub1, sub2)
Base.:*(subs::Subscripts, sub::Subscript) = Subscripts(subs.contents..., sub)
Base.:*(sub::Subscript, subs::Subscripts) = Subscripts(sub, subs.contents...)
Base.:*(subs1::Subscripts, subs2::Subscripts) = Subscripts(subs1.contents..., subs2.contents...)

"""
    expand(subscripts::Subscripts, dimensions::NTuple{N, Int}) where N -> SbExpand

Expand a complete set of subscripts with a given set of variable ranges.
"""
function expand(subscripts::Subscripts, dimensions::NTuple{N, Int}) where N
    @assert dimension(subscripts) == N "expand error: dismatched input dimensions $dimensions."
    dims, dcount, rcount = Vector{Int}(undef, rank(subscripts)), 0, 0
    for i = 1:length(subscripts)
        ipattern = subscripts[i].ipattern
        opattern = subscripts[i].opattern
        for j = 1:dimension(subscripts, i)
            isa(opattern[j], Int) && @assert 0 < opattern[j] <= dimensions[dcount+j] "expand error: opattern($opattern) out of range."
        end
        for j = 1:rank(subscripts, i)
            index, flag = 1, false
            while (index = findnext(isequal(ipattern[j]), opattern, index)) ≠ nothing
                flag || (dims[rcount+j] = dimensions[dcount+index]; flag = true)
                flag && @assert dimensions[dcount+index] == dims[rcount+j] "expand error: dismatched input dimensions."
                index = index + 1
            end
        end
        dcount = dcount + dimension(subscripts, i)
        rcount = rcount + rank(subscripts, i)
    end
    return SbExpand(subscripts, NTuple{rank(subscripts), Int}(dims))
end
struct SbExpand{N, S<:Subscripts, D}
    subscripts::S
    dims::NTuple{D, Int}
    function SbExpand(subscripts::Subscripts, dims::NTuple{D, Int}) where D
        @assert D == rank(subscripts) "SbExpand error: dismatched inputs."
        new{dimension(subscripts), typeof(subscripts), D}(subscripts, dims)
    end
end
Base.eltype(::Type{<:SbExpand{N}}) where {N} = NTuple{N, Int}
Base.IteratorSize(::Type{<:SbExpand}) = Base.SizeUnknown()
function Base.iterate(sbe::SbExpand, state::Int=1)
    while state <= prod(sbe.dims)
        inds = indtosub(sbe.dims, state, corder)
        sbe.subscripts(Val('C'), inds) && return (sbe.subscripts(Val('M'), inds), state+1)
        state = state + 1
    end
    return nothing
end

"""
    Coupling{V, I<:ID{SimpleID}} <: Element{V, I}

The coupling intra/inter interanl degrees of freedom at different lattice points.
"""
abstract type Coupling{V, I<:ID{SimpleID}} <: Element{V, I} end

couplingcenter(::Type{<:Coupling}, i::Int, n::Int, ::Val{1}) = 1
couplingcenter(::Type{<:Coupling}, i::Int, n::Int, ::Val{R}) where R = error("couplingcenter error: no default center for a rank-$R bond.")
@generated function couplingcenters(::Type{C}, centers::NTuple{N, Any}, ::Val{R}) where {C<:Coupling, N, R}
    exprs = [:(isa(centers[$i], Int) ? ((0 < centers[$i] <= R) ? centers[$i] : error("couplingcenters error: center out of range.")) : couplingcenter(C, $i, N, Val(R))) for i = 1:N]
    return Expr(:tuple, exprs...)
end
function couplingcenters(str::AbstractString)
    str=strip(str,' ')
    @assert (str[1] == '(') && (str[end] == ')') "couplingcenters error: wrong input pattern."
    return Tuple(parse(Int, center) for center in split(str[2:end-1], ','))
end

"""
    Couplings(cps::Coupling...)

A pack of couplings intra/inter interanl degrees of freedom at different lattice points.

Alias for `Elements{<:ID{SimpleID}, <:Coupling}`.
"""
const Couplings{I<:ID{SimpleID}, C<:Coupling} = Elements{I, C}
Couplings(cps::Coupling...) = Elements(cps...)

"""
    @couplings cps -> Couplings

Convert an expression/literal to a set of couplings.
"""
macro couplings(cps)
    result = @eval(__module__, $cps)
    return isa(result, Coupling) ? Couplings(result) : isa(result, Couplings) ? result : error("@couplings error: inputs contain strangers that are not coupling/couplings.")
end

"""
    TermFunction <: Function

Abstract type for concrete term functions.
"""
abstract type TermFunction <: Function end

"""
    ==(tf1::TermFunction, tf2::TermFunction) -> Bool

Judge whether two concrete term functions are equivalent to each other.
"""
Base.:(==)(tf1::TermFunction, tf2::TermFunction) = ==(efficientoperations, tf1, tf2)

"""
    isequal(tf1::TermFunction, tf2::TermFunction) -> Bool

Judge whether two concrete term functions are equivalent to each other.
"""
Base.isequal(tf1::TermFunction, tf2::TermFunction) = isequal(efficientoperations, tf1, tf2)

"""
    TermAmplitude(amplitude::Union{Function, Nothing}=nothing)

The function for the amplitude of a term.
"""
struct TermAmplitude{A<:Union{Function, Nothing}} <: TermFunction
    amplitude::A
    TermAmplitude(amplitude::Union{Function, Nothing}=nothing) = new{typeof(amplitude)}(amplitude)
end
(termamplitude::TermAmplitude{Nothing})(args...; kwargs...) = 1
(termamplitude::TermAmplitude{<:Function})(args...; kwargs...) = termamplitude.amplitude(args...; kwargs...)

"""
    TermCouplings(candidate::Coupling)
    TermCouplings(candidate::Couplings)
    TermCouplings(contents::Tuple{<:Tuple{Vararg{Couplings}}, <:Function})
    TermCouplings(candidates::NTuple{N, <:Couplings}, choice::Function) where N

The function for the couplings of a term.
"""
struct TermCouplings{C<:Union{Function, Couplings}} <: TermFunction
    couplings::C
    TermCouplings(coupling::Coupling) = new{Couplings{coupling.id|>typeof, coupling|>typeof}}(Couplings(coupling))
    TermCouplings(couplings::Couplings) = new{couplings|>typeof}(couplings)
    TermCouplings(couplings::Function) = new{couplings|>typeof}(couplings)
end
(termcouplings::TermCouplings{<:Couplings})(args...; kwargs...) = termcouplings.couplings
(termcouplings::TermCouplings{<:Function})(args...; kwargs...) = termcouplings.couplings(args...; kwargs...)

"""
    TermModulate(id::Symbol, modulate::Union{Function, Nothing}=nothing)
    TermModulate(id::Symbol, modulate::Bool)

The function for the modulation of a term.
"""
struct TermModulate{M<:Union{Function, Nothing}, id} <: TermFunction
    modulate::M
    TermModulate(id::Symbol, modulate::Union{Function, Nothing}=nothing) = new{typeof(modulate), id}(modulate)
    TermModulate(id::Symbol, modulate::Bool) = (@assert modulate "TermModulate error: input `modulate` must be `true`."; new{Nothing, id}(nothing))
end
(termmodulate::TermModulate{Nothing, id})(; kwargs...) where {id} = get(kwargs, id, nothing)
(termmodulate::TermModulate{<:Function})(args...; kwargs...) = termmodulate.modulate(args...; kwargs...)

"""
    Term{ST, K, R, I}(value::Any, bondkind::Any, couplings::TermCouplings, amplitude::TermAmplitude, modulate::Union{TermModulate, Nothing}, factor::Any) where {ST, K, R, I}
    Term{ST, K, R}( id::Symbol, value::Any, bondkind::Any;
                    couplings::Union{Function, Coupling, Couplings, TermCouplings},
                    amplitude::Union{Function, Nothing}=nothing,
                    modulate::Union{Function, Bool}=false,
                    ) where {ST, K, R}

A term of a quantum lattice system.
"""
mutable struct Term{statistics, kind, R, id, V, B<:Any, C<:TermCouplings, A<:TermAmplitude, M<:Union{TermModulate, Nothing}}
    value::V
    bondkind::B
    couplings::C
    amplitude::A
    modulate::M
    factor::V
    function Term{ST, K, R, I}(value::Any, bondkind::Any, couplings::TermCouplings, amplitude::TermAmplitude, modulate::Union{TermModulate, Nothing}, factor::Any) where {ST, K, R, I}
        @assert ST in ('F', 'B') "Term error: statistics must be 'F' or 'B'."
        @assert isa(K, Symbol) "Term error: kind must be a Symbol."
        @assert isa(I, Symbol) "Term error: id must be a Symbol."
        new{ST, K, R, I, typeof(value), typeof(bondkind), typeof(couplings), typeof(amplitude), typeof(modulate)}(value, bondkind, couplings, amplitude, modulate, factor)
    end
end
function Term{ST, K, R}(id::Symbol, value::Any, bondkind::Any;
                        couplings::Union{Function, Coupling, Couplings, TermCouplings},
                        amplitude::Union{Function, Nothing}=nothing,
                        modulate::Union{Function, Bool}=false
                        ) where {ST, K, R}
    isa(couplings, TermCouplings) || (couplings = TermCouplings(couplings))
    isa(amplitude, TermAmplitude) || (amplitude = TermAmplitude(amplitude))
    isa(modulate, TermModulate) || (modulate = (modulate === false) ? nothing : TermModulate(id, modulate))
    Term{ST, K, R, id}(value, bondkind, couplings, amplitude, modulate, 1)
end

"""
    statistics(term::Term) -> Char
    statistics(::Type{<:Term{ST}}) where ST -> Char

Get the statistics of a term.
"""
statistics(term::Term) = term |> typeof |> statistics
statistics(::Type{<:Term{ST}}) where {ST} = ST

"""
    kind(term::Term) -> Symbol
    kind(::Type{<:Term{ST, K} where ST}) where K -> Symbol

Get the kind of a term.
"""
kind(term::Term) = term |> typeof |> kind
kind(::Type{<:Term{ST, K} where ST}) where {K} = K

"""
    rank(term::Term) -> Int
    rank(::Type{<:Term{ST, K, R} where {ST, K}}) where R -> Int

Get the rank of a term.
"""
rank(term::Term) = term |> typeof |> rank
rank(::Type{<:Term{ST, K, R} where {ST, K}}) where {R} = R

"""
    id(term::Term) -> Symbol
    id(::Type{<:Term{ST, K, R, I} where {ST, K, R}}) where I -> Symbol

Get the id of a term.
"""
id(term::Term) = term |> typeof |> id
id(::Type{<:Term{ST, K, R, I} where {ST, K, R}}) where {I} = I

"""
    valtype(term::Term)
    valtype(::Type{<:Term{ST, K, R, I, V} where {ST, K, R, I}}) where V

Get the value type of a term.
"""
Base.valtype(term::Term) = term |> typeof |> valtype
Base.valtype(::Type{<:Term{ST, K, R, I, V} where {ST, K, R, I}}) where {V} = V

"""
    abbr(term::Term) -> Symbol
    abbr(::Type{<:Term}) -> Symbol

Get the abbreviation of the kind of a term.
"""
abbr(term::Term) = term |> typeof |> abbr
abbr(::Type{<:Term}) = :tm

"""
    isHermitian(term::Term) -> Bool
    isHermitian(::Type{<:Term}) -> Bool
"""
isHermitian(term::Term) = term |> typeof |> isHermitian
isHermitian(::Type{<:Term}) = error("isHermitian error: not implemented.")

"""
    ==(term1::Term, term2::Term) -> Bool

Judge whether two terms are equivalent to each other.
"""
Base.:(==)(term1::Term, term2::Term) = ==(efficientoperations, term1, term2)

"""
    isequal(term1::Term, term2::Term) -> Bool

Judge whether two terms are equivalent to each other.
"""
Base.isequal(term1::Term, term2::Term) = isequal(efficientoperations, term1, term2)

"""
    show(io::IO, term::Term)

Show a term.
"""
function Base.show(io::IO, term::Term)
    @printf io "%s{%s%s}(id=%s, value=%s, bondkind=%s, factor=%s)" kind(term) rank(term) statistics(term) term|>id decimaltostr(term.value) term.bondkind decimaltostr(term.factor)
end

"""
    repr(term::Term, bond::AbstractBond, config::IDFConfig) -> String

Get the repr representation of a term on a bond with a given config.
"""
function Base.repr(term::Term, bond::AbstractBond, config::IDFConfig)
    cache = String[]
    if term.bondkind == bond|>kind
        value = term.value * term.amplitude(bond) * term.factor
        if abs(value) ≠ 0
            pids = NTuple{length(bond), pidtype(bond)}(point.pid for point in bond)
            interanls = NTuple{length(bond), valtype(config)}(config[pid] for pid in pids)
            for coupling in values(term.couplings(bond))
                (length(expand(coupling, pids, interanls, term|>kind|>Val)) > 0) &&  push!(cache, @sprintf "%s: %s" abbr(term) repr(value*coupling))
            end
        end
    end
    return join(cache, "\n")
end

"""
    replace(term::Term{ST, K, R, I}; kwargs...) where {ST, K, R, I} -> Term

Replace some attributes of a term with key word arguments.
"""
@generated function Base.replace(term::Term{ST, K, R, I}; kwargs...) where {ST, K, R, I}
    exprs = [:(get(kwargs, $name, getfield(term, $name))) for name in QuoteNode.(term|>fieldnames)]
    return :(Term{ST, K, R, I}($(exprs...)))
end

"""
    +(term::Term) -> Term
    -(term::Term) -> Term
    *(term::Term, factor) -> Term
    *(factor, term::Term) -> Term
    /(term::Term, factor) -> Term

Allowed arithmetic operations for a term.
"""
Base.:+(term::Term) = term
Base.:-(term::Term) = term * (-1)
Base.:*(term::Term, factor) = factor * term
Base.:*(factor, term::Term) = replace(term, factor=factor*term.factor)
Base.:/(term::Term, factor) = term * (1/factor)

"""
    one(term::Term) -> Term

Get a unit term.
"""
Base.one(term::Term) = replace(term, value=one(term.value))

"""
    zero(term::Term) -> Term

Get a zero term.
"""
Base.zero(term::Term) = replace(term, value=zero(term.value))

"""
    expand!(operators::Operators, term::Term, bond::AbstractBond, config::IDFConfig,
            table::Union{Table, Nothing}=nothing,
            half::Bool=false,
            coord::Union{Val{true}, Val{false}}=coordpresent
            ) -> Operators
    expand!(operators::Operators, term::Term, bonds::Bonds, config::IDFConfig,
            table::Union{Table, Nothing}=nothing,
            half::Bool=false,
            coord::Union{Val{true}, Val{false}}=coordpresent
            ) -> Operators

Expand the operators of a term on a bond/set-of-bonds with a given config.

The `half` parameter determines the behavior of generating operators, which falls into the following two categories
* `true`: "Hermitian half" of the generated operators
* `false`: "Hermitian whole" of the generated operators
"""
function expand!(   operators::Operators, term::Term, bond::AbstractBond, config::IDFConfig,
                    table::Union{Table, Nothing}=nothing,
                    half::Bool=false,
                    coord::Union{Val{true}, Val{false}}=coordpresent
                    )
    if term.bondkind == bond|>kind
        value = term.value * term.amplitude(bond) * term.factor
        if abs(value) ≠ 0
            apriori = isHermitian(term)
            ptype = otype(term|>typeof, oidtype(config|>typeof|>valtype|>eltype, bond|>typeof, table|>typeof, coord))
            record = ((apriori === nothing) && (length(operators) > 0)) ? Set{ptype|>idtype}() : nothing
            rtype, itype = fieldtype(ptype|>idtype|>eltype, :rcoord), fieldtype(ptype|>idtype|>eltype, :icoord)
            pids = NTuple{length(bond), pidtype(bond)}(point.pid for point in bond)
            rcoords = NTuple{length(bond), SVector{dimension(bond), Float}}(point.rcoord for point in bond)
            icoords = NTuple{length(bond), SVector{dimension(bond), Float}}(point.icoord for point in bond)
            interanls = NTuple{length(bond), valtype(config)}(config[pid] for pid in pids)
            for coupling in values(term.couplings(bond))
                perm = couplingcenters(typeof(coupling), coupling.id.centers, Val(rank(bond)))::NTuple{rank(ptype), Int}
                orcoords, oicoords = termcoords(rtype, rcoords, perm), termcoords(itype, icoords, perm)
                for (coeff, oindexes) in expand(coupling, pids, interanls, term|>kind|>Val) # needs improvement memory allocation 7 times for each fock coupling
                    isa(table, Table) && any(NTuple{rank(ptype), Bool}(!haskey(table, index) for index in oindexes)) && continue
                    id = ID(OID, oindexes, orcoords, oicoords, termseqs(table, oindexes))
                    if apriori === true
                        add!(operators, rawtype(ptype)(convert(ptype|>valtype, value*coeff/(half ? 2 : 1)), id))
                    elseif apriori == false
                        opt = rawtype(ptype)(convert(ptype|>valtype, value*coeff), id)
                        add!(operators, opt)
                        half || add!(operators, opt')
                    else
                        if !((record === nothing) ? haskey(operators, id') : in(id', record))
                            (record === nothing) || push!(record, id)
                            ovalue = valtype(ptype)(value*coeff/termfactor(nothing, id, term|>kind|>Val)) # needs improvement memory allocation 2 times for each
                            opt = rawtype(ptype)(ovalue, id)
                            add!(operators, opt)
                            half || add!(operators, opt')
                        end
                    end
                end
            end
        end
    end
    return operators
end
@generated function expand!(operators::Operators, term::Term, bonds::Bonds, config::IDFConfig,
                            table::Union{Table, Nothing}=nothing,
                            half::Bool=false,
                            coord::Union{Val{true}, Val{false}}=coordpresent)
    exprs = []
    for i = 1:rank(bonds)
        push!(exprs, :(for bond in bonds.bonds[$i] expand!(operators, term, bond, config, table, half, coord) end))
    end
    push!(exprs, :(return operators))
    return Expr(:block, exprs...)
end
termcoords(::Type{<:Nothing}, rcoords::NTuple{N, SVector{M, Float}}, perm::NTuple{R, Int}) where {N, M, R} = NTuple{R, Nothing}(nothing for i = 1:R)
termcoords(::Type{<:SVector}, rcoords::NTuple{N, SVector{M, Float}}, perm::NTuple{R, Int}) where {N, M, R} = NTuple{R, SVector{M, Float}}(rcoords[p] for p in perm)
termseqs(::Nothing, indexes::NTuple{N, <:Index}) where {N} = NTuple{N, Nothing}(nothing for i = 1:N)
@generated termseqs(table::Table{I}, indexes::NTuple{N, I}) where {N, I<:Index} = Expr(:tuple, [:(table[indexes[$i]]) for i = 1:N]...)
termfactor(::Nothing, id::ID{OID}, ::Val{K}) where {K} = isHermitian(id) ? 2 : 1

"""
    expand(term::Term, bond::AbstractBond, config::IDFConfig, table::Union{Table, Nothing}=nothing, half::Bool=false, coord::Union{Val{true}, Val{false}}=coordpresent) -> Operators
    expand(term::Term, bonds::Bonds, config::IDFConfig, table::Union{Table, Nothing}=nothing, half::Bool=false, coord::Union{Val{true}, Val{false}}=coordpresent) -> Operators

Expand the operators of a term on a bond/set-of-bonds with a given config.
"""
function expand(term::Term, bond::AbstractBond, config::IDFConfig, table::Union{Table, Nothing}=nothing, half::Bool=false, coord::Union{Val{true}, Val{false}}=coordpresent)
    optype = otype(term|>typeof, oidtype(config|>typeof|>valtype|>eltype, bond|>typeof, table|>typeof, coord))
    expand!(Operators{idtype(optype), optype}(), term, bond, config, table, half, coord)
end
function expand(term::Term, bonds::Bonds, config::IDFConfig, table::Union{Table, Nothing}=nothing, half::Bool=false, coord::Union{Val{true}, Val{false}}=coordpresent)
    optype = otype(term|>typeof, oidtype(config|>typeof|>valtype|>eltype, bonds|>eltype, table|>typeof, coord))
    expand!(Operators{idtype(optype), optype}(), term, bonds, config, table, half, coord)
end

"""
    update!(term::Term, args...; kwargs...) -> Term

Update the value of a term by its `modulate` function.
"""
function update!(term::Term, args...; kwargs...)
    value = term.modulate(args...; kwargs...)
    (value === nothing) || (term.value = value)
    return term
end

"""
    Parameters{Names}(values::Number...) where Names

A NamedTuple that contain the key-value pairs.
"""
const Parameters{Names} = NamedContainer{Number, Names}
Parameters{Names}(values::Number...) where {Names} = NamedContainer{Names}(values)

"""
    match(params1::Parameters, params2::Parameters, atol=atol, rtol=rtol) -> Bool

Judge whether the second set of parameters matches the first.
"""
@generated function Base.match(params1::Parameters, params2::Parameters, atol=atol, rtol=rtol)
    names = intersect(fieldnames(params1), fieldnames(params2))
    (length(names) == 0) && return true
    name = QuoteNode(names[1])
    expr = :(isapprox(getfield(params1, $name), getfield(params2, $name), atol=atol, rtol=rtol))
    for i = 2:length(names)
        name = QuoteNode(names[i])
        expr = Expr(:&&, expr, :(isapprox(getfield(params1, $name), getfield(params2, $name), atol=atol, rtol=rtol)))
    end
    return expr
end

"""
    GenOperators(constops::Operators, alterops::NamedContainer{Operators}, boundops::NamedContainer{Operators})
    GenOperators(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}, half::Bool, ::Val{coord}) where coord

A set of operators. This is the core of [`AbstractGenerator`](@ref).
"""
struct GenOperators{C<:Operators, A<:NamedContainer{Operators}, B<:NamedContainer{Operators}}
    constops::C
    alterops::A
    boundops::B
end
@generated function GenOperators(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}, half::Bool, ::Val{coord}) where coord
    @assert isa(coord, Bool) "GenOperators error: not supported coord."
    constterms, alterterms = [], []
    for term in fieldtypes(terms)
        (fieldtype(term, :modulate) <: Nothing) && push!(constterms, term)
        (fieldtype(term, :modulate) <: TermModulate) && push!(alterterms, term)
    end
    names = NTuple{fieldcount(terms), Symbol}(id(term) for term in fieldtypes(terms))
    alternames = NTuple{length(alterterms), Symbol}(id(term) for term in alterterms)
    exprs, alterops, boundops = [], [], []
    push!(exprs, quote
        oidtp = oidtype(config|>valtype|>eltype, bonds|>eltype, table|>typeof, coord|>Val)
        choosedterms = (length($constterms) > 0) ? $constterms : $alterterms
        constoptp = otype(choosedterms[1], oidtp)
        constidtp = constoptp |> idtype
        for i = 2:length(choosedterms)
            tempoptp = otype(choosedterms[i], oidtp)
            constoptp = promote_type(constoptp, tempoptp)
            constidtp = promote_type(constidtp, tempoptp|>idtype)
        end
        constops = Operators{constidtp, constoptp}()
        innerbonds = filter(acrossbonds, bonds, Val(:exclude))
        boundbonds = filter(acrossbonds, bonds, Val(:include))
    end)
    for i = 1:fieldcount(terms)
        push!(boundops, :(expand(one(terms[$i]), boundbonds, config, table, half, coord|>Val)))
        if fieldtype(fieldtype(terms, i), :modulate) <: TermModulate
            push!(alterops, :(expand(one(terms[$i]), innerbonds, config, table, half, coord|>Val)))
        else
            push!(exprs, :(expand!(constops, terms[$i], innerbonds, config, table, half, coord|>Val)))
        end
    end
    alterops = Expr(:tuple, alterops...)
    boundops = Expr(:tuple, boundops...)
    push!(exprs, quote
        alterops = NamedContainer{$alternames}($alterops)
        boundops = NamedContainer{$names}($boundops)
        return GenOperators(constops, alterops, boundops)
    end)
    return Expr(:block, exprs...)
end

"""
    ==(genops1::GenOperators, genops2::GenOperators) -> Bool

Judge whether two sets of operators are equivalent to each other.
"""
Base.:(==)(genops1::GenOperators, genops2::GenOperators) = ==(efficientoperations, genops1, genops2)

"""
    isequal(genops1::GenOperators, genops2::GenOperators) -> Bool

Judge whether two sets of operators are equivalent to each other.
"""
Base.isequal(genops1::GenOperators, genops2::GenOperators) = isequal(efficientoperations, genops1, genops2)

"""
    eltype(ops::GenOperators)
    eltype(::Type{<:GenOperators{S, D, B}}) where {S<:Operators, D<:NamedContainer{Operators}, B<:NamedContainer{Operators}}

Get the eltype of a set of operators, which is defined to be the common operator type of all operators it contains.
"""
Base.eltype(ops::GenOperators) = ops |> typeof |> eltype
@generated function Base.eltype(::Type{<:GenOperators{S, D, B}}) where {S<:Operators, D<:NamedContainer{Operators}, B<:NamedContainer{Operators}}
    optp = S |> valtype
    (fieldcount(D) > 0) && (optp = promote_type(optp, mapreduce(valtype, promote_type, fieldtypes(D))))
    (fieldcount(B) > 0) && (optp = promote_type(optp, mapreduce(valtype, promote_type, fieldtypes(B))))
    return optp
end

"""
    idtype(ops::GenOperators)
    idtype(::Type{<:GenOperators{S, D, B}}) where {S<:Operators, D<:NamedContainer{Operators}, B<:NamedContainer{Operators}}

Get the idtype of a set of operators, which is defined to be the common idtype of all operators it contains.
"""
idtype(ops::GenOperators) = ops |> typeof |> idtype
@generated function idtype(::Type{<:GenOperators{S, D, B}}) where {S<:Operators, D<:NamedContainer{Operators}, B<:NamedContainer{Operators}}
    idtp = S |> keytype
    (fieldcount(D) > 0) && (idtp = promote_type(idtp, mapreduce(keytype, promote_type, fieldtypes(D))))
    (fieldcount(B) > 0) && (idtp = promote_type(idtp, mapreduce(keytype, promote_type, fieldtypes(B))))
    return idtp
end

"""
    empty!(ops::GenOperators) -> GenOperators

Empty a set of operators.
"""
@generated function Base.empty!(ops::GenOperators)
    exprs = [:(empty!(ops.constops))]
    for i = 1:fieldcount(fieldtype(ops, :alterops)) push!(exprs, :(empty!(ops.alterops[$i]))) end
    for i = 1:fieldcount(fieldtype(ops, :boundops)) push!(exprs, :(empty!(ops.boundops[$i]))) end
    push!(exprs, :(return ops))
    return Expr(:block, exprs...)
end

"""
    empty(ops::GenOperators) -> GenOperators

Get an empty copy of a set of operators.
"""
@generated function Base.empty(ops::GenOperators)
    exprs = [:(constops = empty(ops.constops))]
    alterops, boundops = [], []
    for i = 1:fieldcount(fieldtype(ops, :alterops)) push!(alterops, :(empty(ops.alterops[$i]))) end
    for i = 1:fieldcount(fieldtype(ops, :boundops)) push!(boundops, :(empty(ops.boundops[$i]))) end
    alterops = Expr(:tuple, alterops...)
    boundops = Expr(:tuple, boundops...)
    push!(exprs, quote
        alterops = NamedContainer{fieldnames(fieldtype(ops|>typeof, :alterops))}($alterops)
        boundops = NamedContainer{fieldnames(fieldtype(ops|>typeof, :boundops))}($boundops)
        return GenOperators(constops, alterops, boundops)
    end)
    return Expr(:block, exprs...)
end

"""
    expand!(operators::Operators, ops::GenOperators, boundary::Boundary; kwargs...) -> Operators

Expand the operators with the given boundary twist and term coefficients.
"""
@generated function expand!(operators::Operators, ops::GenOperators, boundary::Boundary; kwargs...)
    exprs = [:(add!(operators, ops.constops))]
    for name in QuoteNode.(fieldnames(fieldtype(ops, :alterops)))
        push!(exprs, :(value = get(kwargs, $name, nothing)))
        push!(exprs, :(for opt in values(getfield(ops.alterops, $name)) add!(operators, opt*value) end))
    end
    for name in QuoteNode.(fieldnames(fieldtype(ops, :boundops)))
        push!(exprs, :(value = get(kwargs, $name, nothing)))
        push!(exprs, :(for opt in values(getfield(ops.boundops, $name)) add!(operators, boundary(opt)*value) end))
    end
    push!(exprs, :(return operators))
    return Expr(:block, exprs...)
end

"""
    reset!(genops::GenOperators, terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}, half::Bool, ::Val{coord}) where coord -> GenOperators

Reset a set of operators by new terms, bonds, config, table, etc..
"""
@generated function reset!(genops::GenOperators, terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}, half::Bool, ::Val{coord}) where coord
    exprs = []
    push!(exprs, quote
        empty!(genops)
        innerbonds = filter(acrossbonds, bonds, Val(:exclude))
        boundbonds = filter(acrossbonds, bonds, Val(:include))
    end)
    for (i, term) in enumerate(fieldtypes(terms))
        name = QuoteNode(term|>id)
        push!(exprs, :(expand!(getfield(genops.boundops, $name), one(terms[$i]), boundbonds, config, table, half, coord|>Val)))
        if fieldtype(term, :modulate) <: Nothing
            push!(exprs, :(expand!(genops.constops, terms[$i], innerbonds, config, table, half, coord|>Val)))
        else
            push!(exprs, :(expand!(getfield(genops.alterops, $name), one(terms[$i]), innerbonds, config, table, half, coord|>Val)))
        end
    end
    push!(exprs, :(return genops))
    return Expr(:block, exprs...)
end

"""
    AbstractGenerator{coord, TS<:NamedContainer{Term}, BS<:Bonds, C<:IDFConfig, T<:Union{Nothing, Table}, B<:Boundary, OS<:GenOperators}

Abstract generator.

By protocol, a concrete generator must have the following attributes:
* `terms::TS`: the terms contained in a generator
* `bonds::BS`: the bonds on which the terms are defined
* `config::C`: the configuration of the interanl degrees of freedom
* `table::T`: the index-sequence table, `nothing` for not using such a table
* `half::Bool`: true for generating an Hermitian half of the operators and false for generating the whole
* `boundary::B`: boundary twist for the generated operators, `nothing` for no twist
* `operators::OS`: the generated operators
"""
abstract type AbstractGenerator{coord, TS<:NamedContainer{Term}, BS<:Bonds, C<:IDFConfig, T<:Union{Nothing, Table}, B<:Boundary, OS<:GenOperators} end

"""
    ==(gen1::AbstractGenerator, gen2::AbstractGenerator) -> Bool

Judge whether generators are equivalent to each other.
"""
Base.:(==)(gen1::AbstractGenerator, gen2::AbstractGenerator) = ==(efficientoperations, gen1, gen2)

"""
    isequal(gen1::AbstractGenerator, gen2::AbstractGenerator) -> Bool

Judge whether generators are equivalent to each other.
"""
Base.isequal(gen1::AbstractGenerator, gen2::AbstractGenerator) = isequal(efficientoperations, gen1, gen2)

"""
    Parameters(gen::AbstractGenerator) -> Parameters

Get the parameters of the terms of a generator.
"""
@generated function Parameters(gen::AbstractGenerator)
    names = fieldnames(fieldtype(gen, :terms))
    values = [:(gen.terms[$i].value) for i = 1:fieldcount(fieldtype(gen, :terms))]
    return :(Parameters{$names}($(values...)))
end

"""
    expand!(operators::Operators, gen::AbstractGenerator) -> Operators

Expand the operators of a generator.
"""
expand!(operators::Operators, gen::AbstractGenerator) = expand!(operators, gen.operators, gen.boundary; Parameters(gen)...)

"""
    expand(gen::AbstractGenerator) -> Operators
    expand(gen::AbstractGenerator{coord}, name::Symbol) where coord -> Operators
    expand(gen::AbstractGenerator{coord}, i::Int) where coord -> Operators
    expand(gen::AbstractGenerator{coord}, name::Symbol, i::Int) where coord -> Operators

Expand the operators of a generator:
1) the total operators;
2) the operators of a specific term;
3) the operators on a specific bond;
4) the operators of a specific term on a specific bond.
"""
function expand(gen::AbstractGenerator)
    expand!(Operators{idtype(gen.operators), eltype(gen.operators)}(), gen)
end
function expand(gen::AbstractGenerator{coord}, name::Symbol) where coord
    term = getfield(gen.terms, name)
    optp = otype(term|>typeof, oidtype(gen.config|>valtype|>eltype, gen.bonds|>eltype, gen.table|>typeof, coord))
    result = Operators{idtype(optp), optp}()
    if fieldtype(term|>typeof, :modulate) <: TermModulate
        for opt in getfield(gen.operators.alterops, name)|>values add!(result, opt*term.value) end
        for opt in getfield(gen.operators.boundops, name)|>values add!(result, gen.boundary(opt)*term.value) end
    else
        expand!(result, term, gen.bonds, gen.config, gen.table, gen.half, coord)
    end
    return result
end
@generated function expand(gen::AbstractGenerator{coord}, i::Int) where coord
    exprs = [:(result = Operators{idtype(gen.operators), eltype(gen.operators)}())]
    for i = 1:fieldcount(fieldtype(gen, :terms))
        push!(exprs, :(expand!(result, gen.terms[$i], gen.bonds[i], gen.config, gen.table, gen.half, coord)))
    end
    push!(exprs, :(return result))
    return Expr(:block, exprs...)
end
function expand(gen::AbstractGenerator{coord}, name::Symbol, i::Int) where coord
    expand(getfield(gen.terms, name), gen.bonds[i], gen.config, gen.table, gen.half, coord)
end

"""
    update!(gen::AbstractGenerator; kwargs...) -> typeof(gen)

Update the coefficients of the terms in a generator.
"""
@generated function update!(gen::AbstractGenerator; kwargs...)
    exprs = [:(update!(gen.boundary; kwargs...))]
    for i = 1:fieldcount(fieldtype(gen, :terms))
        (fieldtype(fieldtype(fieldtype(gen, :terms), i), :modulate) <: TermModulate) && push!(exprs, :(update!(gen.terms[$i]; kwargs...)))
    end
    push!(exprs, :(return gen))
    return Expr(:block, exprs...)
end

"""
    Generator(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}=nothing, half::Bool=true, boundary::Boundary=Boundary())
    Generator{coord}(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}=nothing, half::Bool=true, boundary::Boundary=Boundary()) where coord

A generator of operators based on terms, configuration of internal degrees of freedom, table of indices and boundary twist.
"""
struct Generator{coord, TS<:NamedContainer{Term}, BS<:Bonds, C<:IDFConfig, T<:Union{Nothing, Table}, B<:Boundary, OS<:GenOperators} <: AbstractGenerator{coord, TS, BS, C, T, B, OS}
    terms::TS
    bonds::BS
    config::C
    table::T
    half::Bool
    boundary::B
    operators::OS
    function Generator{C}(terms::NamedContainer{Term}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}, half::Bool, boundary::Boundary, operators::GenOperators) where C
        new{C, typeof(terms), typeof(bonds), typeof(config), typeof(table), typeof(boundary), typeof(operators)}(terms, bonds, config, table, half, boundary, operators)
    end
end
function Generator(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}=nothing, half::Bool=true, boundary::Boundary=Boundary())
    Generator{coordpresent}(namedterms(terms), bonds, config, table, half, boundary, GenOperators(terms, bonds, config, table, half, coordpresent))
end
function Generator{coord}(terms::Tuple{Vararg{Term}}, bonds::Bonds, config::IDFConfig, table::Union{Nothing, Table}=nothing, half::Bool=true, boundary::Boundary=Boundary()) where coord
    Generator{coord}(namedterms(terms), bonds, config, table, half, boundary, GenOperators(terms, bonds, config, table, half, coord))
end
@generated function namedterms(terms::Tuple{Vararg{Term}})
    names = NTuple{fieldcount(terms), Symbol}(id(fieldtype(terms, i)) for i = 1:fieldcount(terms))
    return :(NamedContainer{$names}(terms))
end

"""
    empty!(gen::Generator) -> Generator

Empty the :bonds, :config, :table and :operators of a generator.
"""
function Base.empty!(gen::Generator)
    empty!(gen.bonds)
    empty!(gen.config)
    isa(gen.table, Table) && empty!(gen.table)
    empty!(gen.operators)
    return gen
end

"""
    empty(gen::Generator{coord}) where coord -> Generator

Get an empty copy of a generator.
"""
function Base.empty(gen::Generator{coord}) where coord
    Generator{coord}(gen.terms, empty(gen.bonds), empty(gen.config), (gen.table === nothing) ? nothing : empty(gen.table), gen.half, gen.boundary, empty(gen.operators))
end

"""
    reset!(gen::Generator{coord}, lattice::AbstractLattice) where coord -> Generator

Reset a generator by a new lattice.
"""
function reset!(gen::Generator{coord}, lattice::AbstractLattice) where coord
    reset!(gen.bonds, lattice)
    reset!(gen.config, lattice.pids)
    isa(gen.table, Table) && reset!(gen.table, gen.config)
    reset!(gen.operators, Tuple(gen.terms), gen.bonds, gen.config, gen.table, gen.half, coord)
    return gen
end

end  # module


text:
# Benchmark on non-weighted sampling

# require the BenchmarkLite package
using BenchmarkLite
using StatsBase

import StatsBase: direct_sample!, xmultinom_sample!
import StatsBase: knuths_sample!, fisher_yates_sample!, self_avoid_sample!
import StatsBase: seqsample_a!, seqsample_c!, seqsample_d!

### generic sampling benchmarking

mutable struct SampleProc{Alg} <: Proc end

abstract type WithRep end
abstract type NoRep end

mutable struct Direct <: WithRep end
tsample!(s::Direct, a, x) = direct_sample!(a, x)

mutable struct Xmultinom <: WithRep end
tsample!(s::Xmultinom, a, x) = xmultinom_sample!(a, x)

mutable struct Sample_WRep <: WithRep end
tsample!(s::Sample_WRep, a, x) = sample!(a, x; replace=true, ordered=false)

mutable struct Sample_WRep_Ord <: WithRep end
tsample!(s::Sample_WRep_Ord, a, x) = sample!(a, x; replace=true, ordered=true)

mutable struct Knuths <: NoRep end
tsample!(s::Knuths, a, x) = knuths_sample!(a, x)

mutable struct Fisher_Yates <: NoRep end
tsample!(s::Fisher_Yates, a, x) = fisher_yates_sample!(a, x)

mutable struct Self_Avoid <: NoRep end
tsample!(s::Self_Avoid, a, x) = self_avoid_sample!(a, x)

mutable struct Seq_A <: NoRep end
tsample!(s::Seq_A, a, x) = seqsample_a!(a, x)

mutable struct Seq_C <: NoRep end
tsample!(s::Seq_C, a, x) = seqsample_c!(a, x)

mutable struct Seq_D <: NoRep end
tsample!(s::Seq_D, a, x) = seqsample_d!(a, x)

mutable struct Sample_NoRep <: NoRep end
tsample!(s::Sample_NoRep, a, x) = sample!(a, x; replace=false, ordered=false)

mutable struct Sample_NoRep_Ord <: NoRep end
tsample!(s::Sample_NoRep_Ord, a, x) = sample!(a, x; replace=false, ordered=true)

# config is in the form of (n, k)

Base.string(p::SampleProc{Alg}) where {Alg} = lowercase(string(Alg))

Base.length(p::SampleProc, cfg::Tuple{Int,Int}) = cfg[2]
Base.isvalid(p::SampleProc{<:WithRep}, cfg::Tuple{Int,Int}) = ((n, k) = cfg; n >= 1 && k >= 1)
Base.isvalid(p::SampleProc{<:NoRep}, cfg::Tuple{Int,Int}) = ((n, k) = cfg; n >= k >= 1)

Base.start(p::SampleProc, cfg::Tuple{Int,Int}) = Vector{Int}(cfg[2])
Base.run(p::SampleProc{Alg}, cfg::Tuple{Int,Int}, s::Vector{Int}) where {Alg} = tsample!(Alg(), 1:cfg[1], s)
Base.done(p::SampleProc, cfg, s) = nothing

### benchmarking

const ns = 5 * (2 .^ [0:9])
const ks = 2 .^ [1:16]

## with replacement

const procs1 = Proc[ SampleProc{Direct}(),
                     SampleProc{Sample_WRep}(),
                     SampleProc{Xmultinom}(),
                     SampleProc{Sample_WRep_Ord}() ]

const cfgs1 = vec([(n, k) for k in ks, n in ns])

rtable1 = run(procs1, cfgs1; duration=0.2)
println()

## without replacement

const procs2 = Proc[ SampleProc{Knuths}(),
                     SampleProc{Fisher_Yates}(),
                     SampleProc{Self_Avoid}(),
                     SampleProc{Sample_NoRep}(),
                     SampleProc{Seq_A}(),
                     SampleProc{Seq_C}(),
                     SampleProc{Seq_D}(),
                     SampleProc{Sample_NoRep_Ord}() ]

const cfgs2 = (Int, Int)[]
for n in 5 * (2 .^ [0:11]), k in 2 .^ [1:16]
    if k < n
        push!(cfgs2, (n, k))
    end
end

rtable2 = run(procs2, cfgs2; duration=0.2)
println()

## show results

println("Sampling With Replacement")
println("===================================")
show(rtable1; unit=:mps, cfghead="(n, k)")
println()

println("Sampling Without Replacement")
println("===================================")
show(rtable2; unit=:mps, cfghead="(n, k)")
println()

ext: jl


text:
function withdraw(p::Player)
    print_dict(p, "spacer")
    print_dict(p, "withdraw_0")
    s = readline(p.socket)
    amount = tryparse(Int64, s)
    if amount == nothing || amount <= 0
        print_dict(p, "withdraw_1")
        return
    end
    if captcha(p, rand(5:10))
        print_dict(p, "withdraw_4")
    else
        print_dict(p, "withdraw_5")
        if captcha(p, rand(5:10))
            print_dict(p, "withdraw_4")
        else
            print_dict(p, "withdraw_6")
            return
        end
    end

    p.balance += amount

    if p.balance > 10000
        print_dict(p, "withdraw_2")
        p.balance = 10000
    end
end

function captcha(p::Player, num)
    print_dict(p, "withdraw_3")
    print_dict(p, "spacer")
    cards = []
    sum = let sum, suits
        sum = 0
        suits = ["♠", "♥", "♦", "♣"]
        for i in 1:num
            x = rand(1:14)
            sum += x
            if (x == 10) x = "T" end
            if (x == 11) x = "J" end
            if (x == 12) x = "Q" end
            if (x == 13) x = "K" end
            if (x == 14) x = "A" end
            y = rand(suits)
            card = """
            ┌─────────┐
            │$(x) $(y)   $(y)  │
            │         │
            │ $(y)     $(y) │
            │    $(x)    │
            │ $(y)     $(y) │
            │         │
            │  $(y)   $(y) $(x)│
            └─────────┘"""
            card = split(card, "\n")
            push!(cards, card)
        end
        sum
    end
    result = let result, offset
        result = ["","","","","","","","","","","","","","",""]
        offset = rand(1:7)
        for card in cards
            for i in 1:offset-1
                result[i] *= "           "
            end
            for i in offset+9:15
                result[i] *= "           "
            end
            i = offset
            for line in card
                result[i] *= line
                i += 1
            end

            offset = ((offset + rand(-2:2)) % 7)
            if offset < 1
                offset = 1
            end
        end
        join(result, "\n")
    end
    write(p.socket, "$result\n")
    print_dict(p, "spacer")

    s = readline(p.socket)
    return sum == tryparse(Int, s)
end


text:
module Ex5_12c
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.12c - Solution of an EMMP
println("Example 5.12c")

# enter output and fault vector dimensions
p = 3; mf = 3;
# generate random dimensions for system order and input vectors
nu = Int(floor(1+4*rand())); mu = Int(floor(1+4*rand()));
nd = Int(floor(1+4*rand())); md = Int(floor(1+4*rand()));

# define random Gu(s) and Gd(s) with triplex sensor redundancy
# and Gf(s) for triplex sensor faults
Gu = ones(3,1)*rss(nu,1,mu); # enter Gu(s) in state-space form
Gd = ones(3,1)*rss(nd,1,md); # enter Gd(s) in state-space form
Gf = eye(3);                 # enter Gf(s) for sensor faults
atol = 1.e-7;                # tolerance for rank tests

# build model with faults
#sysf = [Gu Gd Gf]
sysf = fdimodset([Gu Gd Gf], c = 1:mu, d = mu.+(1:md), f = (mu+md).+(1:mf))

# enter reference model
Mr = FDFilterIF(dss([ 0 1 -1; -1 0 1; 1 -1 0]), 0, 0, mf)

# solve an exact model-matching problem using EMMSYN
Q, R, info = emmsyn(sysf,Mr; atol); info

# one step solution
# solve Qbar*Ge = Me, where Ge = [Gu Gd Gf; I 0 0] and Me = [0 0 Mr ].
Ge = [sysf.sys; eye(mu,mu+md+mf)]; Me = [zeros(p,mu+md) Mr.sys];
Qbar = glsol(Ge, Me; atol)[1]

# compare solutions
@test iszero(Q.sys-Qbar; atol)

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
mutable struct SSLConfig
    data::Ptr{Cvoid}
    rng
    chain::CRT
    dbg
    cert
    key
    alpn_protos

    function SSLConfig()
        conf = new()
        conf.data = Libc.malloc(1000)  # 360
        ccall((:mbedtls_ssl_config_init, libmbedtls), Cvoid, (Ptr{Cvoid},), conf.data)
        finalizer(conf->begin
            ccall((:mbedtls_ssl_config_free, libmbedtls), Cvoid, (Ptr{Cvoid},), conf.data)
            Libc.free(conf.data)
        end, conf)
        conf
    end
end

Base.show(io::IO, c::SSLConfig) = print(io, "MbedTLS.SSLConfig()")

mutable struct SSLContext <: IO
    data::Ptr{Cvoid}
    datalock::ReentrantLock
    config::SSLConfig
    isopen::Bool
    bio

    function SSLContext()
        ctx = new()
        ctx.data = Libc.malloc(1000)  # 488
        ctx.datalock = ReentrantLock()
        ccall((:mbedtls_ssl_init, libmbedtls), Cvoid, (Ptr{Cvoid},), ctx.data)
        finalizer(ctx->begin
            ccall((:mbedtls_ssl_free, libmbedtls), Cvoid, (Ptr{Cvoid},), ctx.data)
            Libc.free(ctx.data)
        end, ctx)
        ctx
    end
end

macro lockdata(ctx, expr)
    esc(quote
        lock($ctx.datalock)
        @assert $ctx.datalock.reentrancy_cnt == 1
        try
            $expr
        finally
            unlock($ctx.datalock)
        end
    end)
end

function config_defaults!(config::SSLConfig; endpoint=MBEDTLS_SSL_IS_CLIENT,
    transport=MBEDTLS_SSL_TRANSPORT_STREAM, preset=MBEDTLS_SSL_PRESET_DEFAULT)
    @err_check ccall((:mbedtls_ssl_config_defaults, libmbedtls), Cint,
        (Ptr{Cvoid}, Cint, Cint, Cint),
        config.data, endpoint, transport, preset)
end

function authmode!(config::SSLConfig, auth)
    ccall((:mbedtls_ssl_conf_authmode, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Cint),
        config.data, auth)
end

function rng!(config::SSLConfig, f_rng::Ptr{Cvoid}, rng)
    ccall((:mbedtls_ssl_conf_rng, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any),
        config.data, f_rng, rng)
end

function rng!(config::SSLConfig, rng::AbstractRNG)
    config.rng = rng
    rng!(config, c_rng[], rng)
end

function ca_chain!(config::SSLConfig, chain=crt_parse_file(joinpath(dirname(@__FILE__), "../deps/cacert.pem")))
    config.chain = chain
    ccall((:mbedtls_ssl_conf_ca_chain, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        config.data, chain.data, C_NULL)
end

function own_cert!(config::SSLConfig, cert::CRT, key::PKContext)
    config.cert = cert
    config.key = key
    @err_check ccall((:mbedtls_ssl_conf_own_cert, libmbedtls), Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        config.data, cert.data, key.data)
end

function setup!(ctx::SSLContext, conf::SSLConfig)
    @lockdata ctx begin
        ctx.config = conf
        @err_check ccall((:mbedtls_ssl_setup, libmbedtls), Cint,
            (Ptr{Cvoid}, Ptr{Cvoid}),
            ctx.data, conf.data)
    end
end

function set_bio!(ssl_ctx::SSLContext, ctx, f_send::Ptr{Cvoid}, f_recv::Ptr{Cvoid})
    @lockdata ssl_ctx begin
        ccall((:mbedtls_ssl_set_bio, libmbedtls), Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            ssl_ctx.data, ctx, f_send, f_recv, C_NULL)
    end
end

function f_send(c_ctx, c_msg, sz)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    !isopen(jl_ctx.bio) && return Cint(MBEDTLS_ERR_NET_CONN_RESET)
    return Cint(unsafe_write(jl_ctx.bio, c_msg, sz))
end

function f_recv(c_ctx, c_msg, sz)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    n = bytesavailable(jl_ctx.bio)
    if n == 0
        return isopen(jl_ctx.bio) ? Cint(MBEDTLS_ERR_SSL_WANT_READ) :
                    jl_ctx.isopen ? Cint(MBEDTLS_ERR_NET_CONN_RESET) :
                                    Cint(n)
    end
    n = min(sz, n)
    unsafe_read(jl_ctx.bio, c_msg, n)
    return Cint(n)
end

function set_bio!(ssl_ctx::SSLContext, jl_ctx::T) where {T<:IO}
    ssl_ctx.bio = jl_ctx
    set_bio!(ssl_ctx, pointer_from_objref(ssl_ctx), c_send[], c_recv[])
    nothing
end

function dbg!(conf::SSLConfig, f::Ptr{Cvoid}, p)
    ccall((:mbedtls_ssl_conf_dbg, libmbedtls), Cvoid,
        (Ptr{Cvoid}, Ptr{Cvoid}, Any),
        conf.data, f, p)
end

function f_dbg(f, level, filename, number, msg)
    f(level, unsafe_string(filename), number, unsafe_string(msg))
    nothing
end

function dbg!(conf::SSLConfig, f)
    conf.dbg = f
    dbg!(conf, c_dbg[], f)
    nothing
end

@enum(DebugThreshold,
    NONE = 0,
    ERROR,
    STATE_CHANGE,
    INFO,
    VERBOSE)

function set_dbg_level(level)
    ccall((:mbedtls_debug_set_threshold, libmbedtls), Cvoid,
        (Cint,), Cint(level))
    nothing
end

Base.wait(ctx::SSLContext) = (eof(ctx.bio); nothing)
                             # eof blocks if the receive buffer is empty

function handshake(ctx::SSLContext)
    while true
        n = @lockdata ctx begin
            ccall((:mbedtls_ssl_handshake, libmbedtls), Cint,
                  (Ptr{Cvoid},), ctx.data)
        end
        if n == 0
            break
        end
        if n != MBEDTLS_ERR_SSL_WANT_READ
            mbed_err(n)
        end
        wait(ctx)
    end
    ctx.isopen = true

    @async while isopen(ctx)
        # Ensure that libuv is reading data from the socket in case the peer
        # has sent a close_notify message on an otherwise idle connection.
        # https://tools.ietf.org/html/rfc5246#section-7.2.1
        Base.start_reading(ctx.bio)
        try
            wait(ctx.bio.readnotify)
        catch e
            if e isa Base.IOError
                # Ignore read errors (IOError ECONNRESET)
                # https://github.com/JuliaWeb/MbedTLS.jl/issues/148
            else
                rethrow(e)
            end
        end
        yield()
    end

    return
end

function set_alpn!(conf::SSLConfig, protos)
    conf.alpn_protos = protos
    @err_check ccall((:mbedtls_ssl_conf_alpn_protocols, libmbedtls), Cint,
                     (Ptr{Cvoid}, Ptr{Ptr{Cchar}}), conf.data, protos)
    nothing
end

function alpn_proto(ctx::SSLContext)
    rv = ccall((:mbedtls_ssl_get_alpn_protocol, libmbedtls), Ptr{Cchar},
               (Ptr{Cvoid},), ctx.data)
    unsafe_string(rv)
end

function Base.unsafe_write(ctx::SSLContext, msg::Ptr{UInt8}, N::UInt)
    nw = 0
    while nw < N
        ret = @lockdata ctx begin
            ccall((:mbedtls_ssl_write, libmbedtls), Cint,
                  (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                  ctx.data, msg, N - nw)
        end
        ret < 0 && mbed_err(ret)
        nw += ret
        msg += ret
    end
    return Int(nw)
end

Base.write(ctx::SSLContext, msg::UInt8) = write(ctx, Ref(msg))

function Base.unsafe_read(ctx::SSLContext, buf::Ptr{UInt8}, nbytes::UInt; err=true)
    nread::UInt = 0
    while nread < nbytes
        n = @lockdata ctx begin
            ccall((:mbedtls_ssl_read, libmbedtls), Cint,
                   (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                   ctx.data, buf + nread, nbytes - nread)
        end
        if n == MBEDTLS_ERR_SSL_PEER_CLOSE_NOTIFY || n == 0
            close(ctx)
            err ? throw(EOFError()) : return nread
        elseif n == MBEDTLS_ERR_SSL_WANT_READ
            wait(ctx)
        elseif n < 0
            mbed_err(n)
        else
            nread += n
        end
    end
end

Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes=length(buf)) = readbytes!(ctx, buf, UInt(nbytes))
function Base.readbytes!(ctx::SSLContext, buf::Vector{UInt8}, nbytes::UInt)
    nr = unsafe_read(ctx, pointer(buf), nbytes; err=false)
    if nr !== nothing
        resize!(buf, nr::UInt)
    else
        nr = nbytes
    end
    return Int(nr::UInt)
end

Base.readavailable(ctx::SSLContext) = read(ctx, bytesavailable(ctx))

function Base.eof(ctx::SSLContext)
    bytesavailable(ctx)>0 && return false
    return eof(ctx.bio) && bytesavailable(ctx) == 0
end

function Base.close(ctx::SSLContext)
    @lockdata ctx begin
        if isopen(ctx.bio)
            try
                # This is ugly, but a harmless broken pipe exception will be
                # thrown if the peer closes the connection without responding
                ccall((:mbedtls_ssl_close_notify, libmbedtls),
                       Cint, (Ptr{Cvoid},), ctx.data)
            catch
            end
            close(ctx.bio)
        end
        ctx.isopen = false
    end
    nothing
end

function Base.isopen(ctx::SSLContext)

    if !ctx.isopen || !isopen(ctx.bio)
        return false
    end

    decrypt_available_bytes(ctx)

    return ctx.isopen && isopen(ctx.bio)
end

function decrypt_available_bytes(ctx::SSLContext)

    # Zero-byte read causes MbedTLS to call f_recv (always non-blocking)
    # and decrypt any bytes that are already in the LibuvStream read buffer.
    # https://esp32.com/viewtopic.php?t=1101#p4884
    n = @lockdata ctx begin
        ccall((:mbedtls_ssl_read, libmbedtls), Cint,
              (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
                 ctx.data,     C_NULL,       0)
    end
    if n == MBEDTLS_ERR_SSL_PEER_CLOSE_NOTIFY
        close(ctx)
    elseif n == MBEDTLS_ERR_SSL_WANT_READ
        # ignore
    elseif n < 0
        mbed_err(n)
    end
end

function get_peer_cert(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_peer_cert, libmbedtls), Ptr{Cvoid}, (Ptr{Cvoid},), ctx.data)
    return CRT(data)
end

function get_version(ctx::SSLContext)
    if isdefined(ctx, :config)
        data = ccall((:mbedtls_ssl_get_version, libmbedtls), Ptr{UInt8}, (Ptr{Cvoid},), ctx.data)
        return unsafe_string(data)
    else
        throw(ArgumentError("`ctx` hasn't been initialized with an MbedTLS.SSLConfig; run `MbedTLS.setup!(ctx, conf)`"))
    end
end

function get_ciphersuite(ctx::SSLContext)
    data = ccall((:mbedtls_ssl_get_ciphersuite, libmbedtls), Ptr{UInt8}, (Ptr{Cvoid},), ctx.data)
    return unsafe_string(data)
end

function Base.bytesavailable(ctx::SSLContext)

    decrypt_available_bytes(ctx)

    @lockdata ctx begin

        # Now that the bufferd bytes have been processed, find out how many
        # decrypted bytes are available.
        return Int(ccall((:mbedtls_ssl_get_bytes_avail, libmbedtls),
                         Csize_t, (Ptr{Cvoid},), ctx.data))
    end
end

function hostname!(ctx::SSLContext, hostname)
    @err_check ccall((:mbedtls_ssl_set_hostname, libmbedtls), Cint,
      (Ptr{Cvoid}, Cstring), ctx.data, hostname)
end

Sockets.getsockname(ctx::SSLContext) = Sockets.getsockname(ctx.bio)

const c_send = Ref{Ptr{Cvoid}}(C_NULL)
const c_recv = Ref{Ptr{Cvoid}}(C_NULL)
const c_dbg = Ref{Ptr{Cvoid}}(C_NULL)
function __sslinit__()
    c_send[] = @cfunction(f_send, Cint, (Ptr{Cvoid}, Ptr{UInt8}, Csize_t))
    c_recv[] = @cfunction(f_recv, Cint, (Ptr{Cvoid}, Ptr{UInt8}, Csize_t))
    c_dbg[] = @cfunction(f_dbg, Cvoid, (Any, Cint, Ptr{UInt8}, Cint, Ptr{UInt8}))
end


text:
module Models

using Base: Bool, Symbol
using ..FastAI

using BSON
using Flux
using Zygote
using DataDeps
using InlineTest

include("layers.jl")
include("blocks.jl")
include("xresnet.jl")
include("unet.jl")

export xresnet18, xresnet50, UNetDynamic

end


text:
__precompile__()

module EnhancedGJK

using StaticArrays
import GeometryTypes
const gt = GeometryTypes
using CoordinateTransformations: Transformation,
                                 transform_deriv,
                                 IdentityTransformation
using LinearAlgebra
using Statistics: mean
using Meshing
using ForwardDiff
# using DiffResults
# using GeometryBasics: Mesh, TriangleFace, GLTriangleFace, Point3f0

export CollisionCache,
       gjk!,  # signed volume method
       gjk,  # signed volume method
       body_simplex,  # return simplex points on body in body frame
       gjk_original!, # johnson method
       gjk_original,  # johnson method
       GJKResult,
       NeighborMesh,
       ReferenceDistance,
       closest_point_in_world,
       closest_point_in_body,
       separation_distance,
       simplex_penetration_distance,
       linear_combination,
       implicit_surface,
       construct_surface,
       construct_centers,
       eval_surface,
       surface_gradient!,
       fzero,
       make_mesh,
       make_geom,
       make_mesh_sdf

include("tagged_points.jl")
include("simplices.jl")
include("johnson_distance.jl")
include("neighbor_mesh.jl")
include("traits.jl")
include("gjk.jl")
include("reference_distance.jl")
include("surface.jl")

end # module


text:
using Main.COSMO
using Test, LinearAlgebra, SparseArrays, Random

function create_diagonal_extractor(Tv, n::Ti) where {Ti}
    # Creates a matrix A that slices out the diagonal entries Xii of
    # a upper-triangular-vectorized square matrix x=vec(X) of dimension n
    A = spzeros(Tv, Ti, n, Ti(n*(n + 1)/2))
    idx = 0
    for i = 1:n
        idx += i
        A[i, idx] = one(Tv)
    end
    return A
end

function create_closest_correlation_model(C::AbstractMatrix{Tv}, settings=nothing;
    psdcone_type=COSMO.PsdConeTriangle) where Tv
    # Closest Correlation Matrix problem
    # min   1/2 ||X-C||^2
    # s.t.    Xii = 1
    #         X ⪴ 0

    @assert psdcone_type in [COSMO.PsdConeTriangle COSMO.PsdConeTriangleLanczos]
    n = size(C, 1)
    @assert size(C, 2) == size(C, 1)

    # Contraint #1: Diagonal equal to one
    A1 = create_diagonal_extractor(Tv, n)
    b1 = -ones(Tv, n)
    constraints = [COSMO.Constraint(A1, b1, COSMO.ZeroSet)]

    # Contraint #2: decision variable is a PSD
    n2 =  Int(n*(n + 1)/2)
    A2 = SparseMatrixCSC(one(Tv)*I, n2, n2)
    b2 = zeros(Tv, n2)
    push!(constraints, COSMO.Constraint(A2, b2, psdcone_type))

    P = SparseMatrixCSC(Tv(1/2)*I, n2, n2)
    q = -COSMO.extract_upper_triangle(C)

    model = COSMO.Model()
    if settings != nothing
        COSMO.assemble!(model, P, q, constraints, settings=settings)
    else
        COSMO.assemble!(model, P, q, constraints)
    end

    return model
end


text:
# io/z-list.jl

include("caller_name.jl")
include("fld-read.jl")
include("fld-write.jl")
include("ir_dump.jl")
include("prompt.jl")
include("shows.jl")


text:
function condition12(varDel::Int64, varAdd1::Int64, varAdd2::Int64, sol::Solution, prob::Problem)::Bool

    test1 = true
    test2 = true

    for var = 1:prob.nbVar
        if !(prob.listOccur[varAdd1, var] == 0 || sol.x[var] == 0 ||  var == varDel)
            test1 = false
        end
        if !(prob.listOccur[varAdd2, var] == 0 || sol.x[var] == 0 || var == varDel)
            test2 = false
        end
    end

    return test1 && test2 && prob.listOccur[varAdd1, varAdd2] == 0
end

function condition21(varDel1::Int64, varDel2::Int64, varAdd::Int64, sol::Solution, prob::Problem)

    test = true

    for var = 1:prob.nbVar

        if !(prob.listOccur[varAdd, var] == 0 || sol.x[var] == 0 || var == varDel1 || var == varDel2)
            test = false
        end

    end

    return test
end

function condition11(varDel::Int64, varAdd::Int64, sol::Solution, prob::Problem)::Bool

    test = true

    for var = 1:prob.nbVar
        if !(prob.listOccur[varAdd, var] == 0 || sol.x[var] == 0 || var == varDel)
            test = false
        end
    end

    return test
end

function condition01(varAdd::Int64, sol::Solution, prob::Problem)::Bool
    return sol.estLibre[varAdd]
end

function exchange12(sol::Solution, prob::Problem)

    varsMax::Tuple{Int64, Int64, Int64} = (-1, -1, -1)
    zMax = sol.z
    xMax = copy(sol.x)
    estLibreMax = copy(sol.estLibre)

    change::Bool = false

    for varDel in sol.var1

        for varAdd1 in sol.var0
            for varAdd2 in sol.var0
                if condition12(varDel, varAdd1, varAdd2, sol, prob)
                    if sol.z + prob.C[varAdd1] + prob.C[varAdd2] - prob.C[varDel] > zMax

                        zMax = sol.z + prob.C[varAdd1] + prob.C[varAdd2] - prob.C[varDel]
                        varsMax = (varDel, varAdd1, varAdd2)
                        change = true

                        varDel, varAdd1, varAdd2 = varsMax

                        xMax[varDel] = 0
                        xMax[varAdd1] = 1
                        xMax[varAdd2] = 1
                        for var in sol.var0
                            if sol.estLibre[var]
                                estLibreMax[var] = !(prob.listOccur[var, varAdd1] == 1 || prob.listOccur[var, varAdd2] == 1)
                            else
                                if prob.listOccur[var, varDel] == 1
                                    estLibreMax[var] = !(prob.listOccur[var, varAdd1] == 1 || prob.listOccur[var, varAdd2] == 1)
                                end
                            end
                        end
                        estLibreMax[varAdd1] = false
                        estLibreMax[varAdd2] = false
                        estLibreMax[varDel] = !(prob.listOccur[varDel, varAdd1] == 1 || prob.listOccur[varDel, varAdd2] == 1)

                        return Solution(xMax, zMax, estLibreMax), change

                    end
                end
            end
        end
    end
    return sol, false
    # if change
    #     println(varsMax)
    #     varDel, varAdd1, varAdd2 = varsMax
    #
    #     xMax[varDel] = 0
    #     xMax[varAdd1] = 1
    #     xMax[varAdd2] = 1
    #     for var in sol.var0
    #         if sol.estLibre[var]
    #             estLibreMax[var] = !(prob.listOccur[var, varAdd1] == 1 || prob.listOccur[var, varAdd2] == 1)
    #         else
    #             if prob.listOccur[var, varDel] == 1
    #                 estLibreMax[var] = !(prob.listOccur[var, varAdd1] == 1 || prob.listOccur[var, varAdd2] == 1)
    #             end
    #         end
    #     end
    #     estLibreMax[varAdd1] = false
    #     estLibreMax[varAdd2] = false
    #     estLibreMax[varDel] = !(prob.listOccur[varDel, varAdd1] == 1 || prob.listOccur[varDel, varAdd2] == 1)
    # end
    # return Solution(xMax, zMax, estLibreMax), change
end

function exchange21(sol::Solution, prob::Problem)
    varsMax::Tuple{Int64, Int64, Int64} = (-1, -1, -1) # (varDel1, varDel2, varAdd)
    zMax = sol.z
    xMax = copy(sol.x)
    estLibreMax = copy(sol.estLibre)

    change::Bool = false

    for indVarDel1 = 1:length(sol.var1)
        varDel1 = sol.var1[indVarDel1]
        for indVarDel2 = indVarDel1:length(sol.var1)
            varDel2 = sol.var1[indVarDel2]

            for indVarAdd = 1:length(sol.var0)
                varAdd = sol.var0[indVarAdd]
                if condition21(varDel1, varDel2, varAdd, sol, prob)
                    if sol.z + prob.C[varAdd] - prob.C[varDel1] - prob.C[varDel2] > zMax

                        zMax = sol.z + prob.C[varAdd] - prob.C[varDel1] - prob.C[varDel2]
                        varsMax = (varDel1, varDel2, varAdd)
                        change = true

                        xMax[varDel1] = 0
                        xMax[varDel2] = 0
                        xMax[varAdd] = 1
                        for var in sol.var0
                            if sol.estLibre[var]
                                estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1 )
                            else
                                if prob.listOccur[var, varDel1] == 1 || prob.listOccur[var, varDel2] == 1
                                    estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1 )
                                end
                            end
                        end
                        estLibreMax[varAdd] = false
                        estLibreMax[varDel1] = !(prob.listOccur[varDel1, varAdd] == 1 )
                        estLibreMax[varDel2] = !(prob.listOccur[varDel2, varAdd] == 1 )

                        return Solution(xMax, zMax, estLibreMax), change

                    end
                end
            end
        end
    end
    return sol, false
    # if change
    #     varDel1, varDel2, varAdd = varsMax
    #
    #     xMax[varDel1] = 0
    #     xMax[varDel2] = 0
    #     xMax[varAdd] = 1
    #     for var in sol.var0
    #         if sol.estLibre[var]
    #             estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1 )
    #         else
    #             if prob.listOccur[var, varDel1] == 1 || prob.listOccur[var, varDel2] == 1
    #                 estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1 )
    #             end
    #         end
    #     end
    #     estLibreMax[varAdd] = false
    #     estLibreMax[varDel1] = !(prob.listOccur[varDel1, varAdd] == 1 )
    #     estLibreMax[varDel2] = !(prob.listOccur[varDel2, varAdd] == 1 )
    # end
    # return Solution(xMax, zMax, estLibreMax), change
end

function exchange11(sol::Solution, prob::Problem)
    varsMax::Tuple{Int64, Int64} = (-1, -1)
    zMax = sol.z
    xMax = copy(sol.x)
    estLibreMax = copy(sol.estLibre)

    change::Bool = false

    for varDel in sol.var1

        for varAdd in sol.var0
            if condition11(varDel, varAdd, sol, prob)
                if sol.z + prob.C[varAdd] - prob.C[varDel] > zMax

                    zMax = sol.z + prob.C[varAdd] - prob.C[varDel]
                    varsMax = (varDel, varAdd)
                    change = true

                    xMax[varDel] = 0
                    xMax[varAdd] = 1
                    for var in sol.var0
                        if sol.estLibre[var]
                            estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1)
                        else
                            if prob.listOccur[var, varDel] == 1
                                estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1)
                            end
                        end
                    end
                    estLibreMax[varAdd] = false
                    estLibreMax[varDel] = !(prob.listOccur[varDel, varAdd] == 1)

                    return Solution(xMax, zMax, estLibreMax), change
                end
            end
        end
    end
    return sol, false
    # if change
    #     varDel, varAdd = varsMax
    #
    #     xMax[varDel] = 0
    #     xMax[varAdd] = 1
    #     for var in sol.var0
    #         if sol.estLibre[var]
    #             estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1)
    #         else
    #             if prob.listOccur[var, varDel] == 1
    #                 estLibreMax[var] = !(prob.listOccur[var, varAdd] == 1)
    #             end
    #         end
    #     end
    #     estLibreMax[varAdd] = false
    #     estLibreMax[varDel] = !(prob.listOccur[varDel, varAdd] == 1)
    # end
    # return Solution(xMax, zMax, estLibreMax), change
end

function exchange01(sol::Solution, prob::Problem)

    varMax::Int64 = -1
    zMax = sol.z
    xMax = copy(sol.x)
    estLibreMax = copy(sol.estLibre)

    change::Bool = false

    for varAdd in sol.var0
        if condition01(varAdd, sol, prob)
            if sol.z + prob.C[varAdd] > zMax

                zMax = sol.z + prob.C[varAdd]
                varMax = varAdd
                change = true

                xMax[varAdd] = 1
                for var in sol.var0
                    if sol.estLibre[var]
                        estLibreMax[var] = (prob.listOccur[var, varAdd] == 0)
                    end
                end

                estLibreMax[varAdd] = false
                return Solution(xMax, zMax, estLibreMax), change
            end
        end
    end
    return sol, false

    # if change
    #     varAdd = varMax
    #
    #     xMax[varAdd] = 1
    #     for var in sol.var0
    #         if sol.estLibre[var]
    #             estLibreMax[var] = (prob.listOccur[var, varAdd] == 0)
    #         end
    #     end
    #
    #     estLibreMax[varAdd] = false
    # end
    # return Solution(xMax, zMax, estLibreMax), change
end


text:
# A set of simple mesh makers for testing
export simple_line_mesh, simple_quad_mesh, simple_hex_mesh

#=
# Builds a 1D interval mesh.
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_line_mesh(nx, bn, interval)
    # mesh only
    Nv = nx;                    # number of vertex nodes
    xv = zeros(1,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = nx-1;                 # number of elements
    el = zeros(Int, 2, nel);    # element vertex maps
    etypes = ones(Int, nel);    # element types (gmsh number)
    numvert = 2*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    f2n = zeros(Int, 1,Nv);     # face2node mapping
    f2e = zeros(Int, 2,Nv);     # face2element mapping
    e2f = zeros(Int, 2,nel);    # element2face mapping
    normals = ones(1,Nv);      # normals of faces
    bdryID = zeros(Int, Nv);    # BID of faces (0=interior)

    scale = interval[2]-interval[1];
    h = scale/(nx-1);               # uniformly divided

    # vertex nodes are initially ordered lexicographically (that's a 17 letter word!)
    for j=1:nx
        xv[1,j] = interval[1] + (j-1)*h;
    end

    # Elements are ordered the same way
    for ei=1:nel
        el[1,ei] = ei;
        el[2,ei] = ei+1;
        e2f[1,ei] = ei;
        e2f[2,ei] = ei+1;
    end

    # face2node, face2element, normals
    for fi=1:Nv
        f2n[1,fi] = fi;
        f2e[1,fi] = fi-1;
        f2e[2,fi] = fi;
        # normals[fi] = 1; # already set to 1
        # bdryID[fi] = 0; # already set to 0
    end
    f2e[2,Nv] = 0;
    bdryID[1] = 1;
    f2e[:,1] = [1,0];
    normals[1] = -1;
    if bn == 2
        bdryID[Nv] = 2; # for two BIDs
    else
        bdryID[Nv] = 1; # for one BID
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    (refel, grid) = grid_from_mesh(mesh);

    # if bn == 2
    #     add_boundary_ID_to_grid(2, x -> x >= interval[2], grid);
    # end

    return (mesh, refel, grid);
end

#=
# Builds a 2D quad mesh
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_quad_mesh(nxy, bn, interval)
    if length(nxy) == 2
        nx = nxy[1];
        ny = nxy[2];
    else
        nx = nxy;
        ny = nx;
    end
    if length(interval) == 2
        interval = [interval[1], interval[2], interval[1], interval[2]];
    end

    # mesh
    Nv = nx*ny;                 # number of vertex nodes
    xv = zeros(2,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = (nx-1)*(ny-1);        # number of elements
    el = zeros(Int, 4, nel);    # element vertex maps
    etypes = 3*ones(Int, nel);  # element types (gmsh number)
    numvert = 4*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    Nf = nx*(ny-1) + ny*(nx-1); # number of faces
    f2n = zeros(Int, 2,Nf);     # face2node mapping
    f2e = zeros(Int, 2,Nf);     # face2element mapping
    e2f = zeros(Int, 4,nel);    # element2face mapping
    normals = ones(2,Nf);       # normals of faces
    bdryID = zeros(Int, Nf);    # BID of faces (0=interior)

    # if bn == 4
    #     bids = [1,2,3,4]; # x=0, x=1, y=0, y=1
    #     allbids = [1,2,3,4];
    # elseif bn == 3
    #     bids = [1,2,3]; # x=0 , x=1, y=0,1
    #     allbids = [1,2,3,3];
    # elseif bn == 2
    #     bids = [1,2]; # x=0,1 , y=0,1
    #     allbids = [1,1,2,2];
    # else
    #     bids = [1]; # everywhere
    #     allbids = [1,1,1,1];
    # end
    bids = [1]; # everywhere
    allbids = [1,1,1,1];

    scalex = interval[2]-interval[1];
    hx = scalex/(nx-1);         # uniformly divided
    scaley = interval[4]-interval[3];
    hy = scaley/(ny-1);         # uniformly divided

    # vertex nodes are ordered lexicographically
    for j=1:ny
        for i=1:nx
            k = i + (j-1)*nx;
            xv[1, k] = interval[1] + (i-1)*hx;
            xv[2, k] = interval[3] + (j-1)*hy;
        end
    end

    # Elements are ordered the same way
    for j=1:(ny-1)
        for i=1:(nx-1)
            ei = i + (j-1)*(nx-1); # element index

            el[1,ei] = i + (j-1)*nx;
            el[2,ei] = i + (j-1)*nx + 1;
            el[3,ei] = i + (j)*nx + 1;
            el[4,ei] = i + (j)*nx;

            # face2node, face2element, element2face, normals
            f1 = 2*(ei-1)+1 + (j-1); # left face index
            f2 = f1+1; # bottom
            f3 = f1+2; # right
            f4 = j==(ny-1) ? Nf-(nx-1)+i : j*2*(nx-1) + j + i*2; # top

            f2n[1,f1] = el[1,ei];
            f2n[2,f1] = el[4,ei];
            f2n[1,f2] = el[1,ei];
            f2n[2,f2] = el[2,ei];
            f2n[1,f3] = el[2,ei];
            f2n[2,f3] = el[3,ei];
            f2n[1,f4] = el[4,ei];
            f2n[2,f4] = el[3,ei];

            f2e[2,f1] = ei;
            f2e[2,f2] = ei;
            f2e[1,f3] = ei;
            f2e[1,f4] = ei;

            e2f[1,ei] = f1;
            e2f[2,ei] = f2;
            e2f[3,ei] = f3;
            e2f[4,ei] = f4;

            normals[:,f1] = [1,0];
            normals[:,f2] = [0,1];
            normals[:,f3] = [1,0];
            normals[:,f4] = [0,1];
        end
    end

    # boundaries
    for j=1:(ny-1)
        eleft = (j-1)*(nx-1) + 1;
        eright = j*(nx-1);

        fleft = e2f[1,eleft];
        fright = e2f[3,eright];

        bdryID[fleft] = 1; # always 1
        bdryID[fright] = allbids[2];

        # need to change normals, f2e for left side
        normals[:,fleft] = [-1,0];
        f2e[:,fleft] = [eleft, 0];
    end
    for i=1:(nx-1)
        ebottom = i;
        etop = nel - (nx-1) + i;

        fbottom = e2f[2,ebottom];
        ftop = e2f[4,etop];

        bdryID[fbottom] = allbids[3];
        bdryID[ftop] = allbids[4];

        # need to change normals, f2e for bottom side
        normals[:,fbottom] = [0,-1];
        f2e[:,fbottom] = [ebottom, 0];
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    (refel, grid) = grid_from_mesh(mesh);

    if bn == 2
        add_boundary_ID_to_grid(2, (x,y) -> (y <= interval[3]) || (y >= interval[4]), grid);
    elseif bn == 3
        add_boundary_ID_to_grid(2, (x,y) -> (x >= interval[2]), grid);
        add_boundary_ID_to_grid(3, (x,y) -> ((y <= interval[3]) || (y >= interval[4])) && (x > interval[1] && x < interval[2]), grid);
    elseif bn == 4
        add_boundary_ID_to_grid(2, (x,y) -> (x >= interval[2]), grid);
        add_boundary_ID_to_grid(3, (x,y) -> (y <= interval[3] && (x > interval[1] && x < interval[2])), grid);
        add_boundary_ID_to_grid(4, (x,y) -> (y >= interval[4] && (x > interval[1] && x < interval[2])), grid);
    end

    return (mesh, refel, grid);
end

#=
# Builds a 3D hex mesh
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_hex_mesh(nxyz, bn, interval)
    if length(nxyz) == 3
        nx = nxyz[1];
        ny = nxyz[2];
        nz = nxyz[3]
    else
        nx = nxyz;
        ny = nx;
        nz = nx;
    end
    if length(interval) == 2
        interval = [interval[1], interval[2], interval[1], interval[2], interval[1], interval[2]];
    end

    # mesh
    Nv = nx*ny*nz;              # number of vertex nodes
    xv = zeros(3,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = (nx-1)*(ny-1)*(nz-1); # number of elements
    el = zeros(Int, 8, nel);    # element vertex maps
    etypes = 5*ones(Int, nel);  # element types (gmsh number)
    numvert = 8*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    Nf = nz*(nx-1)*(ny-1) + ny*(nx-1)*(nz-1) + nx*(nz-1)*(ny-1); # number of faces
    f2n = zeros(Int, 4,Nf);     # face2node mapping
    f2e = zeros(Int, 2,Nf);     # face2element mapping
    e2f = zeros(Int, 6,nel);    # element2face mapping
    normals = ones(3,Nf);       # normals of faces
    bdryID = zeros(Int, Nf);    # BID of faces (0=interior)

    # if bn == 6
    #     bids = [1,2,3,4,5,6]; # all separate
    #     allbids = [1,2,3,4,5,6];
    # elseif bn == 5
    #     bids = [1,2,3,4,5]; # combine z
    #     allbids = [1,2,3,4,5,5];
    # elseif bn == 4
    #     bids = [1,2,3,4]; # combine y and z
    #     allbids = [1,2,3,3,4,4];
    # elseif bn == 3
    #     bids = [1,2,3]; # combine x,y,z
    #     allbids = [1,1,2,2,3,3];
    # elseif bn == 2
    #     bids = [1,2]; # x=0, other
    #     allbids = [1,2,2,2,2,2];
    # else
    #     bids = [1]; # everywhere
    #     allbids = [1,1,1,1,1,1];
    # end

    bids = [1]; # everywhere
    allbids = [1,1,1,1,1,1];

    scalex = interval[2]-interval[1];
    hx = scalex/(nx-1);         # uniformly divided
    scaley = interval[4]-interval[3];
    hy = scaley/(ny-1);         # uniformly divided
    scalez = interval[6]-interval[5];
    hz = scalez/(nz-1);         # uniformly divided

    # vertex nodes are ordered lexicographically
    for k=1:nz
        for j=1:ny
            for i=1:nx
                ni = i + (j-1)*nx + (k-1)*nx*ny;
                xv[1, ni] = interval[1] + (i-1)*hx;
                xv[2, ni] = interval[3] + (j-1)*hy;
                xv[3, ni] = interval[5] + (k-1)*hz;
            end
        end
    end

    # Elements are ordered lexicographically
    # Elemental vertices are ordered according to GMSH order
    nextface = 1;
    for k=1:(nz-1)
        for j=1:(ny-1)
            for i=1:(nx-1)
                ei = i + (j-1)*(nx-1) + (k-1)*(ny-1)*(nx-1); # element index

                el[1,ei] = i + (j-1)*nx + (k-1)*nx*ny;
                el[2,ei] = i + (j-1)*nx + (k-1)*nx*ny + 1;
                el[3,ei] = i + (j)*nx + (k-1)*nx*ny + 1; #GMSH
                el[4,ei] = i + (j)*nx + (k-1)*nx*ny; #GMSH
                el[5,ei] = i + (j-1)*nx + (k)*nx*ny;
                el[6,ei] = i + (j-1)*nx + (k)*nx*ny + 1;
                el[7,ei] = i + (j)*nx + (k)*nx*ny + 1; #GMSH
                el[8,ei] = i + (j)*nx + (k)*nx*ny; #GMSH

                # face2node, face2element, element2face, normals
                if i==1
                    f1 = nextface;
                    nextface = nextface+1;
                else
                    f1 = e2f[2,ei-1];# left face index
                end
                f2 = nextface; # right
                nextface = nextface+1;

                if j==1
                    f3 = nextface;
                    nextface = nextface+1;
                else
                    f3 = e2f[4,ei-(nx-1)]; # bottom
                end
                f4 = nextface; # top
                nextface = nextface+1;

                if k==1
                    f5 = nextface;
                    nextface = nextface+1;
                else
                    f5 = e2f[6,ei-((nx-1)*(ny-1))]; # front
                end
                f6 = nextface; # back
                nextface = nextface+1;

                f2n[:,f1] = el[[1,4,5,8],ei]; # modified to match GMSH
                f2n[:,f2] = el[[2,3,6,7],ei];
                f2n[:,f3] = el[[1,2,5,6],ei];
                f2n[:,f4] = el[[4,3,7,8],ei];
                f2n[:,f5] = el[[1,2,3,4],ei];
                f2n[:,f6] = el[[5,6,7,8],ei];

                f2e[2,f1] = ei;
                f2e[1,f2] = ei;
                f2e[2,f3] = ei;
                f2e[1,f4] = ei;
                f2e[2,f5] = ei;
                f2e[1,f6] = ei;

                e2f[1,ei] = f1;
                e2f[2,ei] = f2;
                e2f[3,ei] = f3;
                e2f[4,ei] = f4;
                e2f[5,ei] = f5;
                e2f[6,ei] = f6;

                normals[:,f1] = [1,0,0];
                normals[:,f2] = [1,0,0];
                normals[:,f3] = [0,1,0];
                normals[:,f4] = [0,1,0];
                normals[:,f5] = [0,0,1];
                normals[:,f6] = [0,0,1];
            end
        end
    end

    # boundaries
    for k=1:(nz-1)
        for j=1:(ny-1)
            eleft = (k-1)*(nx-1)*(ny-1) + (j-1)*(nx-1) + 1;
            eright = (k-1)*(nx-1)*(ny-1) + j*(nx-1);

            fleft = e2f[1,eleft];
            fright = e2f[2,eright];

            bdryID[fleft] = 1; # always 1
            bdryID[fright] = allbids[2];

            # need to change normals, f2e for left side
            normals[:,fleft] = [-1,0,0];
            f2e[:,fleft] = [eleft, 0];
        end
    end

    for k=1:(nz-1)
        for i=1:(nx-1)
            ebottom = (k-1)*(nx-1)*(ny-1) + i;
            etop = (k-1)*(nx-1)*(ny-1) + (nx-1)*(ny-2) + i;

            fbottom = e2f[3,ebottom];
            ftop = e2f[4,etop];

            bdryID[fbottom] = allbids[3];
            bdryID[ftop] = allbids[4];

            # need to change normals, f2e for bottom side
            normals[:,fbottom] = [0,-1, 0];
            f2e[:,fbottom] = [ebottom, 0];
        end
    end

    for j=1:(ny-1)
        for i=1:(nx-1)
            efront = (j-1)*(nx-1) + i;
            eback = (nz-2)*(nx-1)*(ny-1) + efront;

            ffront = e2f[5,efront];
            fback = e2f[6,eback];

            bdryID[ffront] = allbids[5];
            bdryID[fback] = allbids[6];

            # need to change normals, f2e for bottom side
            normals[:,ffront] = [0,0,-1];
            f2e[:,ffront] = [efront, 0];
        end
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    (refel, grid) = grid_from_mesh(mesh);

    tiny = 1e-13;
    if bn == 6
        # bids = [1,2,3,4,5,6]; # all separate
        add_boundary_ID_to_grid(2, (x,y,z) -> (x >= interval[2]-tiny), grid);
        add_boundary_ID_to_grid(3, (x,y,z) -> (y <= interval[3]+tiny), grid);
        add_boundary_ID_to_grid(4, (x,y,z) -> (y >= interval[4]-tiny), grid);
        add_boundary_ID_to_grid(5, (x,y,z) -> (z <= interval[5]+tiny), grid);
        add_boundary_ID_to_grid(6, (x,y,z) -> (z >= interval[6]-tiny), grid);
    elseif bn == 5
        # bids = [1,2,3,4,5]; # combine z
        add_boundary_ID_to_grid(2, (x,y,z) -> (x >= interval[2]-tiny), grid);
        add_boundary_ID_to_grid(3, (x,y,z) -> (y <= interval[3]+tiny), grid);
        add_boundary_ID_to_grid(4, (x,y,z) -> (y >= interval[4]-tiny), grid);
        add_boundary_ID_to_grid(5, (x,y,z) -> (z <= interval[5]+tiny) || (z >= interval[6]-tiny), grid);
    elseif bn == 4
        # bids = [1,2,3,4]; # combine y and z
        add_boundary_ID_to_grid(2, (x,y,z) -> (x >= interval[2]), grid);
        add_boundary_ID_to_grid(3, (x,y,z) -> ((y <= interval[3]+tiny) || (y >= interval[4]-tiny)), grid);
        add_boundary_ID_to_grid(4, (x,y,z) -> ((z <= interval[5]+tiny) || (z >= interval[6]-tiny)), grid);
    elseif bn == 3
        # bids = [1,2,3]; # combine x,y,z
        add_boundary_ID_to_grid(2, (x,y,z) -> (y <= interval[3]+tiny) || (y >= interval[4]-tiny), grid);
        add_boundary_ID_to_grid(3, (x,y,z) -> (z <= interval[5]+tiny) || (z >= interval[6]-tiny), grid);
    elseif bn == 2
        # bids = [1,2]; # x=0, other
        add_boundary_ID_to_grid(2, (x,y,z) -> (x >= interval[2]-tiny) || (y <= interval[3]+tiny) || (y >= interval[4]-tiny) || (z <= interval[5]+tiny) || (z >= interval[6]-tiny), grid);
    end

    return (mesh, refel, grid)
end