text:
module QuantumOperators

using Printf: @printf, @sprintf
using ...Prerequisites: atol, rtol
using ...Prerequisites.NamedVectors: NamedVector
using ...Prerequisites.Traits: efficientoperations, getcontent, contentorder
using ...Prerequisites.Traits: rawtype, fulltype, parametertype, parameterpairs, reparameter, promoteparameters
using ...Prerequisites.CompositeStructures: CompositeDict

import ...Interfaces: id, value, rank, add!, sub!, mul!, div!, ⊗, ⋅, permute
import ...Prerequisites.Traits: contentnames, dissolve, isparameterbound, parameternames

export SingularID, ID, QuantumOperator, OperatorProd, OperatorSum, Scalar, idtype, sequence
export Transformation, Identity, Numericalization, MatrixRepresentation, matrix, matrix!

"""
    SingularID <: NamedVector

A singular id is the building block of the id system of quantum operators.
"""
abstract type SingularID <: NamedVector end

"""
    ID{I<:SingularID, N}

The (composite) id system of quantum operators.

Type alias for `NTuple{N, I} where {N, I<:SingularID}`.
"""
const ID{I<:SingularID, N} = NTuple{N, I}
@inline Base.show(io::IO, cid::Tuple{SingularID, Vararg{SingularID}}) = @printf io "ID(%s)" join(cid, ", ")

"""
    ID(ids::SingularID...)
    ID(ids::NTuple{N, SingularID}) where N

Get the composite id from singular ids.
"""
@inline ID(ids::SingularID...) = ids
@inline ID(ids::NTuple{N, SingularID}) where {N} = ids

"""
    ID(::Type{SID}, attrs::Vararg{NTuple{N}, M}) where {SID<:SingularID, N, M}

Get the composite id from the components of singular ids.
"""
@inline @generated function ID(::Type{SID}, attrs::Vararg{NTuple{N, Any}, M}) where {SID<:SingularID, N, M}
    exprs = []
    for i = 1:N
        args = [:(attrs[$j][$i]) for j = 1:M]
        push!(exprs, :(SID($(args...))))
    end
    return :(ID($(exprs...)))
end

"""
    propertynames(::Type{I}) where I<:ID{SingularID} -> Tuple{Vararg{Symbol}}

Get the property names of a composite id.
"""
@inline @generated function Base.propertynames(I::ID{SingularID})
    exprs = [QuoteNode(Symbol(name, 's')) for name in I|>eltype|>fieldnames]
    return Expr(:tuple, exprs...)
end

"""
    getproperty(cid::ID{SingularID}, name::Symbol)

Get the property of a composite id.
"""
@inline Base.getproperty(cid::ID{SingularID}, name::Symbol) = idgetproperty(cid, Val(name), Val(cid|>propertynames))
@inline @generated function idgetproperty(cid::ID{SingularID}, ::Val{name}, ::Val{names}) where {name, names}
    index = findfirst(isequal(name), names)::Int
    exprs = [:(getfield(cid[$i], $index)) for i = 1:fieldcount(cid)]
    return Expr(:tuple, exprs...)
end

"""
    promote_type(::Type{Tuple{}}, I::Type{<:ID{SingularID, N}}) where N
    promote_type(I::Type{<:ID{SingularID, N}}, ::Type{Tuple{}}) where N

Define the promote rule for ID types.
"""
@inline Base.promote_type(::Type{Tuple{}}, I::Type{<:Tuple{SingularID, Vararg{SingularID}}}) = ID{I|>eltype}
@inline Base.promote_type(I::Type{<:Tuple{SingularID, Vararg{SingularID}}}, ::Type{Tuple{}}) = ID{I|>eltype}

"""
    isless(::Type{<:SingularID}, cid₁::ID{SingularID}, cid₂::ID{SingularID}) -> Bool

Compare two ids and judge whether the first is less than the second.

The comparison rule are as follows:
1. ids with smaller ranks are always less than those with higher ranks;
2. if two ids are of the same rank, the comparison goes just like that between tuples.
"""
@inline function Base.isless(::Type{<:SingularID}, cid₁::ID{SingularID}, cid₂::ID{SingularID})
    r₁, r₂ = cid₁|>rank, cid₂|>rank
    (r₁ < r₂) ? true : (r₁ > r₂) ? false : isless(cid₁, cid₂)
end

"""
    rank(id::ID{SingularID}) -> Int
    rank(::Type{<:ID{SingularID}}) -> Any
    rank(::Type{<:ID{SingularID, N}}) where N -> Int

Get the rank of a composite id.
"""
@inline rank(id::ID{SingularID}) = rank(typeof(id))
@inline rank(::Type{<:ID{SingularID}}) = Any
@inline rank(::Type{<:ID{SingularID, N}}) where N = N

"""
    *(sid₁::SingularID, sid₂::SingularID) -> ID{SingularID}
    *(sid::SingularID, cid::ID{SingularID}) -> ID{SingularID}
    *(cid::ID{SingularID}, sid::SingularID) -> ID{SingularID}
    *(cid₁::ID{SingularID}, cid₂::ID{SingularID}) -> ID{SingularID}

Get the product of the id system.
"""
@inline Base.:*(sid₁::SingularID, sid₂::SingularID) = ID(sid₁, sid₂)
@inline Base.:*(sid::SingularID, cid::ID{SingularID}) = ID(sid, cid...)
@inline Base.:*(cid::ID{SingularID}, sid::SingularID) = ID(cid..., sid)
@inline Base.:*(cid₁::ID{SingularID}, cid₂::ID{SingularID}) = ID(cid₁..., cid₂...)

"""
    QuantumOperator{I, V} <:CompositeDict{I, V}

The abstract type of any quantum operator.
"""
abstract type QuantumOperator{I, V} <:CompositeDict{I, V} end
@inline Base.:(==)(m₁::QuantumOperator, m₂::QuantumOperator) = ==(efficientoperations, m₁, m₂)
@inline Base.isequal(m₁::QuantumOperator, m₂::QuantumOperator) = isequal(efficientoperations, m₁, m₂)
@inline parameternames(::Type{<:QuantumOperator}) = (:id, :value)
@inline isparameterbound(::Type{<:QuantumOperator}, ::Val{:id}, ::Type{I}) where {I<:ID{SingularID}} = !isconcretetype(I)
@inline isparameterbound(::Type{<:QuantumOperator}, ::Val{:value}, ::Type{V}) where V = !isconcretetype(V)

"""
    idtype(m::QuantumOperator)
    idtype(::Type{T}) where {T<:QuantumOperator}

The type of the id of a `QuantumOperator`.
"""
@inline idtype(m::QuantumOperator) = idtype(typeof(m))
@inline @generated idtype(::Type{T}) where {T<:QuantumOperator} = parametertype(supertype(T, :QuantumOperator), :id)

"""
    valtype(m::QuantumOperator)
    valtype(::Type{T}) where {T<:QuantumOperator}

Get the type of the value of a `QuantumOperator`.
"""
@inline Base.valtype(m::QuantumOperator) = valtype(typeof(m))
@inline @generated Base.valtype(::Type{T}) where {T<:QuantumOperator} = parametertype(supertype(T, :QuantumOperator), :value)

"""
    rank(m::QuantumOperator) -> Int
    rank(::Type{M}) where M<:QuantumOperator -> Int

Get the rank of a `QuantumOperator`.
"""
@inline rank(m::QuantumOperator) = rank(typeof(m))
@inline rank(::Type{M}) where {M<:QuantumOperator} = rank(idtype(M))

"""
    OperatorProd{V, I<:ID{SingularID}} <: QuantumOperator{I, V}

The entity that represent the product of several quantum operators.

Basically, a concrete subtype should contain two predefined contents:
- `value::V`: the coefficient of the product
- `id::I`: the total id of the product
"""
abstract type OperatorProd{V, I<:ID{SingularID}} <: QuantumOperator{I, V} end
@inline contentnames(::Type{<:OperatorProd}) = (:value, :id)
@inline parameternames(::Type{<:OperatorProd}) = (:value, :id)
@inline isparameterbound(::Type{<:OperatorProd}, ::Val{:value}, ::Type{V}) where V = false

@inline newvalue(m::OperatorProd, v::Number) = v
"""
    replace(m::OperatorProd, v::Number) -> OperatorProd

Replace the value of an `OperatorProd`.
"""
@inline Base.replace(m::OperatorProd, v::Number) = rawtype(typeof(m))(dissolve(m, newvalue, (v,))...)
@inline dissolve(m::OperatorProd, ::Val{:value}, ::typeof(newvalue), args::Tuple, kwargs::NamedTuple) = newvalue(m, args...; kwargs...)

"""
    value(m::OperatorProd) -> valtype(m)

Get the value of an `OperatorProd`.
"""
@inline value(m::OperatorProd) = getcontent(m, :value)

"""
    id(m::OperatorProd) -> idtype(m)

Get the id of an `OperatorProd`.
"""
@inline id(m::OperatorProd) = getcontent(m, :id)

"""
    Scalar{V}

Scalar quantum operator.
"""
const Scalar{V} = OperatorProd{V, Tuple{}}

"""
    isapprox(m₁::OperatorProd, m₂::OperatorProd; atol::Real=atol, rtol::Real=rtol) -> Bool

Compare two `OperatorProd`s and judge whether they are approximate to each other.
"""
@inline function Base.isapprox(m₁::OperatorProd, m₂::OperatorProd; atol::Real=atol, rtol::Real=rtol)
    isapprox(efficientoperations, contentorder(typeof(m₁), :value)|>Val, dissolve(m₁), dissolve(m₂); atol=atol, rtol=rtol)::Bool
end

"""
    replace(m::OperatorProd; kwargs...) -> typeof(m)

Return a copy of a concrete `OperatorProd` with some of the field values replaced by the keyword arguments.
"""
@inline Base.replace(m::OperatorProd; kwargs...) = replace(efficientoperations, m; kwargs...)

"""
    promote_rule(::Type{M₁}, ::Type{M₂}) where {M₁<:OperatorProd, M₂<:OperatorProd}

Define the promote rule for `OperatorProd` types.
"""
@inline function Base.promote_rule(::Type{M₁}, ::Type{M₂}) where {M₁<:OperatorProd, M₂<:OperatorProd}
    M₁<:M₂ && return M₂
    M₂<:M₁ && return M₁
    r₁, r₂ = rank(M₁), rank(M₂)
    M = r₂==0 ? rawtype(M₁) : r₁==0 ? rawtype(M₂) : typejoin(rawtype(M₁), rawtype(M₂))
    return fulltype(M, promoteparameters(parameterpairs(M₁), parameterpairs(M₂)))
end

"""
    promote_type(::Type{M}, ::Type{V}, ::Val{:*}) where {M<:OperatorProd, V<:Number}
    promote_type(::Type{V}, ::Type{M}, ::Val{:*}) where {M<:OperatorProd, V<:Number}

Define the promote rule for the multiplication between an `OperatorProd` and a scalar.
"""
@inline Base.promote_type(::Type{M}, ::Type{V}, ::Val{:*}) where {M<:OperatorProd, V<:Number} = promote_type(V, M, Val(:*))
@inline function Base.promote_type(::Type{V}, ::Type{M}, ::Val{:*}) where {M<:OperatorProd, V<:Number}
    return reparameter(M, :value, promote_type(valtype(M), V))
end

"""
    getindex(m::OperatorProd, i) -> OperatorProd

Overloaded `[]`.
"""
@inline Base.getindex(m::OperatorProd, i) = rawtype(typeof(m))(dissolve(m, getindex, (i,))...)
@inline dissolve(m::OperatorProd, ::Val{:value}, ::typeof(getindex), ::Tuple{Any}, ::NamedTuple) = one(valtype(m))
@inline dissolve(m::OperatorProd, ::Val{:id}, ::typeof(getindex), i::Tuple{Any}, ::NamedTuple) = ID(id(m)[first(i)])

"""
    length(m::OperatorProd) -> Int

Get the length of an `OperatorProd`.
"""
@inline Base.length(m::OperatorProd) = rank(m)
@inline Base.firstindex(m::OperatorProd) = 1
@inline Base.lastindex(m::OperatorProd) = rank(m)

"""
    one(::Type{M}) where M<:OperatorProd
    one(m::OperatorProd)

Get the identity quantum operator.
"""
@inline function Base.one(::Type{M}) where M<:OperatorProd
    rtype = rawtype(M)
    vtype = isconcretetype(valtype(M)) ? valtype(M) : Int
    @assert fieldnames(rtype) == (:value, :id) "one error: not supported type($(nameof(rtype)))."
    return rtype(one(vtype), ID())
end
@inline Base.one(m::OperatorProd) = rawtype(typeof(m))(dissolve(m, one)...)
@inline dissolve(m::OperatorProd, ::Val{:value}, ::typeof(one), ::Tuple, ::NamedTuple) = one(valtype(m))
@inline dissolve(::OperatorProd, ::Val{:id}, ::typeof(one), ::Tuple, ::NamedTuple) = ID()

"""
    convert(::Type{M}, m::Scalar) where {M<:Scalar}
    convert(::Type{M}, m::Number) where {M<:Scalar}
    convert(::Type{M}, m::OperatorProd) where {M<:OperatorProd}

1) Convert a scalar quantum operator from one type to another;
2) Convert a scalar to a scalar quantum operator;
3) Convert a quantum operator from one type to another.
"""
@inline Base.convert(::Type{M}, m::Scalar) where {M<:Scalar} = typeof(m)<:M ? m : one(M)*value(m)
@inline Base.convert(::Type{M}, m::Number) where {M<:Scalar} = one(M)*m
@inline function Base.convert(::Type{M}, m::OperatorProd) where {M<:OperatorProd}
    (typeof(m) <: M) && return m
    @assert convertible(M, m) "convert error: $(typeof(m)) cannot be converted to $M."
    return replace(m, convert(valtype(M), value(m)))
end
function convertible(::Type{M}, m::OperatorProd) where {M<:OperatorProd}
    !(rawtype(typeof(m)) <: rawtype(M)) && return false
    S = supertype(typeof(m), nameof(M))
    for (i, name) in enumerate(parameternames(M))
        (name != :value) && !(parametertype(S, i) <: parametertype(M, i)) && return false
    end
    return true
end

"""
    split(m::OperatorProd) -> Tuple{Any, Vararg{OperatorProd}}

Split an `OperatorProd` into the coefficient and a sequence of rank-1 `OperatorProd`s.
"""
@inline @generated function Base.split(m::OperatorProd)
    exprs = [:(value(m))]
    for i = 1:rank(m)
        push!(exprs, :(m[$i]))
    end
    return Expr(:tuple, exprs...)
end

"""
    OperatorSum{I<:ID{SingularID}, M<:OperatorProd} <: QuantumOperator{I, M}

The sum of `OperatorProd`s.

Similar items are automatically merged with the aid of the id system.
"""
struct OperatorSum{I<:ID{SingularID}, M<:OperatorProd} <: QuantumOperator{I, M}
    contents::Dict{I, M}
    OperatorSum(contents::Dict{<:ID{SingularID}, <:OperatorProd}) = new{keytype(contents), valtype(contents)}(contents)
end
@inline function Base.promote_rule(::Type{MS₁}, ::Type{MS₂}) where {MS₁<:OperatorSum, MS₂<:OperatorSum}
    I = promote_type(idtype(MS₁), idtype(MS₂))
    M = promote_type(valtype(MS₁), valtype(MS₂))
    return OperatorSum{I, M}
end
function Base.show(io::IO, ms::OperatorSum)
    @printf io "OperatorSum with %s %s:\n" length(ms) nameof(valtype(ms))
    for m in values(ms)
        @printf io "  %s\n" m
    end
end
@inline Base.zero(ms::OperatorSum) = zero(typeof(ms))
@inline Base.zero(::Type{MS}) where {MS<:OperatorSum} = MS()

"""
    OperatorSum(ms::OperatorProd...)
    OperatorSum{I, M}(ms::OperatorProd...) where {I<:ID{SingularID}, M<:OperatorProd}

Get the sum of `OperatorProd`s with similar items merged.
"""
@inline OperatorSum(ms::OperatorSum) = OperatorSum{idtype(ms), valtype(ms)}(ms)
@inline OperatorSum(ms::OperatorProd...) = OperatorSum{idtype(eltype(ms)), eltype(ms)}(ms...)
@inline  OperatorSum{I, M}(ms::OperatorSum) where {I<:ID{SingularID}, M<:OperatorProd} = OperatorSum(Dict{I, M}(ms.contents))
function OperatorSum{I, M}(ms::OperatorProd...) where {I<:ID{SingularID}, M<:OperatorProd}
    contents = Dict{I, M}()
    for m in ms
        abs(value(m))==0 || (contents[id(m)] = m)
    end
    return OperatorSum(contents)
end

"""
    repr(ms::OperatorSum) -> String

Get the repr representation of a sum of `OperatorProd`s.
"""
function Base.repr(ms::OperatorSum)
    cache = [@sprintf("OperatorSum with %s %s:", length(ms), nameof(valtype(ms)))]
    for m in ms|>values
        push!(cache, @sprintf("  %s", repr(m)))
    end
    return join(cache, "\n")
end

"""
    isapprox(ms₁::OperatorSum, ms₂::OperatorSum; atol::Real=atol, rtol::Real=rtol) -> Bool

Compare two `OperatorSum`s and judge whether they are approximate to each other.
"""
@inline function Base.isapprox(ms₁::OperatorSum, ms₂::OperatorSum; atol::Real=atol, rtol::Real=rtol)
    for (k, m) in ms₁
        isapprox(value(m), 0, atol=atol, rtol=rtol) && continue
        haskey(ms₂, k) || return false
        isapprox(m, ms₂[k]) || return false
    end
    for (k, m) in ms₂
        isapprox(value(m), 0,  atol=atol, rtol=rtol) && continue
        haskey(ms₁, k) || return false
        isapprox(m, ms₁[k]) || return false
    end
    return true
end

"""
    add!(ms::OperatorSum) -> typeof(ms)
    add!(ms::OperatorSum, m::Number) -> typeof(ms)
    add!(ms::OperatorSum, m::OperatorProd) -> typeof(ms)
    add!(ms::OperatorSum, mms::OperatorSum) -> typeof(ms)

Get the in-place addition of quantum operators.
"""
@inline add!(ms::OperatorSum) = ms
@inline add!(ms::OperatorSum, m::Number) = add!(ms, one(valtype(ms))*m)
@inline function add!(ms::OperatorSum, m::OperatorProd)
    m = convert(ms|>valtype, m)
    old = get(ms, id(m), nothing)
    new = isnothing(old) ? m : replace(old, value(old)+value(m))
    abs(value(new))==0.0 ? delete!(ms, id(m)) : (ms[id(m)] = new)
    return ms
end
@inline function add!(ms::OperatorSum, mms::OperatorSum)
    for m in mms|>values
        add!(ms, m)
    end
    return ms
end

"""
    sub!(ms::OperatorSum) -> typeof(ms)
    sub!(ms::OperatorSum, m::Number) -> typeof(ms)
    sub!(ms::OperatorSum, m::OperatorProd) -> typeof(ms)
    sub!(ms::OperatorSum, mms::OperatorSum) -> typeof(ms)

Get the in-place subtraction of quantum operators.
"""
@inline sub!(ms::OperatorSum) = ms
@inline sub!(ms::OperatorSum, m::Number) = add!(ms, one(valtype(ms))*(-m))
@inline function sub!(ms::OperatorSum, m::OperatorProd)
    m = convert(ms|>valtype, m)
    old = get(ms, id(m), nothing)
    new = isnothing(old) ? -m : replace(old, value(old)-value(m))
    abs(value(new))==0.0 ? delete!(ms, id(m)) : (ms[id(m)] = new)
    return ms
end
@inline function sub!(ms::OperatorSum, mms::OperatorSum)
    for m in mms|>values
        sub!(ms, m)
    end
    return ms
end

"""
    mul!(ms::OperatorSum, factor::Scalar) -> OperatorSum
    mul!(ms::OperatorSum, factor::Number) -> OperatorSum

Get the in-place multiplication of an `OperatorSum` with a scalar.
"""
@inline mul!(ms::OperatorSum, factor::Scalar) = mul!(ms, value(factor))
function mul!(ms::OperatorSum, factor::Number)
    @assert isa(one(ms|>valtype|>valtype)*factor, ms|>valtype|>valtype) "mul! error: mismatched type, $(ms|>valtype) and $(factor|>typeof)."
    for m in values(ms)
        ms[id(m)] = replace(m, value(m)*factor)
    end
    return ms
end

"""
    div!(ms::OperatorSum, factor::Scalar) -> OperatorSum
    div!(ms::OperatorSum, factor::Number) -> OperatorSum

Get the in-place division of an `OperatorSum` with a scalar.
"""
@inline div!(ms::OperatorSum, factor::Scalar) = div!(ms, value(factor))
function div!(ms::OperatorSum, factor::Number)
    @assert isa(one(ms|>valtype|>valtype)/factor, ms|>valtype|>valtype) "div! error: mismatched type, $(ms|>valtype) and $(factor|>typeof)."
    for m in values(ms)
        ms[id(m)] = replace(m, value(m)/factor)
    end
    return ms
end

"""
    +(m::OperatorProd) -> typeof(m)
    +(m::OperatorProd, factor::Number) -> OperatorSum
    +(factor::Number, m::OperatorProd) -> OperatorSum
    +(m₁::OperatorProd, m₂::OperatorProd) -> OperatorSum
    +(ms::OperatorSum) -> typeof(ms)
    +(ms::OperatorSum, factor::Number) -> OperatorSum
    +(factor::Number, ms::OperatorSum) -> OperatorSum
    +(ms::OperatorSum, m::OperatorProd) -> OperatorSum
    +(m::OperatorProd, ms::OperatorSum) -> OperatorSum
    +(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `+` between quantum operators.
"""
@inline Base.:+(m::OperatorProd) = m
@inline Base.:+(ms::OperatorSum) = ms
@inline Base.:+(factor::Number, m::OperatorProd) = m + one(m)*factor
@inline Base.:+(m::OperatorProd, factor::Number) = m + one(m)*factor
@inline Base.:+(factor::Number, m::Scalar) = replace(m, value(m)+factor)
@inline Base.:+(m::Scalar, factor::Number) = replace(m, value(m)+factor)
@inline Base.:+(m₁::Scalar, m₂::Scalar) = replace(m₁, value(m₁)+value(m₂))
@inline Base.:+(factor::Number, ms::OperatorSum) = ms + one(valtype(ms))*factor
@inline Base.:+(ms::OperatorSum, factor::Number) = ms + one(valtype(ms))*factor
@inline function Base.:+(m₁::OperatorProd, m₂::OperatorProd)
    I = promote_type(m₁|>idtype, m₂|>idtype)
    M = promote_type(typeof(m₁), typeof(m₂))
    return OperatorSum{I, M}(m₁, m₂)
end
@inline Base.:+(m::OperatorProd, ms::OperatorSum) = ms + m
@inline function Base.:+(ms::OperatorSum, m::OperatorProd)
    I = promote_type(ms|>idtype, m|>idtype)
    M = promote_type(valtype(ms), typeof(m))
    return add!(OperatorSum{I, M}(ms), m)
end
@inline function Base.:+(ms₁::OperatorSum, ms₂::OperatorSum)
    I = promote_type(ms₁|>idtype, ms₂|>idtype)
    M = promote_type(valtype(ms₁), valtype(ms₂))
    return add!(OperatorSum{I, M}(ms₁), ms₂)
end

"""
    *(factor::Number, m::OperatorProd) -> OperatorProd
    *(m::OperatorProd, factor::Number) -> OperatorProd
    *(m₁::OperatorProd, m₂::OperatorProd) -> OperatorProd
    *(factor::Number, ms::OperatorSum) -> OperatorSum
    *(ms::OperatorSum, factor::Number) -> OperatorSum
    *(m::OperatorProd, ms::OperatorSum) -> OperatorSum
    *(ms::OperatorSum, m::OperatorProd) -> OperatorSum
    *(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `*` between quantum operators or a quantum operator and a scalar.
"""
@inline Base.:*(factor::Scalar, m::OperatorProd) = m * value(factor)
@inline Base.:*(m::OperatorProd, factor::Scalar) = m * value(factor)
@inline Base.:*(m₁::Scalar, m₂::Scalar) = replace(m₁, value(m₁)*value(m₂))
@inline Base.:*(factor::Number, m::OperatorProd) = m * factor
@inline Base.:*(m::OperatorProd, factor::Number) = replace(m, factor*value(m))
@inline Base.:*(factor::Scalar, ms::OperatorSum) = ms * value(factor)
@inline Base.:*(ms::OperatorSum, factor::Scalar) = ms * value(factor)
@inline Base.:*(factor::Number, ms::OperatorSum) = ms * factor
@inline function Base.:*(ms::OperatorSum, factor::Number)
    abs(factor)==0 && return empty(ms)
    result = OperatorSum{idtype(ms), promote_type(valtype(ms), typeof(factor), Val(:*))}()
    for (id, m) in ms
        result[id] = m * factor
    end
    return result
end
@inline Base.:*(m::OperatorProd, ms::OperatorSum) = OperatorSum((m * mm for mm in ms|>values)...)
@inline Base.:*(ms::OperatorSum, m::OperatorProd) = OperatorSum((mm * m for mm in ms|>values)...)
@inline Base.:*(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum((m₁ * m₂ for m₁ in ms₁|>values for m₂ in ms₂|>values)...)
@inline function Base.:*(m₁::OperatorProd, m₂::OperatorProd)
    @assert((m₁|>typeof|>nameof == m₂|>typeof|>nameof) && (m₁|>typeof|>fieldcount == m₂|>typeof|>fieldcount == 2),
            "\"*\" error: not implemented between $(m₁|>typeof|>nameof) and $(m₂|>typeof|>nameof)."
            )
    return rawtype(typeof(m₁))(value(m₁)*value(m₂), id(m₁)*id(m₂))
end

"""
    -(m::OperatorProd) -> typeof(m)
    -(m::OperatorProd, factor::Number) -> OperatorSum
    -(factor::Number, m::OperatorProd) -> OperatorSum
    -(m₁::OperatorProd, m₂::OperatorProd) -> OperatorSum
    -(ms::OperatorSum) -> typeof(ms)
    -(ms::OperatorSum, factor::Number) -> OperatorSum
    -(factor::Number, ms::OperatorSum) -> OperatorSum
    -(m::OperatorProd, ms::OperatorSum) -> OperatorSum
    -(ms::OperatorSum, m::OperatorProd) -> OperatorSum
    -(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `-` between quantum operators.
"""
@inline Base.:-(m::OperatorProd) = m * (-1)
@inline Base.:-(ms::OperatorSum) = ms * (-1)
@inline Base.:-(factor::Number, m::OperatorProd) = one(m)*factor - m
@inline Base.:-(m::OperatorProd, factor::Number) = m + one(m)*(-factor)
@inline Base.:-(factor::Number, m::Scalar) = replace(m, factor-value(m))
@inline Base.:-(m::Scalar, factor::Number) = replace(m, value(m)-factor)
@inline Base.:-(m₁::Scalar, m₂::Scalar) = replace(m₁, value(m₁)-value(m₂))
@inline Base.:-(factor::Number, ms::OperatorSum) = one(valtype(ms))*factor - ms
@inline Base.:-(ms::OperatorSum, factor::Number) = ms + one(valtype(ms))*(-factor)
@inline function Base.:-(m₁::OperatorProd, m₂::OperatorProd)
    I = promote_type(m₁|>idtype, m₂|>idtype)
    M = promote_type(typeof(m₁), typeof(m₂))
    return OperatorSum{I, M}(m₁, -m₂)
end
@inline function Base.:-(m::OperatorProd, ms::OperatorSum)
    I = promote_type(m|>idtype, ms|>idtype)
    M = promote_type(typeof(m), valtype(ms))
    return sub!(OperatorSum{I, M}(m), ms)
end
@inline function Base.:-(ms::OperatorSum, m::OperatorProd)
    I = promote_type(ms|>idtype, m|>idtype)
    M = promote_type(valtype(ms), typeof(m))
    return sub!(OperatorSum{I, M}(ms), m)
end
@inline function Base.:-(ms₁::OperatorSum, ms₂::OperatorSum)
    I = promote_type(ms₁|>idtype, ms₂|>idtype)
    M = promote_type(valtype(ms₁), valtype(ms₂))
    return sub!(OperatorSum{I, M}(ms₁), ms₂)
end

"""
    /(m::OperatorProd, factor::Number) -> OperatorProd
    /(m::OperatorProd, factor::Scalar) -> OperatorProd
    /(ms::OperatorSum, factor::Number) -> OperatorSum
    /(ms::OperatorSum, factor::Scalar) -> OperatorSum

Overloaded `/` between a quantum operator and a scalar.
"""
Base.:/(m::OperatorProd, factor::Number) = m * (one(valtype(m))/factor)
Base.:/(m::OperatorProd, factor::Scalar) = m * (one(valtype(m))/value(factor))
Base.:/(ms::OperatorSum, factor::Number) = ms * (one(valtype(valtype(ms)))/factor)
Base.:/(ms::OperatorSum, factor::Scalar) = ms * (one(valtype(valtype(ms)))/value(factor))

"""
    //(m::OperatorProd, factor::Integer) -> OperatorProd
    //(m::OperatorProd, factor::Scalar) -> OperatorProd
    //(ms::OperatorSum, factor::Integer) ->  OperatorSum
    //(ms::OperatorSum, factor::Scalar) -> OperatorSum

Overloaded `//` between a quantum operator and a scalar.
"""
Base.://(m::OperatorProd, factor::Integer) = m * (1//factor)
Base.://(m::OperatorProd, factor::Scalar) = m * (1//value(factor))
Base.://(ms::OperatorSum, factor::Number) = ms * (1//factor)
Base.://(ms::OperatorSum, factor::Scalar) = ms * (1//value(factor))

"""
    ^(m::OperatorProd, n::Integer) -> OperatorProd
    ^(ms::OperatorSum, n::Integer) -> OperatorSum

Overloaded `^` between a quantum operator and an integer.
"""
Base.:^(m::OperatorProd, n::Integer) = (@assert n>0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->m, Val(n))))
Base.:^(ms::OperatorSum, n::Integer) = (@assert n>0 "^ error: non-positive integers are not allowed."; prod(ntuple(i->ms, Val(n))))

"""
    ⊗(m::OperatorProd, ms::OperatorSum) -> OperatorSum
    ⊗(ms::OperatorSum, m::OperatorProd) -> OperatorSum
    ⊗(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `⊗` between quantum operators.
"""
⊗(m::OperatorProd, ms::OperatorSum) = OperatorSum((m ⊗ mm for mm in ms|>values)...)
⊗(ms::OperatorSum, m::OperatorProd) = OperatorSum((mm ⊗ m for mm in ms|>values)...)
⊗(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum((m₁ ⊗ m₂ for m₁ in ms₁|>values for m₂ in ms₂|>values)...)

"""
    ⋅(m::OperatorProd, ms::OperatorSum) -> OperatorSum
    ⋅(ms::OperatorSum, m::OperatorProd) -> OperatorSum
    ⋅(ms₁::OperatorSum, ms₂::OperatorSum) -> OperatorSum

Overloaded `⋅` between quantum operators.
"""
⋅(m::OperatorProd, ms::OperatorSum) = OperatorSum((m ⋅ mm for mm in ms|>values)...)
⋅(ms::OperatorSum, m::OperatorProd) = OperatorSum((mm ⋅ m for mm in ms|>values)...)
⋅(ms₁::OperatorSum, ms₂::OperatorSum) = OperatorSum((m₁ ⋅ m₂ for m₁ in ms₁|>values for m₂ in ms₂|>values)...)

"""
    sequence(m::OperatorProd, table) -> NTuple{rank(m), Int}

Get the sequence of the id of a quantum operator according to a table.
"""
@generated sequence(m::OperatorProd, table) = Expr(:tuple, [:(get(table, id(m)[$i], nothing)) for i = 1:rank(m)]...)

"""
    replace(m::OperatorProd, pairs::Pair{<:SingularID, <:Union{OperatorProd, OperatorSum}}...) -> OperatorProd/OperatorSum
    replace(ms::OperatorSum, pairs::Pair{<:SingularID, <:Union{OperatorProd, OperatorSum}}...) -> OperatorSum

Replace the rank-1 components of a quantum operator with new quantum operators.
"""
function Base.replace(m::OperatorProd, pairs::Pair{<:SingularID, <:Union{OperatorProd, OperatorSum}}...)
    replacedids = NTuple{length(pairs), idtype(m)|>eltype}(pair.first for pair in pairs)
    ms = split(m)
    result = ms[1]
    for i = 1:rank(m)
        index = findfirst(isequal(id(m)[i]), replacedids)
        result = result * (isa(index, Int) ? pairs[index].second : ms[i+1])
    end
    return result
end
function Base.replace(ms::OperatorSum, pairs::Pair{<:SingularID, <:Union{OperatorProd, OperatorSum}}...)
    result = operatorsumtype(pairs...)()
    for m in values(ms)
        add!(result, replace(m, pairs...))
    end
    return result
end
@generated function operatorsumtype(ms::Vararg{Pair{<:SingularID, <:Union{OperatorProd, OperatorSum}}, N}) where N
    M = Union{}
    for i = 1:N
        E = fieldtype(ms[i], 2)
        M = promote_type(E<:OperatorSum ? valtype(E) : E, M)
    end
    M = reparameter(M, :id, ID{eltype(idtype(M))})
    return OperatorSum{idtype(M), M}
end

"""
    permute(id₁::SingularID, id₂::SingularID) -> Tuple{Vararg{OperatorProd}}

Permutation rule of two singular ids.
"""
permute(::SingularID, ::SingularID) = error("permute error: not implemented for $(nameof(T)).")

"""
    permute!(result::OperatorSum, m::OperatorProd, table) -> OperatorSum
    permute!(result::OperatorSum, ms::OperatorSum, table) -> OperatorSum

Permute the singular ids of a quantum operator to the descending order according to a table, and store the permuted quantum operators in result.

!!! note
    To use this function, the user must implement a method of `permute`, which computes the result of the permutation of two ids and takes the following interface:
    ```julia
    permute(id₁::SingularID, id₂::SingularID) -> Union{OperatorProd, OperatorSum}
    ```
    Here, `id₁` and `id₂` are two arbitrary singular ids contained in `id(m)`.
"""
function Base.permute!(result::OperatorSum, m::OperatorProd, table)
    cache = valtype(result)[m]
    while length(cache) > 0
        current = pop!(cache)
        pos = operatorprodcommuteposition(sequence(current, table))
        if isa(pos, Nothing)
            add!(result, current)
        else
            left = current[1:pos-1] * value(m)
            right = current[pos+2:end]
            for middle in permute(id(current)[pos], id(current)[pos+1])
                temp = left * middle * right
                (temp === nothing) || push!(cache, temp)
            end
        end
    end
    return result
end
function Base.permute!(result::OperatorSum, ms::OperatorSum, table)
    for m in values(ms)
        permute!(result, m, table)
    end
    return result
end
function operatorprodcommuteposition(sequences)
    pos = 1
    while pos < length(sequences)
        (sequences[pos] < sequences[pos+1]) && return pos
        pos += 1
    end
    return nothing
end

"""
    permute(m::OperatorProd, table) -> OperatorSum
    permute(ms::OperatorSum, table) -> OperatorSum

Permute the singular ids of a quantum operator to the descending order according to a table.
"""
@inline function permute(m::OperatorProd, table)
    M = reparameter(typeof(m), :id, ID{m|>idtype|>eltype})
    permute!(OperatorSum{idtype(M), M}(), m, table)
end
@inline function permute(ms::OperatorSum, table)
    M = reparameter(valtype(ms), :id, ID{ms|>idtype|>eltype})
    permute!(OperatorSum{idtype(M), M}(), ms, table)
end

"""
    Transformation <: Function

Abstract transformation that could transform the quantum operators.
"""
abstract type Transformation <: Function end
@inline Base.:(==)(transformation₁::Transformation, transformation₂::Transformation) = ==(efficientoperations, transformation₁, transformation₂)
@inline Base.isequal(transformation₁::Transformation, transformation₂::Transformation) = isequal(efficientoperations, transformation₁, transformation₂)
@inline Base.valtype(transformation::Transformation, m::OperatorProd) = valtype(typeof(transformation), typeof(m))
@inline Base.valtype(transformation::Transformation, ms::OperatorSum) = valtype(typeof(transformation), typeof(ms))
@inline Base.valtype(T::Type{<:Transformation}, MS::Type{<:OperatorSum}) = OperatorSum{idtype(valtype(T, valtype(MS))), valtype(T, valtype(MS))}

"""
    (transformation::Transformation)(ms::OperatorSum) -> OperatorSum

Get the transformed quantum operators.
"""
function (transformation::Transformation)(ms::OperatorSum)
    result = valtype(transformation, ms)()
    for m in values(ms)
        add!(result, transformation(m))
    end
    return result
end

"""
    Identity <: Transformation

The identity transformation.
"""
struct Identity <: Transformation end
@inline Base.valtype(::Type{Identity}, M::Type{<:OperatorProd}) = M
@inline (i::Identity)(m::OperatorProd) = m

"""
    Numericalization{T<:Number} <: Transformation

The numericalization transformation, which converts the value of a quantum operator to a number of type `T`.
"""
struct Numericalization{T<:Number} <: Transformation end
@inline Base.valtype(num::Numericalization) = valtype(typeof(num))
@inline Base.valtype(::Type{<:Numericalization{T}}) where {T<:Number} = T
@inline Base.valtype(T::Type{<:Numericalization}, M::Type{<:OperatorProd}) = reparameter(M, :value, valtype(T))
@inline (n::Numericalization)(m::OperatorProd) = convert(valtype(n, m), m)

"""
    MatrixRepresentation <: Transformation

The matrix representation.
"""
abstract type MatrixRepresentation <: Transformation end

"""
    matrix

Generic matrix representation.
"""
function matrix end

"""
    matrix!

Generic matrix representation.
"""
function matrix! end

end #module


text:
function MOIU.map_indices(index_map::Function, change::DecisionCoefficientChange)
    return DecisionCoefficientChange(index_map(change.decision), change.new_coefficient)
end
MOIU.map_indices(index_map::Function, change::DecisionStateChange) = change
MOIU.map_indices(index_map::Function, change::KnownValuesChange) = change

function MOIU.modify_function(f::AffineDecisionFunction, change::Union{MOI.ScalarConstantChange, MOI.ScalarCoefficientChange})
    return typeof(f)(MOIU.modify_function(f.variable_part, change),
                     copy(f.decision_part))
end

function MOIU.modify_function(f::QuadraticDecisionFunction{T}, change::Union{MOI.ScalarConstantChange, MOI.ScalarCoefficientChange}) where T
    return typeof(f)(
        MOIU.modify_function(f.variable_part, change),
        copy(f.decision_part),
        copy(f.cross_terms))
end

function MOIU.modify_function(f::VectorAffineDecisionFunction, change::Union{MOI.VectorConstantChange, MOI.MultirowChange})
    return typeof(f)(MOIU.modify_function(f.variable_part, change),
                     copy(f.decision_part))
end

function MOIU.modify_function(f::AffineDecisionFunction, change::DecisionCoefficientChange)
    return typeof(f)(copy(f.variable_part),
                     MOIU.modify_function(f.decision_part,
                                          MOI.ScalarCoefficientChange(change.decision, change.new_coefficient)))
end

function MOIU.modify_function(f::QuadraticDecisionFunction{T}, change::DecisionCoefficientChange) where T
    return typeof(f)(
        copy(f.variable_part),
        MOIU.modify_function(f.decision_part,
                             MOI.ScalarCoefficientChange(change.decision, change.new_coefficient)),
        copy(f.cross_terms))
end

function MOIU.modify_function(f::VectorAffineDecisionFunction, change::DecisionMultirowChange)
    return typeof(f)(copy(f.variable_part),
                     MOIU.modify_function(f.decision_part,
                                          DecisionMultirowChange(change.decision, change.new_coefficients)))
end

function MOIU.modify_function(f::Union{SingleDecision, AffineDecisionFunction, QuadraticDecisionFunction, VectorAffineDecisionFunction}, change::Union{DecisionStateChange, KnownValuesChange})
    # Nothing to do here, handled in bridges
    return f
end

# Can rely on scalarize bridge if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.ScalarizeBridge{T, AffineDecisionFunction{T}},
                    change::DecisionModification) where T
    for constraint in bridge.scalar_constraints
        MOI.modify(model, constraint, change)
    end
    return nothing
end

# Can rely on vectorize bridge if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.VectorizeBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.vector_constraint, change)
    return nothing
end

# Can rely on flipsign bridges if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.FlipSignBridge,
                    change::DecisionModification)
    MOI.modify(model, bridge.constraint, change)
    return nothing
end

# Can rely on norm-bridges if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.NormInfinityBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.constraint, change)
    return nothing
end
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.NormOneBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.nn_index, change)
    return nothing
end

function MOI.add_constraint(uf::MOIU.UniversalFallback, f::SingleDecision, s::MOI.AbstractScalarSet)
    F = typeof(f)
    S = typeof(s)
    # Add constraint as usual but force index
    constraints = get!(uf.constraints, (F, S)) do
        OrderedDict{CI{F, S}, Tuple{F, S}}()
    end::OrderedDict{CI{F, S}, Tuple{F, S}}
    ci = CI{typeof(f), typeof(s)}(f.decision.value)
    constraints[ci] = (f, s)
    return ci
end

function MOI.modify(uf::MOIU.UniversalFallback, ::CI{F,S}, ::Union{DecisionStateChange,KnownValuesChange}) where {F, S}
    return nothing
end

function MOI.delete(uf::MOIU.UniversalFallback, ci::CI{F, S}) where {F <: SingleDecision, S <: MOI.AbstractScalarSet}
    decision_ci = CI{MOI.SingleVariable, SingleDecisionSet{Float64}}(ci.value)
    if !MOI.is_valid(uf, ci) && MOI.is_valid(uf, decision_ci)
        throw(MOI.InvalidIndex(ci))
    end
    # Delete constraint as usual
    if MOI.is_valid(uf, ci)
        delete!(uf.constraints[(F, S)], ci)
        delete!(uf.con_to_name, ci)
        uf.name_to_con = nothing
        for d in values(uf.conattr)
            delete!(d, ci)
        end
    end
    # Update SingleDecisionSet for printing
    if MOI.is_valid(uf, decision_ci)
        set = MOI.get(uf, MOI.ConstraintSet(), decision_ci)
        new_set = SingleDecisionSet(set.stage, set.decision, NoSpecifiedConstraint(), set.is_recourse)
        MOI.set(uf, MOI.ConstraintSet(), decision_ci, new_set)
    end
    return nothing
end

function MOI.get(model::MOIU.CachingOptimizer,
                 attr::DecisionIndex,
                 ci::CI)
    return MOI.get(model.optimizer, attr, model.model_to_optimizer_map[ci])
end

function MOI.get(b::MOIB.AbstractBridgeOptimizer,
                 attr::DecisionIndex, ci::MOI.ConstraintIndex)
    return MOIB.call_in_context(b, ci, bridge -> MOI.get(b, attr, bridge))
end

ext: jl


text:
# MobileNetv2

"""
    mobilenetv2(width_mult, configs; max_width = 1280, nclasses = 1000)

Create a MobileNetv2 model.
([reference](https://arxiv.org/abs/1801.04381)).

# Arguments
- `width_mult`: Controls the number of output feature maps in each block
                (with 1.0 being the default in the paper)
- `configs`: A "list of tuples" configuration for each layer that details:
  - `t`: The expansion factor that controls the number of feature maps in the bottleneck layer
  - `c`: The number of output feature maps
  - `n`: The number of times a block is repeated
  - `s`: The stride of the convolutional kernel
- `max_width`: The maximum number of feature maps in any layer of the network
- `nclasses`: The number of output classes
"""
function mobilenetv2(width_mult, configs; max_width = 1280, nclasses = 1000)
  # building first layer
  inplanes = _round_channels(32 * width_mult, width_mult == 0.1 ? 4 : 8)
  layers = []
  append!(layers, conv_bn((3, 3), 3, inplanes, stride = 2))

  # building inverted residual blocks
  for (t, c, n, s) in configs
    outplanes = _round_channels(c * width_mult, width_mult == 0.1 ? 4 : 8)
    for i in 1:n
      push!(layers, invertedresidual(3, inplanes, inplanes * t, outplanes, relu6;
                                     stride = i == 1 ? s : 1))
      inplanes = outplanes
    end
  end

  # building last several layers
  outplanes = (width_mult > 1) ? _round_channels(max_width * width_mult, width_mult == 0.1 ? 4 : 8) :
                                 max_width

  return Chain(Chain(layers...,
                     conv_bn((1, 1), inplanes, outplanes, relu6, bias = false)...),
               Chain(AdaptiveMeanPool((1, 1)), MLUtils.flatten, Dense(outplanes, nclasses)))
end

# Layer configurations for MobileNetv2
const mobilenetv2_configs = [
  # t, c, n, s
  (1, 16, 1, 1),
  (6, 24, 2, 2),
  (6, 32, 3, 2),
  (6, 64, 4, 2),
  (6, 96, 3, 1),
  (6, 160, 3, 2),
  (6, 320, 1, 1)
]

# Model definition for MobileNetv2
struct MobileNetv2
  layers
end

"""
    MobileNetv2(width_mult = 1.0; pretrain = false, nclasses = 1000)

Create a MobileNetv2 model with the specified configuration.
([reference](https://arxiv.org/abs/1801.04381)).
Set `pretrain` to `true` to load the pretrained weights for ImageNet.

# Arguments
- `width_mult`: Controls the number of output feature maps in each block
                (with 1.0 being the default in the paper;
                 this is usually a value between 0.1 and 1.4)
- `pretrain`: Whether to load the pre-trained weights for ImageNet
- `nclasses`: The number of output classes

See also [`Metalhead.mobilenetv2`](#).
"""
function MobileNetv2(width_mult::Number = 1; pretrain = false, nclasses = 1000)
  layers = mobilenetv2(width_mult, mobilenetv2_configs; nclasses = nclasses)
  pretrain && loadpretrain!(layers, string("MobileNetv2"))

  MobileNetv2(layers)
end

@functor MobileNetv2

(m::MobileNetv2)(x) = m.layers(x)

backbone(m::MobileNetv2) = m.layers[1]
classifier(m::MobileNetv2) = m.layers[2:end]

# MobileNetv3

"""
    mobilenetv3(width_mult, configs; max_width = 1024, nclasses = 1000)

Create a MobileNetv3 model.
([reference](https://arxiv.org/abs/1905.02244)).

# Arguments
- `width_mult`: Controls the number of output feature maps in each block
                (with 1.0 being the default in the paper;
                 this is usually a value between 0.1 and 1.4)
- `configs`: a "list of tuples" configuration for each layer that details:
  - `k::Int` - The size of the convolutional kernel
  - `c::Float` - The multiplier factor for deciding the number of feature maps in the hidden layer
  - `t::Int` - The number of output feature maps for a given block
  - `r::Int` - The reduction factor (`>= 1` or `nothing` to skip) for squeeze and excite layers
  - `use_hs::Bool` - Whether to use Hard-Swish activation function
  - `s::Int` - The stride of the convolutional kernel
- `max_width`: The maximum number of feature maps in any layer of the network
- `nclasses`: the number of output classes
"""
function mobilenetv3(width_mult, configs; max_width = 1024, nclasses = 1000)
  # building first layer
  inplanes = _round_channels(16 * width_mult, 8)
  layers = []
  append!(layers, conv_bn((3, 3), 3, inplanes, hardswish; stride = 2))
  explanes = 0
  # building inverted residual blocks
  for (k, t, c, r, use_hs, s) in configs
    # inverted residual layers
    outplanes = _round_channels(c * width_mult, 8)
    explanes = _round_channels(inplanes * t, 8)
    activation = use_hs ? hardswish : relu
    push!(layers, invertedresidual(k, inplanes, explanes, outplanes, activation;
                                   stride = s, reduction = r))
    inplanes = outplanes
  end

  # building last several layers
  output_channel = max_width
  output_channel = width_mult > 1.0 ? _round_channels(output_channel * width_mult, 8) : output_channel
  classifier = (
    Dense(explanes, output_channel, hardswish),
    Dropout(0.2),
    Dense(output_channel, nclasses),
  )

  return Chain(Chain(layers...,
                     conv_bn((1, 1), inplanes, explanes, hardswish, bias = false)...),
               Chain(AdaptiveMeanPool((1, 1)), MLUtils.flatten, classifier...))
end

# Configurations for small and large mode for MobileNetv3
mobilenetv3_configs = Dict(
  :small => [
    # k, t, c, SE, HS, s
    (3, 1, 16, 4, false, 2),
    (3, 4.5, 24, nothing, false, 2),
    (3, 3.67, 24, nothing, false, 1),
    (5, 4, 40, 4, true, 2),
    (5, 6, 40, 4, true, 1),
    (5, 6, 40, 4, true, 1),
    (5, 3, 48, 4, true, 1),
    (5, 3, 48, 4, true, 1),
    (5, 6, 96, 4, true, 2),
    (5, 6, 96, 4, true, 1),
    (5, 6, 96, 4, true, 1),
  ],
  :large => [
    # k, t, c, SE, HS, s
    (3, 1, 16, nothing, false, 1),
    (3, 4, 24, nothing, false, 2),
    (3, 3, 24, nothing, false, 1),
    (5, 3, 40, 4, false, 2),
    (5, 3, 40, 4, false, 1),
    (5, 3, 40, 4, false, 1),
    (3, 6, 80, nothing, true, 2),
    (3, 2.5, 80, nothing, true, 1),
    (3, 2.3, 80, nothing, true, 1),
    (3, 2.3, 80, nothing, true, 1),
    (3, 6, 112, 4, true, 1),
    (3, 6, 112, 4, true, 1),
    (5, 6, 160, 4, true, 2),
    (5, 6, 160, 4, true, 1),
    (5, 6, 160, 4, true, 1)
  ]
)

# Model definition for MobileNetv3
struct MobileNetv3
  layers
end

"""
    MobileNetv3(mode::Symbol = :small, width_mult::Number = 1; pretrain = false, nclasses = 1000)

Create a MobileNetv3 model with the specified configuration.
([reference](https://arxiv.org/abs/1905.02244)).
Set `pretrain = true` to load the model with pre-trained weights for ImageNet.

# Arguments
- `mode`: :small or :large for the size of the model (see paper).
- `width_mult`: Controls the number of output feature maps in each block
                (with 1.0 being the default in the paper;
                 this is usually a value between 0.1 and 1.4)
- `pretrain`: whether to load the pre-trained weights for ImageNet
- `nclasses`: the number of output classes

See also [`Metalhead.mobilenetv3`](#).
"""
function MobileNetv3(mode::Symbol = :small, width_mult::Number = 1; pretrain = false, nclasses = 1000)
  @assert mode in [:large, :small] "`mode` has to be either :large or :small"

  max_width = (mode == :large) ? 1280 : 1024
  layers = mobilenetv3(width_mult, mobilenetv3_configs[mode]; max_width = max_width, nclasses = nclasses)
  pretrain && loadpretrain!(layers, string("MobileNetv3", mode))
  MobileNetv3(layers)
end

@functor MobileNetv3

(m::MobileNetv3)(x) = m.layers(x)

backbone(m::MobileNetv3) = m.layers[1]
classifier(m::MobileNetv3) = m.layers[2:end]


text:
base_dir = dirname(dirname(pathof(PowerSystems)))

sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))

@testset "Check bus index" begin
    @test sort([b.number for b in collect(get_components(Bus, sys))]) == [1, 2, 3, 4, 10]
    @test sort(collect(Set([
        b.arc.from.number for b in collect(get_components(Branch, sys))
    ]))) == [1, 2, 3, 4]
    @test sort(collect(Set([
        b.arc.to.number for b in collect(get_components(Branch, sys))
    ]))) == [2, 3, 4, 10]

    # TODO: add test for loadzones testing MAPPING_BUSNUMBER2INDEX

end

@testset "Test unique bus numbers" begin
    number = 100
    bus1 = Bus(;
        number = number,
        name = "bus100",
        bustype = BusTypes.PV,
        angle = 1.0,
        voltage = 1.0,
        voltagelimits = (min = -1.0, max = 1.0),
        basevoltage = 1.0,
        area = nothing,
        load_zone = nothing,
    )
    bus2 = Bus(;
        number = number,
        name = "bus101",
        bustype = BusTypes.PV,
        angle = 1.0,
        voltage = 1.0,
        voltagelimits = (min = -1.0, max = 1.0),
        basevoltage = 1.0,
        area = nothing,
        load_zone = nothing,
    )

    add_component!(sys, bus1)
    @test_throws ArgumentError add_component!(sys, bus2)
end


text:
""" Main file """

# Implementation of a bi-objective branch and bound for monodimensional knapsack problem.
# The GeneticAlgorithms module provides two EMOA for enriching the initial lower bound.
# The branchAndBound! function in main.jl performs the enumeration to compute all non dominated solutions of the pareto front.

#------------------------------ STEPS ---------------------------------#

# The algorithm works as follows:

# First step, run the dichotomy methoc in order to compute the X_SE1 set, the set of all non dominated and supported solutions.
# Second step, run the primal heuristic - in this case one of the genetic algorithms available -
#   using the lower bound computed at step 1 as seeding solutions.
# Third step, start the branch and bound with the enriched lower bound obtained at step 2.

# All these steps are visible in the function main(...)

#------------------------------- COMPONENTS ---------------------------#

# Several configurations exist that decide which tools will be used during the execution.
# In order to change these tools, modify the default constructor of the object `Components`.

# The first dichotomy method solves a sequence of scalarized problems. The solver used to solve these scalarizations,
#   can be GLPK (called by JuMP) or COMBO. COMBO is an efficient knapsack solver. The binary for COMBO is located in deps/

# The method for computing the upperbound also solves a sequence of scalarized problems. The solver used can be GLPK or COMBO
#   and thus provides a feasible discrete solution. But it can also be a linear relaxation providing an unfeasible continuous solution.

# The nadirShift field of Components is a boolean activating or not an optimization on the fathoming of subproblems.
#   We made this an option in order to measure the improvements of this feature. Leave it on.

# The primal_heuristic field of Components is a boolean activating or not the use of the genetic algorithm.

# All functions are public

#--------------------------- LOAD LIBRARIES ---------------------------#

println("\n\$ Loading libraries")

using Printf
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "Printf", "")
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "DataStructures", "(LinkedList)")
using DataStructures
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "Random", "")
using Random
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "JuMP", "(FrameWork for GLKP solver)")
using JuMP
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "GLPK", "(Solver)")
using GLPK
@printf("%-6s %-7s %-21s %s\n", "[INFO]", "using", "GeneticAlgorithms", "(Primal Heuristic)")
include("GeneticAlgorithms/main.jl")

#--------------------------- GLOBAL VARIABLES ---------------------------#

const comboPath = joinpath(@__DIR__,"..","deps","libcombo.so")

# JUMP (calling GLPK) or COMBO (the efficient KP solver)
@enum FirstDicho JUMP COMBO
# EXACT_JUMP (calling GLPK), EXACT_COMBO calling deps/libcombo.so, RELAX_LIN_CLASSIC solving linear relaxation
@enum MethodUB EXACT_JUMP EXACT_COMBO RELAX_LIN_CLASSIC

# Algorithm components selected for the execution
mutable struct Components
    firstDicho::FirstDicho
    methodUB::MethodUB
    nadirsShift::Bool
    primal_heuristic::Bool

    Components() = new(COMBO, EXACT_COMBO, true, true) # best config possible (see experimentations)
end

# Configuration
mutable struct Config
    debug::Bool
    # verbose
    # ...

    Config() = new(false) # by default we are in Release mode
end

CONFIG     = Config()
COMPONENTS = Components()

#--------------------------- CODE ---------------------------#

println("\n\$ Loading source files")

# Elementary
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "structs", "(core)")
include("structs.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "GAInterface", "(GeneticAlgorithms)")
include("GAInterface.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "utils", "(core)")
include("utils.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "assignment", "(core)")
include("assignment.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "debug", "(core)")
include("debug.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "files", "(core)")
include("files.jl")

# Algorithms
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "solve1OKP", "(core)")
include("solve1OKP.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "common", "(core)")
include("common.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "lowerBound", "(core)")
include("lowerBound.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "upperBound", "(core)")
include("upperBound.jl")
@printf("%-6s %-7s %-21s %-s\n", "[INFO]", "include", "subProblemStudy", "(core)")
include("subProblemStudy.jl")

println("")

"""
Returns the new nadir points to study and the time spent on computation. Modifies the LB in place.
"""
function branchAndBound!(prob::BiOKP, LB::LowerBound, assignment::Assignment, nadirsToStudy::Vector{PairOfSolution}; verbose = false, parentToChild = ParentToChild(), num::Int = 1, iterator = nothing, timeMax = nothing, start = nothing)

    CONFIG.debug && DEBUG_LB(prob, LB)
    CONFIG.debug && DEBUG_nadirs(prob, nadirsToStudy)
    CONFIG.debug && DEBUG_parentToChild(prob, parentToChild)

    # extreme scenario : exceed time limit
	if timeMax != nothing && time() - start > timeMax
		return nadirsToStudy, false
	end

	testTime = false

	iterator != nothing && (iterator.value += 1)

    # Study the subproblem
    prunedType, dominatedNadirs, nonDomNadirs, UB = SPS_partially_implicit(prob, assignment, nadirsToStudy, parentToChild = parentToChild)

    # update the LB if necessary
    if prunedType == NONE || prunedType == OPTIMALITY
		newNadirsToStudy = LB_update!(prob, LB, dominatedNadirs, UB)
    end

	if prunedType == NONE && assignment.lastAssigned < prob.nbVar
		canAddVar = (assignment.weight + prob.weights[assignment.lastAssigned + 1] <= prob.maxWeight)

        # if the variable can be added without breaking the capacity constraint
		if canAddVar
			ASSIGNMENT_addVar!(assignment, prob)
			parentToChild = UB_computeParentToChild(prob, UB, assignment.lastAssigned, 1)

            newNadirsToStudy, testTime = branchAndBound!(prob,
                                                LB,
												assignment,
												newNadirsToStudy,
												parentToChild = parentToChild,
												num = num + 1,
												iterator = iterator,
												timeMax = timeMax,
												start = start)
		end

		ASSIGNMENT_removeVar!(assignment, prob, canAddVar)
		parentToChild = UB_computeParentToChild(prob, UB, assignment.lastAssigned, 0)

        newNadirsToStudy, testTime = branchAndBound!(prob,
											LB,
											assignment,
											newNadirsToStudy,
											parentToChild = parentToChild,
											num = num + 1,
											iterator = iterator,
											timeMax = timeMax,
											start = start)

		ASSIGNMENT_backtrack!(assignment)

        CONFIG.debug && DEBUG_LB(prob, LB)
        CONFIG.debug && DEBUG_nadirs(prob, newNadirsToStudy)
        CONFIG.debug && DEBUG_nadirs(prob, nonDomNadirs)
        CONFIG.debug && DEBUG_parentToChild(prob, parentToChild)

		return mergeNadirs(prob, newNadirsToStudy, nonDomNadirs), (timeMax == nothing || time() - start < timeMax)
	elseif prunedType == OPTIMALITY

        CONFIG.debug && DEBUG_LB(prob, LB)
        CONFIG.debug && DEBUG_nadirs(prob, newNadirsToStudy)
        CONFIG.debug && DEBUG_nadirs(prob, nonDomNadirs)
        CONFIG.debug && DEBUG_parentToChild(prob, parentToChild)

		return mergeNadirs(prob,newNadirsToStudy, nonDomNadirs), (timeMax == nothing || time() - start < timeMax)
	else

        CONFIG.debug && DEBUG_LB(prob, LB)
        CONFIG.debug && DEBUG_nadirs(prob, nadirsToStudy)
        CONFIG.debug && DEBUG_parentToChild(prob, parentToChild)

		return nadirsToStudy, (timeMax == nothing || time() - start < timeMax)
	end
end

"""
Returns the `LB` for the given problem. Also specifies the number of subproblems studied during execution.

The problem `prob` is assumed to have integer constants and exactly two objectives.

`WARNING` : solving medium to large instances can lead to excessively long resolutions.

# EXAMPLES

```jldoctest
julia> p = FILES_readProblem("Inst_10_1")
BiOKP(10, 2, [50.0 50.0 … 22.0 7.0; 42.0 20.0 … 1.0 16.0], [13.0, 41.0, 29.0, 27.0, 21.0, 30.0, 43.0, 42.0, 14.0, 40.0], 120.0, true)

julia> main(p)
(31,
[93.0, 157.0]
[136.0, 142.0]
[148.0, 141.0]
[153.0, 137.0]
[164.0, 123.0]
[177.0, 118.0]
[189.0, 117.0]
[191.0, 87.0]
)
```

# DIDACTIC

```jldoctest
julia> main()
(19,
[10.0, 15.0]
[13.0, 14.0]
[21.0, 12.0]
[22.0, 10.0]
[30.0, 7.0]
)
```
"""
function main(prob::BiOKP; timeMax = nothing, start = nothing)

	@assert prob.nbObj == 2 "This Branch and Bound supports only a Bio-objective problem"

    # secure resolution parameters
    if !prob.isInteger
        @assert false "The current algorithm only support instances with integer parameters"
        if COMPONENTS.methodUB == EXACT_COMBO
            COMPONENTS.methodUB = EXACT_JUMP
        end
        COMPONENTS.nadirsShift = false
    end

	assignment = Assignment(prob)

	compt = Iterator()

    # compute the set of non dominated and supported solutions (also called lower bound)
    LB = LB_dichoSearch(prob, assignment)

    # if primal heuristic is enabled, then the EMOA is called to add new solutions to the lower bound
    if COMPONENTS.primal_heuristic && prob.nbVar >= 50 # don't call the GAs on very small instances
        momkp        = GAINTERFACE_exportProb(prob)
        seeding_sols = GAINTERFACE_exportLB(prob, LB)
        # by default time spent on GAs is equal to prob.nbVar/10
        front = GeneticAlgorithms.solve(momkp, max_t = Float64(prob.nbVar)/10, seeding_sols = seeding_sols)
        LB           = GAINTERFACE_importLB(prob, front)
    end

	nadirs = LB_getNadirs(prob, LB)

    # launch the enumeration
	branchAndBound!(prob, LB, assignment, nadirs, parentToChild = ParentToChild(), iterator = compt, timeMax = timeMax, start = start)

    CONFIG.debug && DEBUG_LB(prob,LB)

	return compt.value, LB
end

function main(fname::String)
	prob = BiOKP(fname)
	return main(prob)
end

function main()
	prob = BiOKP()
	return main(prob)
end


text:
#=
# Variable struct and related functions
=#

mutable struct Variable
    symbol::Symbol          # symbol used in expressions referring to this variable
    symvar::Array{Basic,1}  # Array of symbolic layer variable symbols(Basic symbols)
    index::Int              # index in the Femshop list of variables
    type::String            # constants for SCALAR, VECTOR, etc.
    location::String        # constant for NODAL or CELL
    values::Array{Float64}  # a C x N array, C is number of components (SCALAR=1, etc.)

    # These are mainly for ARRAY type vars. Should there be a separate struct for those?
    indexer                  # indexer object or array of them
    total_components::Int   # Number of components
    #transposed::Bool        # If true, the values are actually trasposed -> N x C array (better for large C?)

    # not currently used, may be removed
    dependson               # a list of variables that this depends on
    ready::Bool             # Is this variable's values ready? Can dependent variables use it?
end


text:
# test Metropolis
@testset "metropolis" begin
    for tc in testcases
        @show tc.dist
        pdf = make_lpdf(tc.dist)
        @inferred pdf(tc.theta0)
        thetas, accept_ratio, logdensities, blobs =
            metropolis(pdf, tc.sample_ppdf, tc.theta0;
                       niter=tc.niter,
                       hasblob=tc.hasblob,
                       use_progress_meter=false,
                       tc.otherkws...);
        !tc.hasblob && @test blobs==nothing
        @test length(thetas)==tc.niter÷2
        @test length(logdensities)==tc.niter÷2
        @test 0.15<accept_ratio<0.45
        test_mean_std(thetas, tc, tc.tolm)
        test_blobs(tc.blob_truths, blobs)
        #test_logdensities(tc.truths, logdensities, pdf)
    end
end


text:
# TEMPLATE STEP 1
# 1. Create an environment for CpelNano in your home directory (only once)
# 2. Install CpelNano in the environment using (only once)
#   add https://github.com/jordiabante/CpelNano.jl.git#NatureGeneticsModel
# 3. Copy the content of this file in a file in your home directory
# 4. Replace the environment path for the correct one where the arrow points
# 5. Define the proper paths

# Load dependencies
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate("$HOME/julia_env/CpelNano/") # <= create environment
@everywhere using CpelNano

# IO
outdir = ""
dataDir = ""
bam = "$(dataDir)/"
fasta = "$(dataDir)/"

# IO
config = CpelNano.CpelNanoConfig();
config.out_prefix = "prefix"
config.max_em_iters = 100
config.out_dir = outdir
config.max_em_init = 10
config.min_cov = 5.0
config.pe = true

# Make true if you want to activate informME mode
config.informme_mode = false

# Analyze BAM file
CpelNano.analyze_bam(bam,fasta,config)


text:
using BinDeps

basedir = dirname(@__FILE__)

program = "exiftool"

if is_unix()
    file = "Image-ExifTool-10.49"
    extension = ".tar.gz"
    binary_name = target = program
end

if is_windows()
    file = "exiftool-10.49"
    extension = ".zip"
    binary_name = "$program.exe"
    target = "exiftool(-k).exe"
end

filename = file*extension
url = "http://www.sno.phy.queensu.ca/~phil/exiftool/$filename"

if is_unix()
    run(
        @build_steps begin
        FileDownloader(url, joinpath(basedir, "downloads", filename))
        CreateDirectory(joinpath(basedir, "src"))
        FileUnpacker(joinpath(basedir, "downloads", filename), joinpath(basedir, "src"), "")
        end
       )
    mv(joinpath(basedir, "src", file), joinpath(basedir, "src", "exiftool"), remove_destination = true)
end

if is_windows()
    run(
        @build_steps begin
        FileDownloader(url, joinpath(basedir, "downloads", filename))
        CreateDirectory(joinpath(basedir, "src"))
        FileUnpacker(joinpath(basedir, "downloads", filename), joinpath(basedir, "src"), target)
        CreateDirectory(joinpath(basedir, "src", "exiftool"))
        end
       )
    mv(joinpath(basedir, "src", target), joinpath("src", "exiftool", binary_name), remove_destination = true)
end


text:
#========================================================================================#
#	Laboratory 10
#
# Working with arrays.
#
# Author: Niall Palfreyman, 22/03/2022
#========================================================================================#
[
	Activity(
		"""
		Arrays play a very fundamental role in Julia, so we need to learn how to deal with them
		efficiently. In this laboratory we learn how to index into arrays in Julia - that is,
		how to use indices to access and manipulate the entries in an array. We will also look
		at various ways of applying code to the elements in an array.

		Julia offers four kinds of indexing: SUBSCRIPTING, LINEAR indexing, MULTIPLE indexing
		and LOGICAL indexing. All four are extremely useful for scientific programming. To
		investigate these forms of indexing, first create a (4x4) matrix by reshaping a range:

			m = reshape(1:16,(4,4))

		First notice the order in which Julia has placed the numbers 1:16. Do they run down the
		columns or along the rows?
		""",
		"",
		x -> occursin("col",lowercase(x))
	),
	Activity(
		"""
		So: SUBSCRIPTING! To access a matrix element with subscripts, we use square brackets
		enclosing two subscripts i (row) and j (column) like this: m[i,j]. In Julia, the row
		index ALWAYS comes before the column index! Display the element m[3,4[] now, then tell
		me its value:
		""",
		"Just enter m[3,4] at the Julia prompt",
		x -> x==15
	),
	Activity(
		"""
		We have already used the colon operator to create ranges of numbers, for example 1:3 can
		be collected into a vector [1,2,3]. We can also use ranges to select regions of an array.
		Use m[1:3,2:4] to display the top-right (3×3) region of m, then tell me the result
		(Remember you can use reply(ans) to tell me the result of the previous calculation):
		""",
		"",
		x -> x == [5 9 13;6 10 14;7 11 15]
	),
	Activity(
		"""
		1:2:7 creates a range of numbers from 1 to 7 in steps of 2: [1,3,5,7]. Reshape 1:81 into a
		(9×9) matrix, then extract from it the (3×3) matrix obtained by taking the first, then
		every third, row and column of your (9×9) matrix. Then tell me the result:
		""",
		"",
		x -> x == [1 28 55;4 31 58;7 34 61]
	),
	Activity(
		"""
		We can select a whole row or column by writing ":". Display and tell me the third row of m:
		""",
		"",
		x -> x == [3 7 11 15]
	),
	Activity(
		"""
		To assign values to regions of a matrix, we select the desired matrix range and assign to it
		a value. However, we should first be aware of something important. Tell me the type of m:
		""",
		"",
		x -> x <: Base.ReshapedArray
	),
	Activity(
		"""
		Notice that m only LOOKS like an ordinary matrix, but is actually a reshaped range, so Julia
		will not let us change its value by writing to it. To see this, try out the following:

			m[2:3,3:4] .= 1

		Then tell me the name of the function that the error message recommends we use:
		""",
		"",
		x -> occursin("collect",lowercase(x))
	),
	Activity(
		"""
		Let's take the error message's advice:

			m = collect(m)
			m[2:3,3:4] .= 1

		Now tell me the new value of m:
		""",
		"",
		x -> x == [1 5 9 13;2 6 1 1;3 7 1 1;4 8 12 16]
	),
	Activity(
		"""
		And what is the type of m now?
		""",
		"",
		x -> x <: Matrix
	),
	Activity(
		"""
		What kind of error do you get if you try to add an element outside the size of the array m:

			m[5,1] = 1
		""",
		"",
		x -> occursin("bounds",lowercase(x))
	),
	Activity(
		"""
		Although we can't add elements outside m's bounds, we can extend the size of m by adding
		new rows or columns: Blank ' ' adds columns (hcat: horizontal concatenation), and
		semicolon ';' adds rows (vcat: vertical concatenation). Tell me the result of this line:

			p = ones(4)
			pp = [m p]
		""",
		"",
		x -> size(x) == (4,5)
	),
	Activity(
		"""
		And what is the result of this line of code?

			q = ones(4)'
			qq = [m;q]
		""",
		"",
		x -> size(x) == (5,4)
	),
	Activity(
		"""
		OK, now let's look at LINEAR indexing! Internally, Julia represents all arrays as vectors =
		it simply makes this vector look like a matrix to us. Create the following matrix:

			A = [1 2 3;4 5 6]

		If we just enter A at the Julia prompt, we see a (2×3) matrix, but if you enter A[:], you
		will see a list of all elements of A in a linear order. Does this order first run down the
		columns or along the rows?
		""",
		"",
		x -> occursin("col",lowercase(x))
	),
	Activity(
		"""
		We can use the linear ordering to index the elements of A. What is the value of A[5]?
		""",
		"",
		x -> x==3
	),
	Activity(
		"""
		We can also use linear indexing to change the elements of A. What is the contents of A
		after entering this line of code?

			A[5] = 99
		""",
		"",
		x -> [1 2 99;4 5 6]
	),
	Activity(
		"""
		Arrays can have MULTIPLE indices: Vectors have one index and Matrices have two, but we
		can have 3, 4, 5, ...! Of course we need to be a bit careful with the size of these arrays:
		A size (10) Vector contains 10 elements; a size (10,10) Matrix contains 100 elements. How
		many elements will there be in a size (10,10,10,10) array?
		""",
		"",
		x -> x == 1e4
	),
	Activity(
		"""
		We can create multiply indexed arrays using zeros(), ones(), rand(), randn() or reshape().
		A 3-dimensional array might represent 3-dimensional data, such as a chemical concentration
		at various locations in a cell, or it might represent the time-series of elements of a
		matrix A(t). In this case, A[3,5,4] might represent the element [3,5] at time t=4 of the
		time-series.

		Enter R = rand(2,3,4); at the Julia prompt, and study the following expressions:

			R[1,:,:]
			R[]:,1,:]
			R[]:,:,1]
			size(R)

		How many elements does R contain?
		""",
		"Use the length function",
		x -> x == 24
	),
	Activity(
		"""
		Here is a (4,4) magic square:

			A = [1 15 14 4;10 11 8 5;7 6 9 12;16 2 3 13]

		Use Julia's sum() function to find the sum of A's elements along any row, column or
		diagonal:
		""",
		"",
		x -> x==34
	),
	Activity(
		"""
		We can generate new magic squares by swapping columns of an existing magic square. Study
		the matrix B = A(:,[1,3,2,4]). Is this also a magic square (remember that rows, columns and
		diagonals of a magic square must all sum to the same number)? Which substructures of A were
		permuted by using the Vector [1,3,2,4] as a permutation index?
		""",
		"",
		x -> occursin("col",lowercase(x))
	),
	Activity(
		"""
		Use a permutation index to swap ROWS 2 and 3 of A, and tell me the result:
		""",
		"",
		x -> x == [1 15 14 4;7 6 9 12;10 11 8 5;16 2 3 13]
	),
	Activity(
		"""
		LOGICAL indexing uses an array of logical values as an index to another array. For example,
		suppose we want to reduce to zero all elements less than 3 in the Vector v = [1,2,3,4,5].
		One way would be to access each individual element in v, check whether it is less than 3,
		and set it to zero. But this would be very inefficient, because accessing each individual
		element costs time. Instead, we can use LOGICAL indexing to change the entire vector v in
		one sweep ...

		First, create the Vector v. Then tell me the result of entering this line:

			d = (v .< 3)
		""",
		"",
		x -> x == [1,1,0,0,0]
	),
	Activity(
		"""
		What is the type of the elements of d?
		""",
		"",
		x -> x == Bool
	),
	Activity(
		"""
		Notice that d is a Vector of Bools with the same length as v. We can use d to index
		elements of v:

			v[d]

		How many elements of v does d pick out?
		""",
		"",
		x -> x == 2
	),
	Activity(
		"""
		Now enter this line:

			v[d] .= 0

		You will see that this zeros out all elements of v that are less than 3. If you recreate
		v, you can even condense this entire process into one step:

			v = collect(1:5)
			v[v.<3] .= 0

		Now set to zero all numbers in the vector -5:2:5 which are greater than 2, and tell me
		the number of non-zero elements in your result:
		""",
		"",
		x -> x == 4
	),
	Activity(
		"""
		Before continuing, let's look at several different ways to apply some code to all elements
		in an array. If the code is just one function, this is easy: we just use the broadcast
		dot (.):

			v = 1:7
			isodd.(v)
		""",
		"",
		x -> x == isodd.(1:7)
	),
	Activity(
		"""
		If the code is a little more complicated, we might map() an anonymous function over
		the array:

			map( x->(sin(x) >= 0), v)
		""",
		"",
		x -> x == map( y->(sin(y) >= 0), 1:7)
	),
	Activity(
		"""
		And finally, if the code is particularly complicated, we can use a do statement that
		allows us to define a complicated mapping over all elements of the array:

			map(v) do x
				if x < 4
					isodd(x)
				else
					iseven(x)
				end
			end
		""",
		"",
		x -> x == (map(1:7) do y if y<4 isodd(y) else iseven(y) end end)
	),
	Activity(
		"""
		OK, now the last two activities in this laboratory give you practice in applying indexing,
		broadcasting and mapping to problems that often arise in signal-processing. Have fun! :)

		Use logical indexing to generate a list of all odd multiples of 3 in the range 1:50 :
		""",
		"Use isodd(), rem(), &, and remember to use broadcasts (.) and brackets",
		x -> x == [3,9,15,21,27,33,39,45]
	),
	Activity(
		"""
		This function decides whether or not its argument n is a prime number:

			function isprime(n::Int)
				if n < 2 return false end
				if n in 2:3 return true end
				for i in 2:floor(Int,sqrt(n))
					if rem(n,i) == 0 return false end
				end
				true
			end

		Use the isprime function to generate a list of twenty numbers from 1 to 20, in which
		all prime numbers AND all multiples of 3 are zeroed out:
		""",
		"Use rem(), isprime() and |",
		x -> x == [1,0,0,4,0,0,0,8,0,10,0,0,0,14,0,16,0,0,0,20]
	),
]


text:
mutable struct Annotation
    cascade::Cascade
    label::Vector{Label}
end

Annotation( ; cascade, label) = Annotation(cascade, label)


text:
using EconPDEs

Base.@kwdef struct GarleanuPanageasModel

  # utility function
  γA::Float64 = 1.5
  ψA::Float64 = 0.7
  γB::Float64 = 10.0
  ψB::Float64 = 0.05
  ρ::Float64 = 0.001
  δ::Float64 = 0.02

  # proportion a
  νA::Float64 = 0.01

  # consumption
  μ::Float64 = 0.02
  σ::Float64 = 0.041

  # earning function
  B1::Float64 = 30.72
  δ1::Float64 = 0.0525
  B2::Float64 = -30.29
  δ2::Float64 = 0.0611
  ω::Float64 = 0.92
end

function (m::GarleanuPanageasModel)(state::NamedTuple, y::NamedTuple)
  (; γA, ψA, γB, ψB, ρ, δ, νA, μ, σ, B1, δ1, B2, δ2, ω) = m
  (; x) = state
  (; pA, pAx_up, pAx_down, pAxx, pB, pBx_up, pBx_down, pBxx, ϕ1, ϕ1x_up, ϕ1x_down, ϕ1xx, ϕ2, ϕ2x_up, ϕ2x_down, ϕ2xx) = y

  scale = δ / (δ + δ1) * B1 + δ / (δ + δ2) * B2
  B1 = B1 / scale
  B2 = B2 / scale
  # pA is wealth / consumption ratio of agent A
  # pB is wealth / consumption ratio of agent B
  # ϕ1 is value of claim that promises B_1ωexp(-(δ+δ1)(s-t))C_s/C_t for s ≥ t
  # ϕ2 is value of claim that promises B_2ωexp(-(δ+δ2)(s-t))C_s/C_t for s ≥ t
  # Market price of risk κ
  pAx, pBx, ϕ1x, ϕ2x = pAx_up, pBx_up, ϕ1x_up, ϕ2x_up
  iter = 0
  @label start
  Γ = 1 / (x / γA + (1 - x) / γB)
  p = x * pA + (1 - x) * pB
  σx = σ * x * (Γ / γA - 1) / (1 + Γ * x * (1 - x) / (γA * γB) * ((1 - γB * ψB) / (ψB - 1) * (pBx / pB) - (1 - γA * ψA) / (ψA - 1) * (pAx / pA)))
  σpA = pAx / pA * σx
  σpB = pBx / pB * σx
  σϕ1 = ϕ1x / ϕ1 * σx
  σϕ2 = ϕ2x / ϕ2 * σx
  κ = Γ * (σ - x * (1 - γA * ψA) / (γA * (ψA - 1)) * σpA - (1 - x) * (1 - γB * ψB) / (γB * (ψB - 1)) * σpB)
  σCA = κ / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * σpA
  σCB = κ / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * σpB
  # A.16 Equation in Garleanu Panageas has a typo
  mcA = κ^2 * (1 + ψA) / (2 * γA) + (1 - ψA * γA) / (γA * (ψA - 1)) * κ * σpA - (1 - γA * ψA) / (2 * γA * (ψA - 1)) * σpA^2
  mcB = κ^2 * (1 + ψB) / (2 * γB) + (1 - ψB * γB) / (γB * (ψB - 1)) * κ * σpB - (1 - γB * ψB) / (2 * γB * (ψB - 1)) * σpB^2
  r =  ρ + 1 / (ψA * x  + ψB * (1 - x))  * (μ - x * mcA - (1 - x) * mcB - δ * ((νA / pA + (1 - νA) / pB) * (ϕ1 + ϕ2) - 1))
  μCA = ψA * (r - ρ) + mcA
  μCB = ψB * (r - ρ) + mcB
  μx = x * (μCA - δ - μ) + δ * νA / pA * (ϕ1 + ϕ2) - σ * σx
  if (iter == 0) & (μx <= 0)
    iter += 1
    pAx, pBx, ϕ1x, ϕ2x = pAx_down, pBx_down, ϕ1x_down, ϕ2x_down
    @goto start
  end

  μpA = pAx / pA * μx + 0.5 * pAxx / pA * σx^2
  μpB = pBx / pB * μx + 0.5 * pBxx / pB * σx^2
  μϕ1 = ϕ1x / ϕ1 * μx + 0.5 * ϕ1xx / ϕ1 * σx^2
  μϕ2 = ϕ2x / ϕ2 * μx + 0.5 * ϕ2xx / ϕ2 * σx^2

  # Market Pricing
  pAt = - pA * (1 / pA + (μCA - δ) + μpA + σCA * σpA - r - κ * (σpA + σCA))
  pBt = - pB * (1 / pB + (μCB - δ) + μpB + σCB * σpB - r - κ * (σpB + σCB))
  ϕ1t = - ϕ1 * (B1 * ω / ϕ1 + (μ - δ - δ1) + μϕ1 + σ * σϕ1 - r - κ * (σϕ1 + σ))
  ϕ2t = - ϕ2 * (B2 * ω / ϕ2 + (μ - δ - δ2) + μϕ2 + σ * σϕ2 - r - κ * (σϕ2 + σ))

  return (; pAt, pBt, ϕ1t, ϕ2t)
end

m = GarleanuPanageasModel()
xn = 200
stategrid = OrderedDict(:x => range(0.0, 1.0, length = xn))
yend = OrderedDict(:pA => ones(xn), :pB => ones(xn), :ϕ1 => ones(xn), :ϕ2 => ones(xn))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5


text:
#==============================================================================#
# SageMakerRuntime.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/runtime.sagemaker-2017-05-13.normal.json
#==============================================================================#

__precompile__()

module SageMakerRuntime

using AWSCore

"""
    using AWSSDK.SageMakerRuntime.invoke_endpoint
    invoke_endpoint([::AWSConfig], arguments::Dict)
    invoke_endpoint([::AWSConfig]; EndpointName=, Body=, <keyword arguments>)

    using AWSCore.Services.runtime_sagemaker
    runtime_sagemaker([::AWSConfig], "POST", "/endpoints/{EndpointName}/invocations", arguments::Dict)
    runtime_sagemaker([::AWSConfig], "POST", "/endpoints/{EndpointName}/invocations", EndpointName=, Body=, <keyword arguments>)

# InvokeEndpoint Operation

After you deploy a model into production using Amazon SageMaker hosting services, your client applications use this API to get inferences from the model hosted at the specified endpoint.

For an overview of Amazon SageMaker, see [How It Works](http://docs.aws.amazon.com/sagemaker/latest/dg/how-it-works.html)

Amazon SageMaker strips all POST headers except those supported by the API. Amazon SageMaker might add additional headers. You should not rely on the behavior of headers outside those enumerated in the request syntax.

# Arguments

## `EndpointName = ::String` -- *Required*
The name of the endpoint that you specified when you created the endpoint using the [CreateEndpoint](http://docs.aws.amazon.com/sagemaker/latest/dg/API_CreateEndpoint.html) API.

## `Body = blob` -- *Required*
Provides input data, in the format specified in the `ContentType` request header. Amazon SageMaker passes all of the data in the body to the model.

## `*header:* Content-Type = ::String`
The MIME type of the input data in the request body.

## `*header:* Accept = ::String`
The desired MIME type of the inference in the response.

# Returns

`InvokeEndpointOutput`

# Exceptions

`InternalFailure`, `ServiceUnavailable`, `ValidationError` or `ModelError`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/runtime.sagemaker-2017-05-13/InvokeEndpoint)
"""
@inline invoke_endpoint(aws::AWSConfig=default_aws_config(); args...) = invoke_endpoint(aws, args)

@inline invoke_endpoint(aws::AWSConfig, args) = AWSCore.Services.runtime_sagemaker(aws, "POST", "/endpoints/{EndpointName}/invocations", args)

@inline invoke_endpoint(args) = invoke_endpoint(default_aws_config(), args)

end # module SageMakerRuntime

#==============================================================================#
# End of file
#==============================================================================#


text:
export
# From external modules
  # LinearAlgebra
  tr,

# Modules
  LinearAlgebra,
  NDTensors,

# NDTensors module
  # Types
  Block,
  Spectrum,
  # Methods
  eigs,
  entropy,
  truncerror,
  # Deprecated
  addblock!,

# ITensors.jl
  index_id_rng,

# argsdict/argsdict.jl
  argsdict,

# decomp.jl
  eigen,
  factorize,
  polar,
  qr,
  svd,

# global_variables.jl
  # Methods

  # Macros
  @disable_warn_order,
  @reset_warn_order,
  @set_warn_order,

# index.jl
  # Types
  Index,
  IndexVal,
  # Methods
  dag,
  dim,
  dir,
  hasid,
  hasind,
  hasplev,
  hasqns,
  id,
  ind,
  isindequal,
  noprime,
  plev,
  prime,
  removetags,
  removeqns,
  replacetags,
  replacetags!,
  setdir,
  setprime,
  setspace,
  settags,
  sim,
  space,
  splitblocks,
  tags,
  val,

# indexset.jl
  # Types
  IndexSet,
  Order,
  # Methods
  allhastags,
  anyhastags,
  dims,
  firstintersect,
  firstsetdiff,
  getfirst,
  mapprime,
  maxdim,
  mindim,
  permute,
  pop,
  popfirst,
  push,
  pushfirst,
  replaceind,
  replaceinds,
  replaceprime,
  swapinds,
  setindex,
  swapind,
  swapinds,
  swapprime,
  swaptags,

# itensor.jl
  # Types
  ITensor,
  # Methods
  ⊙,
  addtags!,
  apply,
  array,
  axpy!,
  blockoffsets,
  combinedind,
  combiner,
  commonind,
  commoninds,
  complex!,
  delta,
  dense,
  δ,
  diagITensor,
  dot,
  eachnzblock,
  firstind,
  filterinds,
  hadamard_product,
  hascommoninds,
  hasind,
  hasinds,
  hassameinds,
  ind,
  inds,
  insertblock!,
  itensor,
  mul!,
  matrix,
  mapprime!,
  noncommonind,
  noncommoninds,
  norm,
  normalize!,
  noprime!,
  nnzblocks,
  nzblocks,
  nzblock,
  nnz,
  order,
  permute,
  prime!,
  product,
  randn!,
  randomITensor,
  removetags!,
  replacetags!,
  replaceind!,
  replaceinds!,
  swapinds!,
  rmul!,
  scale!,
  scalar,
  setelt,
  store,
  setprime!,
  swapprime!,
  settags!,
  swaptags!,
  uniqueinds,
  uniqueind,
  unioninds,
  unionind,
  vector,
  emptyITensor,

# iterativesolvers.jl
  davidson,

# not.jl
  not,

# readwrite.jl
  readcpp,

# tagset.jl
  # Types
  TagSet,
  # Macros
  @ts_str,
  # Methods
  addtags,
  hastags,

# mps/dmrg.jl
  dmrg,

# mps/abstractmps.jl
  # Macros
  @preserve_ortho,
  # Methods
  add,
  contract,
  common_siteind,
  common_siteinds,
  dag!,
  findfirstsiteind,
  findfirstsiteinds,
  findsite,
  findsites,
  firstsiteind,
  firstsiteinds,
  logdot,
  loginner,
  lognorm,
  movesite,
  movesites,
  ortho_lims,
  siteinds,

# mps/mpo.jl
  # Types
  MPO,
  # Methods
  error_contract,
  maxlinkdim,
  orthogonalize,
  orthogonalize!,
  randomMPO,
  truncate,
  truncate!,
  unique_siteind,
  unique_siteinds,

# mps/mps.jl
  # Types
  MPS,
  # Methods
  ⋅,
  inner,
  isortho,
  linkdim,
  linkdims,
  linkind,
  linkinds,
  productMPS,
  randomMPS,
  replacebond,
  replacebond!,
  sample,
  sample!,
  siteind,
  siteinds,
  replace_siteinds!,
  replace_siteinds,
  swapbondsites,
  totalqn,

# mps/observer.jl
  # Types
  AbstractObserver,
  DMRGObserver,
  DMRGMeasurement,
  NoObserver,
  # Methods
  checkdone!,
  energies,
  measure!,
  measurements,
  truncerrors,

# mps/projmpo.jl
  ProjMPO,
  lproj,
  product,
  rproj,
  noiseterm,
  position!,

# mps/projmposum.jl
  ProjMPOSum,

# mps/projmpo_mps.jl
  ProjMPO_MPS,

# mps/sweeps.jl
  Sweeps,
  cutoff,
  cutoff!,
  get_cutoffs,
  get_maxdims,
  get_mindims,
  get_noises,
  maxdim,
  maxdim!,
  mindim,
  mindim!,
  noise,
  noise!,
  nsweep,
  sweepnext,

# physics/autompo.jl
  AutoMPO,
  add!,

# physics/lattices.jl
  Lattice,
  LatticeBond,
  square_lattice,
  triangular_lattice,

# physics/sitetype.jl
  SiteType,
  @SiteType_str,
  StateName,
  @StateName_str,
  op,
  ops,
  OpName,
  @OpName_str,
  state,
  TagType,
  @TagType_str,
  has_fermion_string,

# qn/qn.jl
  # Types
  QN,
  # Methods
  isactive,
  isfermionic,
  modulus,
  val,

# qn/qnindex.jl
  flux,
  hasqns,
  nblocks,
  qn