text:
module QuantumSystems

using LinearAlgebra: dot, norm
using StaticArrays: SVector
using Printf: @printf, @sprintf
using ..QuantumOperators: ID, OperatorProd, Operator, Operators, LaTeX, latexformat
using ..Spatials: AbstractPID, AbstractBond, Point, Bond, rcoord, pidtype
using ..DegreesOfFreedom: SimpleIID, CompositeIID, SimpleInternal, CompositeInternal
using ..DegreesOfFreedom: Index, CompositeOID, OID, OIDToTuple, Hilbert, Table
using ..DegreesOfFreedom: Subscripts, Subscript, IIDSpace, subscriptexpr, wildcard, diagonal
using ..DegreesOfFreedom: Coupling, Couplings, @couplings, couplinginternals, Term, TermCouplings, TermAmplitude, TermModulate
using ...Essentials: kind, dtype
using ...Interfaces: decompose, dimension
using ...Prerequisites: Float, atol, rtol, delta, decimaltostr, concatenate
using ...Prerequisites.Traits: rawtype, getcontent
using ...Prerequisites.VectorSpaces: VectorSpace, VectorSpaceStyle, VectorSpaceCartesian

import ..QuantumOperators: script, latexname, ishermitian, optype
import ..DegreesOfFreedom: statistics, couplingcenters, abbr, termfactor
import ...Interfaces: rank, ⊗, ⋅, expand, expand!, permute
import ...Prerequisites.VectorSpaces: shape

export majorana, annihilation, creation, flatex, blatex
export FID, Fock, FockCoupling, Onsite, Hopping, Pairing, Hubbard, InterOrbitalInterSpin, InterOrbitalIntraSpin, SpinFlip, PairHopping, Coulomb, FockTerm
export isnormalordered, @σ⁰_str, @σˣ_str, @σʸ_str, @σᶻ_str, @σ⁺_str, @σ⁻_str, @fc_str

export slatex
export SID, Spin, SpinCoupling, SpinTerm
export totalspin, @heisenberg_str, @ising_str, @gamma_str, @dm_str, @sˣ_str, @sʸ_str, @sᶻ_str, @sc_str

export nlatex
export NID, Phonon, PhononCoupling, PhononKinetic, PhononPotential, PhononTerm
export @kinetic_str, @potential_str

# Canonical fermionic/bosonic systems and hardcore bosonic systems
"""
    majorana

Indicate that the nambu index is majorana.
"""
const majorana = 0

"""
    annihilation

Indicate that the nambu index is annihilation.
"""
const annihilation = 1

"""
    creation

Indicate that the nambu index is creation.
"""
const creation = 2

"""
    FID{T, O<:Union{Int, Symbol}, S<:Union{Int, Symbol}, N<:Union{Int, Symbol}} <: SimpleIID

The Fock id.
"""
struct FID{T, O<:Union{Int, Symbol}, S<:Union{Int, Symbol}, N<:Union{Int, Symbol}} <: SimpleIID
    orbital::O
    spin::S
    nambu::N
    function FID{T}(orbital::Union{Int, Symbol}, spin::Union{Int, Symbol}, nambu::Union{Int, Symbol}) where T
        @assert T∈(:f, :b, wildcard) "FID error: wrong statistics."
        isa(nambu, Int) && @assert nambu∈(0, 1, 2) "FID error: wrong input nambu($nambu)."
        new{T, typeof(orbital), typeof(spin), typeof(nambu)}(orbital, spin, nambu)
    end
end
Base.show(io::IO, fid::FID) = @printf io "FID{%s}(%s)" repr(statistics(fid)) join(ntuple(i->getfield(fid, i), Val(fieldcount(typeof(fid)))), ", ")
@inline Base.adjoint(fid::FID{T, <:Union{Int, Symbol}, <:Union{Int, Symbol}, Int}) where T = FID{T}(fid.orbital, fid.spin, fid.nambu==0 ? 0 : 3-fid.nambu)
@inline @generated function Base.replace(fid::FID; kwargs...)
    exprs = [:(get(kwargs, $name, getfield(fid, $name))) for name in QuoteNode.(fieldnames(fid))]
    return :(rawtype(typeof(fid)){statistics(fid)}($(exprs...)))
end
@inline statistics(::Type{<:FID{T}}) where T = T
@inline FID(iid::FID, ::CompositeInternal) = iid

"""
    FID{T}(; orbital::Union{Int, Symbol}=1, spin::Union{Int, Symbol}=1, nambu::Union{Int, Symbol}=annihilation) where T

Create a Fock id.
"""
@inline FID{T}(; orbital::Union{Int, Symbol}=1, spin::Union{Int, Symbol}=1, nambu::Union{Int, Symbol}=annihilation) where T = FID{T}(orbital, spin, nambu)

"""
    Fock{T} <: SimpleInternal{FID{T, Int, Int, Int}}

The Fock internal degrees of freedom.
"""
struct Fock{T} <: SimpleInternal{FID{T, Int, Int, Int}}
    norbital::Int
    nspin::Int
    nnambu::Int
    function Fock{T}(norbital::Int, nspin::Int, nnambu::Int) where T
        @assert T∈(:f, :b) "Fock error: wrong statistics."
        @assert nnambu∈(1, 2) "Fock error: wrong input nnambu($nnambu)."
        new{T}(norbital, nspin, nnambu)
    end
end
@inline Base.eltype(::Type{Fock}) = (FID{T, Int, Int, Int} where T)
@inline shape(fock::Fock) = (1:fock.norbital, 1:fock.nspin, fock.nnambu==1 ? (0:0) : (1:2))
@inline Base.CartesianIndex(fid::FID{T}, fock::Fock{T}) where T = CartesianIndex(fid.orbital, fid.spin, fid.nambu)
@inline FID(index::CartesianIndex{3}, fock::Fock) = FID{statistics(fock)}(index[1], index[2], index[3])
Base.summary(io::IO, fock::Fock) = @printf io "%s-element Fock{%s}" length(fock) repr(statistics(fock))
function Base.show(io::IO, fock::Fock)
    @printf io "%s{%s}(%s)" fock|>typeof|>nameof repr(statistics(fock)) join(("$name=$(getfield(fock, name))" for name in fock|>typeof|>fieldnames), ", ")
end
@inline Base.match(::Type{<:FID{wildcard}}, ::Type{<:Fock{T}}) where T = true
@inline Base.match(::Type{<:FID{T}}, ::Type{<:Fock{T}}) where T = true
@inline Base.match(::Type{<:FID{T₁}}, ::Type{<:Fock{T₂}}) where {T₁, T₂} = false

@inline function CompositeIID(iidspace::IIDSpace{I, V}) where {I<:CompositeIID{<:Tuple{Vararg{FID}}}, V<:CompositeInternal{:⊗, <:Tuple{Vararg{Fock}}}}
    @assert rank(I)==rank(V) "CompositeIID error: mismatched composite iid and composite internal space."
    return CompositeIID(map(
        (order, fid, fock)->FID{statistics(fid)}(fid.orbital, fid.spin, fidnambu(Val(kind(iidspace)), fid.nambu, fock.nnambu, order)),
        ntuple(i->i, Val(rank(I))),
        iidspace.iid.contents,
        iidspace.internal.contents
        ))
end
@inline fidnambu(::Val, ::Symbol, nnambu::Int, order::Int) = nnambu==1 ? majorana : order%2==1 ? creation : annihilation
@inline fidnambu(::Val, nambu::Int, nnambu::Int, ::Int) = ((@assert nnambu==1 ? nambu==0 : 0<nambu<nambu+1 "fidnambu error: nambu out of range."); nambu)

@inline function shape(iidspace::IIDSpace{<:FID, <:Fock})
    @assert isa(iidspace.iid.nambu, Int) "shape error: nambu not determined."
    obrange = fidshape(iidspace.iid.orbital, iidspace.internal.norbital)
    sprange = fidshape(iidspace.iid.spin, iidspace.internal.nspin)
    phrange = iidspace.iid.nambu:iidspace.iid.nambu
    return (obrange, sprange, phrange)
end
@inline fidshape(::Symbol, n::Int) = 1:n
@inline fidshape(v::Int, n::Int) where field = ((@assert 0<v<n+1 "shape error: $field out of range."); v:v)

"""
    Fock{T}(; norbital::Int=1, nspin::Int=2, nnambu::Int=2) where T

Construct a Fock degrees of freedom.
"""
@inline Fock{T}(; norbital::Int=1, nspin::Int=2, nnambu::Int=2) where T = Fock{T}(norbital, nspin, nnambu)

"""
    script(::Val{:site}, index::Index{<:AbstractPID, <:FID}; kwargs...) -> Int
    script(::Val{:orbital}, index::Index{<:AbstractPID, <:FID}; kwargs...) -> Int
    script(::Val{:spint}, index::Index{<:AbstractPID, <:FID}; kwargs...) -> Int
    script(::Val{:spsym}, index::Index{<:AbstractPID, <:FID}; kwargs...) -> String
    script(::Val{:nambu}, index::Index{<:AbstractPID, <:FID}; kwargs...) -> String

Get the required script of a Fock index.
"""
@inline script(::Val{:site}, index::Index{<:AbstractPID, <:FID}; kwargs...) = index.pid.site
@inline script(::Val{:orbital}, index::Index{<:AbstractPID, <:FID}; kwargs...) = index.iid.orbital
@inline script(::Val{:spint}, index::Index{<:AbstractPID, <:FID}; kwargs...) = index.iid.spin
@inline script(::Val{:spsym}, index::Index{<:AbstractPID, <:FID}; kwargs...) = index.iid.spin==1 ? "↓" : index.iid.spin==2 ? "↑" : error("script error: wrong spin.")
@inline script(::Val{:nambu}, index::Index{<:AbstractPID, <:FID}; kwargs...) = index.iid.nambu==creation ? "\\dagger" : ""

"""
    flatex

The default LaTeX format for a fermionic oid.
"""
const flatex = LaTeX{(:nambu,), (:site, :orbital, :spsym)}('c')
@inline latexname(::Type{<:Index{<:AbstractPID, <:FID{:f}}}) = Symbol("Index{AbstractPID, FID{:f}}")
@inline latexname(::Type{<:CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}}) = Symbol("CompositeOID{Index{AbstractPID, FID{:f}}}")
latexformat(Index{<:AbstractPID, <:FID{:f}}, flatex)
latexformat(CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}, flatex)

"""
    blatex

The default LaTeX format for a bosonic oid.
"""
const blatex = LaTeX{(:nambu,), (:site, :orbital, :spsym)}('b')
@inline latexname(::Type{<:Index{<:AbstractPID, <:FID{:b}}}) = Symbol("Index{AbstractPID, FID{:b}}")
@inline latexname(::Type{<:CompositeOID{<:Index{<:AbstractPID, <:FID{:b}}}}) = Symbol("CompositeOID{Index{AbstractPID, FID{:b}}}")
latexformat(Index{<:AbstractPID, <:FID{:b}}, blatex)
latexformat(CompositeOID{<:Index{<:AbstractPID, <:FID{:b}}}, blatex)

"""
    angle(id::OID{<:Index{<:AbstractPID, <:FID}}, vectors::AbstractVector{<:AbstractVector{<:Number}}, values::AbstractVector{<:Number}) -> Complex{<:Number}

Get the twist phase corresponding to a Fock oid.
"""
function Base.angle(id::OID{<:Index{<:AbstractPID, <:FID}}, vectors::AbstractVector{<:AbstractVector{<:Number}}, values::AbstractVector{<:Number})
    datatype = promote_type(eltype(values), Float)
    phase = length(vectors)==1 ? 2*convert(datatype, pi)*dot(decompose(id.icoord, vectors[1]), values) :
            length(vectors)==2 ? 2*convert(datatype, pi)*dot(decompose(id.icoord, vectors[1], vectors[2]), values) :
            length(vectors)==3 ? 2*convert(datatype, pi)*dot(decompose(id.icoord, vectors[1], vectors[2], vectors[3]), values) :
            error("angle error: not supported number of input basis vectors.")
    return id.index.iid.nambu==annihilation ? phase : id.index.iid.nambu==creation ? -phase : error("angle error: not supported Fock index.")
end

"""
    isnormalordered(opt::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID}}}}) -> Bool

Judge whether an operator is normal ordered.
"""
function isnormalordered(opt::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID}}}})
    flag = true
    for i = 1:rank(opt)
        flag && (opt.id[i].index.iid.nambu == annihilation) && (flag = false)
        flag || (opt.id[i].index.iid.nambu == creation) && return false
    end
    return true
end

"""
    permute(id₁::OID{<:Index{<:AbstractPID, <:FID{:f}}}, id₂::OID{<:Index{<:AbstractPID, <:FID{:f}}}) -> Tuple{Vararg{Operator}}

Permute two fermionic oids and get the result.
"""
function permute(id₁::OID{<:Index{<:AbstractPID, <:FID{:f}}}, id₂::OID{<:Index{<:AbstractPID, <:FID{:f}}})
    @assert id₁.index≠id₂.index || id₁.rcoord≠id₂.rcoord || id₁.icoord≠id₂.icoord "permute error: permuted ids should not be equal to each other."
    if id₁.index'==id₂.index && id₁.rcoord==id₂.rcoord && id₁.icoord==id₂.icoord
        return (Operator(1), Operator(-1, ID(id₂, id₁)))
    else
        return (Operator(-1, ID(id₂, id₁)),)
    end
end

"""
    *(  f1::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}}},
        f2::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}}}
        ) -> Union{Nothing, Operator}

Get the multiplication of two fermionic Fock operators.
"""
@inline function Base.:*(
        f1::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}}},
        f2::Operator{<:Number, <:ID{CompositeOID{<:Index{<:AbstractPID, <:FID{:f}}}}}
        )
    rank(f1)>0 && rank(f2)>0 && f1.id[end]==f2.id[1] && return nothing
    return invoke(*, Tuple{OperatorProd, OperatorProd}, f1, f2)
end

"""
    permute(id₁::OID{<:Index{<:AbstractPID, <:FID{:b}}}, id₂::OID{<:Index{<:AbstractPID, <:FID{:b}}}) -> Tuple{Vararg{Operator}}

Permute two bosonic oids and get the result.
"""
function permute(id₁::OID{<:Index{<:AbstractPID, <:FID{:b}}}, id₂::OID{<:Index{<:AbstractPID, <:FID{:b}}})
    @assert id₁.index≠id₂.index || id₁.rcoord≠id₂.rcoord || id₁.icoord≠id₂.icoord "permute error: permuted ids should not be equal to each other."
    if id₁.index'==id₂.index && id₁.rcoord==id₂.rcoord && id₁.icoord==id₂.icoord
        if id₁.index.iid.nambu == creation
            return (Operator(1), Operator(1, ID(id₂, id₁)))
        else
            return (Operator(-1), Operator(1, ID(id₂, id₁)))
        end
    else
        return (Operator(1, ID(id₂, id₁)),)
    end
end

"""
    FockCoupling(value::Number,
        orbitals::Subscript{<:NTuple{N, Union{Int, Symbol}}},
        spins::Subscript{<:NTuple{N, Union{Int, Symbol}}},
        nambus::Union{NTuple{N, Int}, NTuple{N, Symbol}}
        ) where N

Fock coupling.

Type alias for `Coupling{V, I<:ID{FID}, C<:Subscripts}`.
"""
const FockCoupling{V, I<:ID{FID}, C<:Subscripts} = Coupling{V, I, C}
function FockCoupling(value::Number,
        orbitals::Subscript{<:NTuple{N, Union{Int, Symbol}}},
        spins::Subscript{<:NTuple{N, Union{Int, Symbol}}},
        nambus::Union{NTuple{N, Int}, NTuple{N, Symbol}}
        ) where N
    return Coupling(value, ID(FID{wildcard}, orbitals.pattern, spins.pattern, nambus), Subscripts((orbital=orbitals, spin=spins)))
end
function Base.show(io::IO, fc::FockCoupling)
    wildcards = ntuple(i->wildcard, Val(rank(fc)))
    @printf io "FockCoupling{%s}(value=%s" rank(fc) decimaltostr(fc.value)
    fc.cid.orbitals≠wildcards && @printf io ", orbitals=%s" repr(fc.subscripts, 1:length(fc.subscripts), :orbital)
    fc.cid.spins≠wildcards && @printf io ", spins=%s" repr(fc.subscripts, 1:length(fc.subscripts), :spin)
    fc.cid.nambus≠wildcards && @printf io ", nambus=[%s]" join(fc.cid.nambus, " ")
    @printf io ")"
end
function Base.repr(fc::FockCoupling)
    contents = String[]
    wildcards = ntuple(i->wildcard, Val(rank(fc)))
    fc.cid.orbitals≠wildcards && push!(contents, @sprintf "ob%s" repr(fc.subscripts, 1:length(fc.subscripts), :orbital))
    fc.cid.spins≠wildcards && push!(contents, @sprintf "sp%s" repr(fc.subscripts, 1:length(fc.subscripts), :spin))
    fc.cid.nambus≠wildcards && push!(contents, @sprintf "ph[%s]" join(fc.cid.nambus, " "))
    result = decimaltostr(fc.value)
    length(contents)==0 && (result = @sprintf "%s {%s}" result rank(fc))
    length(contents)>0 && (result = @sprintf "%s %s" result join(contents, " ⊗ "))
    return result
end

"""
    FockCoupling{N}(value::Number=1;
        orbitals::Union{NTuple{N, Int}, Subscript}=Subscript(N),
        spins::Union{NTuple{N, Int}, Subscript}=Subscript(N),
        nambus::Union{NTuple{N, Int}, NTuple{N, Symbol}}=ntuple(i->wildcard, Val(N))
        ) where N

Construct a Fock coupling.
"""
function FockCoupling{N}(value::Number=1;
        orbitals::Union{NTuple{N, Int}, Subscript}=Subscript(N),
        spins::Union{NTuple{N, Int}, Subscript}=Subscript(N),
        nambus::Union{NTuple{N, Int}, NTuple{N, Symbol}}=ntuple(i->wildcard, Val(N))
        ) where N
    isa(orbitals, Subscript) || (orbitals = Subscript(orbitals))
    isa(spins, Subscript) || (spins = Subscript(spins))
    return FockCoupling(value, orbitals, spins, nambus)
end

"""
    ⊗(fc₁::FockCoupling, fc₂::FockCoupling) -> FockCoupling

Get the direct product between two Fock couplings.
"""
function ⊗(fc₁::FockCoupling, fc₂::FockCoupling)
    @assert rank(fc₁)==rank(fc₂) "⊗ error: mismatched rank."
    @assert length(fc₁.subscripts)==length(fc₂.subscripts)==1 "⊗ error: not supported."
    wildcards = ntuple(i->wildcard, Val(rank(fc₁)))
    orbitals = fockcouplingchoicerule(first(fc₁.subscripts).orbital, first(fc₂.subscripts).orbital, wildcards, :orbitals)
    spins = fockcouplingchoicerule(first(fc₁.subscripts).spin, first(fc₂.subscripts).spin, wildcards, :spins)
    nambus = fockcouplingchoicerule(fc₁.cid.nambus, fc₂.cid.nambus, wildcards, :nambus)
    return FockCoupling(fc₁.value*fc₂.value, orbitals, spins, nambus)
end
function fockcouplingchoicerule(v₁, v₂, wildcards::Tuple{Vararg{Symbol}}, field::Symbol)
    t₁, t₂ = Tuple(v₁)==wildcards, Tuple(v₂)==wildcards
    @assert t₁||t₂ "fockcouplingchoicerule error: mismatched $field ($v₁ and $v₂)."
    return t₁ ? v₂ : v₁
end

"""
    ⋅(fc₁::FockCoupling, fc₂::FockCoupling) -> FockCoupling

Get the dot product of two rank-2 Fock couplings.

A rank-2 FockCoupling can be considered as a matrix acting on the orbital, spin and nambu spaces.
Therefore, the dot product here is defined as the multiplication between such matrices.
"""
function ⋅(fc₁::FockCoupling, fc₂::FockCoupling)
    @assert rank(fc₁)==rank(fc₂)==2 "⋅ error: input fockcouplings must be rank-2."
    value = fc₁.value * fc₂.value
    attrpairs = Pair{Symbol, NTuple{2, Int}}[]
    for attrname in (:orbitals, :spins, :nambus)
        v₁ = getproperty(fc₁.cid, attrname)
        v₂ = getproperty(fc₂.cid, attrname)
        if isa(v₁, NTuple{2, Int}) && isa(v₂, NTuple{2, Int})
            value = value * delta(v₁[2], v₂[1])
            push!(attrpairs, attrname=>(v₁[1], v₂[2]))
        else
            @assert v₁==(wildcard, wildcard) && v₂==(wildcard, wildcard) "⋅ error: mismatched $attrname ($v₁ and $v₂)."
        end
    end
    return FockCoupling{2}(value; attrpairs...)
end

"""
    σ⁰"sp" -> Couplings
    σ⁰"ob" -> Couplings
    σ⁰"ph" -> Couplings

The Pauli matrix σ⁰, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σ⁰_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σ⁰_str error: wrong input mode."
    mode=="sp" && return FockCoupling{2}(1, spins=(1, 1)) + FockCoupling{2}(1, spins=(2, 2))
    mode=="ob" && return FockCoupling{2}(1, orbitals=(1, 1)) + FockCoupling{2}(1, orbitals=(2, 2))
    return FockCoupling{2}(1, nambus=(annihilation, creation)) + FockCoupling{2}(1, nambus=(creation, annihilation))
end

"""
    σˣ"sp" -> Couplings
    σˣ"ob" -> Couplings
    σˣ"ph" -> Couplings

The Pauli matrix σˣ, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σˣ_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σˣ_str error: wrong input mode."
    mode=="sp" && return FockCoupling{2}(1, spins=(1, 2)) + FockCoupling{2}(1, spins=(2, 1))
    mode=="ob" && return FockCoupling{2}(1, orbitals=(1, 2)) + FockCoupling{2}(1, orbitals=(2, 1))
    return FockCoupling{2}(1, nambus=(annihilation, annihilation)) + FockCoupling{2}(1, nambus=(creation, creation))
end

"""
    σʸ"sp" -> Couplings
    σʸ"ob" -> Couplings
    σʸ"ph" -> Couplings

The Pauli matrix σʸ, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σʸ_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σʸ_str error: wrong input mode."
    mode=="sp" && return FockCoupling{2}(1im, spins=(1, 2)) + FockCoupling{2}(-1im, spins=(2, 1))
    mode=="ob" && return FockCoupling{2}(1im, orbitals=(1, 2)) + FockCoupling{2}(-1im, orbitals=(2, 1))
    return FockCoupling{2}(1im, nambus=(annihilation, annihilation)) + FockCoupling{2}(-1im, nambus=(creation, creation))
end

"""
    σᶻ"sp" -> Couplings
    σᶻ"ob" -> Couplings
    σᶻ"ph" -> Couplings

The Pauli matrix σᶻ, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σᶻ_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σᶻ_str error: wrong input mode."
    mode=="sp" && return FockCoupling{2}(-1, spins=(1, 1)) + FockCoupling{2}(1, spins=(2, 2))
    mode=="ob" && return FockCoupling{2}(-1, orbitals=(1, 1)) + FockCoupling{2}(1, orbitals=(2, 2))
    return FockCoupling{2}(-1, nambus=(annihilation, creation)) + FockCoupling{2}(1, nambus=(creation, annihilation))
end

"""
    σ⁺"sp" -> Couplings
    σ⁺"ob" -> Couplings
    σ⁺"ph" -> Couplings

The Pauli matrix σ⁺, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σ⁺_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σ⁺_str error: wrong input mode."
    mode=="sp" && return Couplings(FockCoupling{2}(1, spins=(2, 1)))
    mode=="ob" && return Couplings(FockCoupling{2}(1, orbitals=(2, 1)))
    return Couplings(FockCoupling{2}(1, nambus=(creation, creation)))
end

"""
    σ⁻"sp" -> Couplings
    σ⁻"ob" -> Couplings
    σ⁻"ph" -> Couplings

The Pauli matrix σ⁻, which can act on the space of spins("sp"), orbitals("ob") or particle-holes("ph").
"""
macro σ⁻_str(mode::String)
    @assert mode∈("sp", "ob", "ph") "@σ⁻_str error: wrong input mode."
    mode=="sp" && return Couplings(FockCoupling{2}(1, spins=(1, 2)))
    mode=="ob" && return Couplings(FockCoupling{2}(1, orbitals=(1, 2)))
    return Couplings(FockCoupling{2}(1, nambus=(annihilation, annihilation)))
end

"""
    fc"..." -> FockCoupling

Construct a FockCoupling from a literal string.
"""
macro fc_str(str::String)
    pos=findfirst(x->x∈('{', 's', 'o', 'p'), str)
    @assert 0<pos<length(str) "@fc_str error: wrong input pattern."
    coeff = eval(Meta.parse(str[1:pos-1]))
    ps = strip(replace(str[pos:end], " ⊗ "=>"⊗"))
    if ps[1]=='{' && ps[end]=='}'
        N = parse(Int, ps[2:end-1])
        return FockCoupling{N}(coeff)
    else
        attrpairs = []
        N = nothing
        for component in split(ps, '⊗')
            attrpair, rank = fccomponent(component)
            isnothing(N) && (N = rank)
            @assert N==rank "@fc_str error: mismatched ranks."
            push!(attrpairs, attrpair)
        end
        @assert length(attrpairs)>0 "@fc_str error: wrong input pattern."
        return Expr(:call, :(FockCoupling{$N}), Expr(:parameters, attrpairs...), coeff)
    end
end
@inline σᵅname(mode) = mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="ph" ? :nambus : error("σᵅname error: wrong input mode.")
function fccomponent(str::AbstractString)
    attrname = σᵅname(str[1:2])
    expr = Meta.parse(str[3:end])
    if attrname==:nambus
        @assert expr.head∈(:hcat, :vect) "fccomponent error: wrong input pattern for nambus."
        N = length(expr.args)
        attrvalue = Tuple(expr.args)
    else
        @assert expr.head∈(:call, :hcat, :vect) "fccomponent error: wrong input pattern for orbitals and spins."
        N = expr.head∈(:hcat, :vect) ? length(expr.args) : length(expr.args[1].args)
        attrvalue = subscriptexpr(expr)
    end
    return Expr(:kw, attrname, attrvalue), N
end

"""
    Onsite(id::Symbol, value;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false,
        )

Onsite term.

Type alias for `Term{:Onsite, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const Onsite{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:Onsite, id, V, Int, C, A, M}
@inline function Onsite(id::Symbol, value;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    Term{:Onsite}(id, value, 0, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:Onsite}) = :st
@inline ishermitian(::Type{<:Onsite}) = nothing

"""
    Hopping(id::Symbol, value, bondkind::Int=1;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false,
        )

Hopping term.

Type alias for `Term{:Hopping, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const Hopping{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:Hopping, id, V, Int, C, A, M}
@inline function Hopping(id::Symbol, value, bondkind::Int=1;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    @assert bondkind≠0 "Hopping error: input bondkind (neighbor) cannot be 0. Use `Onsite` instead."
    Term{:Hopping}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:Hopping}) = :hp
@inline ishermitian(::Type{<:Hopping}) = false
@inline couplingcenters(::FockCoupling, ::Bond, ::Val{:Hopping}) = (1, 2)

"""
    Pairing(id::Symbol, value, bondkind::Int=0;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false,
        )

Pairing term.

Type alias for `Term{:Pairing, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const Pairing{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:Pairing, id, V, Int, C, A, M}
@inline function Pairing(id::Symbol, value, bondkind::Int=0;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    Term{:Pairing}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:Pairing}) = :pr
@inline ishermitian(::Type{<:Pairing}) = false
@inline couplingcenters(::FockCoupling, ::Bond, ::Val{:Pairing}) = (1, 2)
@inline fidnambu(::Val{:Pairing}, ::Symbol, nnambu::Int, order::Int) = ((@assert nnambu==2 "fidnambu error: nnambu must be 2 for Pairing."); annihilation)
function expand!(operators::Operators, term::Pairing, bond::AbstractBond, hilbert::Hilbert; half::Bool=false, table::Union{Nothing, Table}=nothing)
    argtypes = Tuple{Operators, Term, AbstractBond, Hilbert}
    invoke(expand!, argtypes, operators, term, bond, hilbert; half=half, table=table)
    isa(bond, Bond) && invoke(expand!, argtypes, operators, term, reverse(bond), hilbert; half=half, table=table)
    return operators
end

"""
    Hubbard(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)

Hubbard term.

Type alias for `Term{:Hubbard, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const Hubbard{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:Hubbard, id, V, Int, C, A, M}
@inline function Hubbard(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)
    Term{:Hubbard}(id, value, 0, couplings=@couplings(fc"1 sp[2 2 1 1] ⊗ ph[2 1 2 1]"), amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:Hubbard}) = :hb
@inline ishermitian(::Type{<:Hubbard}) = true

"""
    InterOrbitalInterSpin(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)

Interorbital-interspin term.

Type alias for `Term{:InterOrbitalInterSpin, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const InterOrbitalInterSpin{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:InterOrbitalInterSpin, id, V, Int, C, A, M}
@inline function InterOrbitalInterSpin(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)
    Term{:InterOrbitalInterSpin}(id, value, 0;
        couplings=@couplings(fc"1 ob[α α β β](α < β) ⊗ sp[σ₁ σ₁ σ₂ σ₂](σ₁ ≠ σ₂) ⊗ ph[2 1 2 1]"),
        amplitude=amplitude,
        modulate=modulate
        )
end
@inline abbr(::Type{<:InterOrbitalInterSpin}) = :nons
@inline ishermitian(::Type{<:InterOrbitalInterSpin}) = true

"""
    InterOrbitalIntraSpin(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)

Interorbital-intraspin term.

Type alias for `Term{:InterOrbitalIntraSpin, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const InterOrbitalIntraSpin{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:InterOrbitalIntraSpin, id, V, Int, C, A, M}
@inline function InterOrbitalIntraSpin(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)
    Term{:InterOrbitalIntraSpin}(id, value, 0, couplings=@couplings(fc"1 ob[α α β β](α < β) ⊗ ph[2 1 2 1]"), amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:InterOrbitalIntraSpin}) = :noes
@inline ishermitian(::Type{<:InterOrbitalIntraSpin}) = true

"""
    SpinFlip(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)

Spin-flip term.

Type alias for `Term{:SpinFlip, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const SpinFlip{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:SpinFlip, id, V, Int, C, A, M}
@inline function SpinFlip(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)
    Term{:SpinFlip}(id, value, 0, couplings=@couplings(fc"1 ob[α β α β](α < β) ⊗ sp[2 1 1 2] ⊗ ph[2 2 1 1]"), amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:SpinFlip}) = :sf
@inline ishermitian(::Type{<:SpinFlip}) = false

"""
    PairHopping(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)

Pair-hopping term.

Type alias for `Term{:PairHopping, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const PairHopping{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:PairHopping, id, V, Int, C, A, M}
@inline function PairHopping(id::Symbol, value; amplitude::Union{Function, Nothing}=nothing, modulate::Union{Function, Bool}=false)
    Term{:PairHopping}(id, value, 0, couplings=@couplings(fc"1 ob[α α β β](α < β) ⊗ sp[2 1 1 2] ⊗ ph[2 2 1 1]"), amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:PairHopping}) = :ph
@inline ishermitian(::Type{<:PairHopping}) = false

"""
    Coulomb(id::Symbol, value, bondkind::Int=1;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()*FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )

Coulomb term.

Type alias for `Term{:Coulomb, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const Coulomb{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:Coulomb, id, V, Int, C, A, M}
@inline function Coulomb(id::Symbol, value, bondkind::Int=1;
        couplings::Union{Function, Couplings}=@couplings(FockCoupling{2}()*FockCoupling{2}()),
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    @assert bondkind≠0 "Coulomb error: input bondkind cannot be 0. Use `Hubbard/InterOrbitalInterSpin/InterOrbitalIntraSpin/SpinFlip/PairHopping` instead."
    Term{:Coulomb}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:Coulomb}) = :cl
@inline ishermitian(::Type{<:Coulomb}) = nothing
@inline termfactor(id::ID{OID, 4}, ::Val{:Coulomb}) = id[2]'==id[1] && id[4]'==id[3] ? 2 : 1
@inline couplingcenters(::FockCoupling, ::Bond, ::Val{:Coulomb}) = (1, 1, 2, 2)

"""
    FockTerm

Type alias for `Union{Onsite, Hopping, Pairing, Hubbard, InterOrbitalInterSpin, InterOrbitalIntraSpin, SpinFlip, PairHopping, Coulomb}`.
"""
const FockTerm = Union{Onsite, Hopping, Pairing, Hubbard, InterOrbitalInterSpin, InterOrbitalIntraSpin, SpinFlip, PairHopping, Coulomb}

# SU(2) spin systems
const sidtagmap = Dict(1=>'x', 2=>'y', 3=>'z', 4=>'+', 5=>'-')
const sidseqmap = Dict(v=>k for (k, v) in sidtagmap)
const sidajointmap = Dict('x'=>'x', 'y'=>'y', 'z'=>'z', '+'=>'-', '-'=>'+')
const sidrepmap = Dict('x'=>'ˣ', 'y'=>'ʸ', 'z'=>'ᶻ', '+'=>'⁺', '-'=>'⁻', '0'=>'⁰')
const sidreprevmap = Dict(v=>k for (k, v) in sidrepmap)

"""
    SID{S, O<:Union{Int, Symbol}} <: SimpleIID

The spin id.
"""
struct SID{S, O<:Union{Int, Symbol}, T<:Union{Char, Symbol}} <: SimpleIID
    orbital::O
    tag::T
    function SID{S}(orbital::Union{Int, Symbol}, tag::Union{Char, Symbol}) where S
        @assert isa(S, Rational{Int}) && S.den==2 || isa(S, Integer) || S==wildcard "SID error: not supported spin($S)."
        isa(tag, Char) && @assert tag in ('x', 'y', 'z', '+', '-') "SID error: not supported tag($tag)."
        new{S, typeof(orbital), typeof(tag)}(orbital, tag)
    end
end
@inline Base.adjoint(sid::SID) = SID{totalspin(sid)}(sid.orbital, sidajointmap[sid.tag])
Base.show(io::IO, sid::SID) = @printf io "SID{%s}(%s)" totalspin(sid) join(map(repr, ntuple(i->getfield(sid, i), Val(fieldcount(typeof(sid))))), ", ")
@inline @generated function Base.replace(sid::SID; kwargs...)
    exprs = [:(get(kwargs, $name, getfield(sid, $name))) for name in QuoteNode.(fieldnames(sid))]
    return :(rawtype(typeof(sid)){totalspin(sid)}($(exprs...)))
end
@inline totalspin(sid::SID) = totalspin(typeof(sid))
@inline totalspin(::Type{<:SID{S}}) where S = S
@inline statistics(::Type{<:SID}) = :b
@inline SID(iid::SID, ::CompositeInternal) = iid

"""
    SID{S}(tag::Union{Char, Symbol}; orbital::Union{Int, Symbol}=1) where S

Create a spin id.
"""
@inline SID{S}(tag::Union{Char, Symbol}; orbital::Union{Int, Symbol}=1) where S = SID{S}(orbital, tag)

"""
    Matrix(sid::SID{S, <:Union{Int, Symbol}, Char}, dtype::Type{<:Number}=Complex{Float}) where S -> Matrix{dtype}

Get the matrix representation of a sid.
"""
function Base.Matrix(sid::SID{S, <:Union{Int, Symbol}, Char}, dtype::Type{<:Number}=Complex{Float}) where S
    N = Int(2*S+1)
    result = zeros(dtype, (N, N))
    spin = convert(dtype, S)
    for i = 1:N, j = 1:N
        row, col = N+1-i, N+1-j
        m, n = spin+1-i, spin+1-j
        result[row, col] = (sid.tag == 'x') ? (delta(i+1, j)+delta(i, j+1))*sqrt(spin*(spin+1)-m*n)/2 :
            (sid.tag == 'y') ? (delta(i+1, j)-delta(i, j+1))*sqrt(spin*(spin+1)-m*n)/2im :
            (sid.tag == 'z') ? delta(i, j)*m :
            (sid.tag == '+') ? delta(i+1, j)*sqrt(spin*(spin+1)-m*n) :
            delta(i, j+1)*sqrt(spin*(spin+1)-m*n)
    end
    return result
end

"""
    Spin{S} <: SimpleInternal{SID{S, Int, Char}}

The spin internal degrees of freedom.
"""
struct Spin{S} <: SimpleInternal{SID{S, Int, Char}}
    norbital::Int
    function Spin{S}(norbital::Int) where S
        @assert isa(S, Rational{Int}) && S.den==2 || isa(S, Integer) "Spin error: not supported spin($S)."
        new{S}(norbital)
    end
end
@inline Base.eltype(::Type{Spin}) = (SID{S, Int, Char} where S)
@inline shape(sp::Spin) = (1:sp.norbital, 1:length(sidtagmap))
@inline Base.CartesianIndex(sid::SID, ::Spin) = CartesianIndex(sid.orbital, sidseqmap[sid.tag])
@inline SID(index::CartesianIndex{2}, sp::Spin) = SID{totalspin(sp)}(index[1], sidtagmap[index[2]])
Base.summary(io::IO, spin::Spin) = @printf io "%s-element Spin{%s}" length(spin) totalspin(spin)
@inline totalspin(spin::Spin) = totalspin(typeof(spin))
@inline totalspin(::Type{<:Spin{S}}) where S = S
function Base.show(io::IO, spin::Spin)
    @printf io "%s{%s}(%s)" spin|>typeof|>nameof totalspin(spin) join(("$name=$(getfield(spin, name))" for name in spin|>typeof|>fieldnames), ", ")
end
@inline Base.match(::Type{<:SID{wildcard}}, ::Type{<:Spin{S}}) where S = true
@inline Base.match(::Type{<:SID{S}}, ::Type{<:Spin{S}}) where S = true
@inline Base.match(::Type{<:SID{S₁}}, ::Type{<:Spin{S₂}}) where {S₁, S₂} = false
@inline sidrange(::Symbol) = 1:5
@inline sidrange(::Symbol, n::Int) = 1:n
@inline sidrange(tag::Char) = (pos=sidseqmap[tag]; pos:pos)
@inline sidrange(v::Int, n::Int) = ((@assert 0<v<n+1 "shape error: orbital out of range."); v:v)
@inline shape(iidspace::IIDSpace{<:SID, <:Spin}) = (sidrange(iidspace.iid.orbital, iidspace.internal.norbital), sidrange(iidspace.iid.tag))

"""
    Spin{S}(; norbital::Int=1) where S

Construct a spin degrees of freedom.
"""
@inline Spin{S}(; norbital::Int=1) where S = Spin{S}(norbital)

"""
    script(::Val{:site}, index::Index{<:AbstractPID, <:SID}; kwargs...) -> Int
    script(::Val{:orbital}, index::Index{<:AbstractPID, <:SID}; kwargs...) -> Int
    script(::Val{:tag}, index::Index{<:AbstractPID, <:SID}; kwargs...) -> Char

Get the required script of a spin oid.
"""
@inline script(::Val{:site}, index::Index{<:AbstractPID, <:SID}; kwargs...) = index.pid.site
@inline script(::Val{:orbital}, index::Index{<:AbstractPID, <:SID}; kwargs...) = index.iid.orbital
@inline script(::Val{:tag}, index::Index{<:AbstractPID, <:SID}; kwargs...) = index.iid.tag

"""
    slatex

The default LaTeX format for a spin oid.
"""
const slatex = LaTeX{(:tag,), (:site,)}('S')
@inline latexname(::Type{<:Index{<:AbstractPID, <:SID}}) = Symbol("Index{AbstractPID, SID}")
@inline latexname(::Type{<:CompositeOID{<:Index{<:AbstractPID, <:SID}}}) = Symbol("CompositeOID{Index{AbstractPID, SID}}")
latexformat(Index{<:AbstractPID, <:SID}, slatex)
latexformat(CompositeOID{<:Index{<:AbstractPID, <:SID}}, slatex)

"""
    permute(id₁::OID{<:Index{<:AbstractPID, SID{S, Int, Char}}}, id₂::OID{<:Index{<:AbstractPID, SID{S, Int, Char}}}) where S -> Tuple{Vararg{Operator}}

Permute two spin oids and get the result.
"""
function permute(id₁::OID{<:Index{<:AbstractPID, SID{S, Int, Char}}}, id₂::OID{<:Index{<:AbstractPID, SID{S, Int, Char}}}) where S
    @assert id₁.index≠id₂.index || id₁.rcoord≠id₂.rcoord || id₁.icoord≠id₂.icoord "permute error: permuted ids should not be equal to each other."
    if id₁.index.pid==id₂.index.pid && id₁.index.iid.orbital==id₂.index.iid.orbital && id₁.rcoord==id₂.rcoord && id₁.icoord==id₂.icoord
        if id₁.index.iid.tag == 'x'
            id₂.index.iid.tag=='y' && return (Operator(+1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='z' && return (Operator(-1im, ID(permutesoid(id₁, 'y'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='+' && return (Operator(-1, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='-' && return (Operator(+1, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
        elseif id₁.index.iid.tag == 'y'
            id₂.index.iid.tag=='x' && return (Operator(-1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='z' && return (Operator(+1im, ID(permutesoid(id₁, 'x'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='+' && return (Operator(-1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='-' && return (Operator(-1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
        elseif id₁.index.iid.tag == 'z'
            id₂.index.iid.tag=='x' && return (Operator(+1im, ID(permutesoid(id₁, 'y'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='y' && return (Operator(-1im, ID(permutesoid(id₁, 'x'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='+' && return (Operator(+1, ID(id₂)), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='-' && return (Operator(-1, ID(id₂)), Operator(1, ID(id₂, id₁)))
        elseif id₁.index.iid.tag == '+'
            id₂.index.iid.tag=='x' && return (Operator(+1, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='y' && return (Operator(+1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='z' && return (Operator(-1, ID(id₁)), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='-' && return (Operator(+2, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
        elseif id₁.index.iid.tag == '-'
            id₂.index.iid.tag=='x' && return (Operator(-1, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='y' && return (Operator(1im, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='z' && return (Operator(+1, ID(id₁)), Operator(1, ID(id₂, id₁)))
            id₂.index.iid.tag=='+' && return (Operator(-2, ID(permutesoid(id₁, 'z'))), Operator(1, ID(id₂, id₁)))
        end
    else
        return (Operator(1, ID(id₂, id₁)),)
    end
end
@inline permutesoid(id::OID{<:Index{<:AbstractPID, <:SID}}, tag::Char) = replace(id, index=replace(id.index, iid=replace(id.index.iid, tag=tag)))

"""
    SpinCoupling(value::Number, tags::NTuple{N, Char}, orbitals::Subscript{<:NTuple{N, Union{Int, Symbol}}}) where N

Spin coupling.

Type alias for `Coupling{V, I<:ID{SID}, C<:Subscripts}`.
"""
const SpinCoupling{V, I<:ID{SID}, C<:Subscripts} = Coupling{V, I, C}
@inline function SpinCoupling(value::Number, tags::NTuple{N, Char}, orbitals::Subscript{<:NTuple{N, Union{Int, Symbol}}}) where N
    return Coupling(value, ID(SID{wildcard}, orbitals.pattern, tags), Subscripts((orbital=orbitals,)))
end
function Base.show(io::IO, sc::SpinCoupling)
    @printf io "SpinCoupling(value=%s" decimaltostr(sc.value)
    @printf io ", tags=%s" join(NTuple{rank(sc), String}("S"*sidrepmap[tag] for tag in sc.cid.tags), "")
    sc.cid.orbitals≠ntuple(i->wildcard, Val(rank(sc))) && @printf io ", orbitals=%s" repr(sc.subscripts, 1:length(sc.subscripts), :orbital)
    @printf io ")"
end
function Base.repr(sc::SpinCoupling)
    result = [@sprintf "%s %s" decimaltostr(sc.value) join(NTuple{rank(sc), String}("S"*sidrepmap[tag] for tag in sc.cid.tags), "")]
    sc.cid.orbitals≠ntuple(i->wildcard, Val(rank(sc))) && push!(result, @sprintf "ob%s" repr(sc.subscripts, 1:length(sc.subscripts), :orbital))
    return join(result, " ")
end

"""
    SpinCoupling(value::Number, tags::NTuple{N, Char}; orbitals::Union{NTuple{N, Int}, Subscript}=Subscript(N)) where N

Construct a spin coupling.
"""
function SpinCoupling(value::Number, tags::NTuple{N, Char}; orbitals::Union{NTuple{N, Int}, Subscript}=Subscript(N)) where N
    isa(orbitals, Subscript) || (orbitals = Subscript(orbitals))
    return SpinCoupling(value, tags, orbitals)
end

"""
    sc"..." -> SpinCoupling

Construct a SpinCoupling from a literal string.
"""
macro sc_str(str::String)
    fpos = findfirst(r"S[⁺⁻ˣʸᶻ⁰]", str).start
    lpos = 1 + lastindex(str) - findfirst(r"[⁺⁻ˣʸᶻ⁰]S", str|>reverse).start
    coeff = eval(Meta.parse(str[firstindex(str):prevind(str, fpos)]))
    tags = Tuple(sidreprevmap[tag] for tag in replace(str[thisind(str, fpos):thisind(str, lpos)], "S"=>""))
    orbitals = scorbitals(strip(str[nextind(str, lpos):end]), length(tags))
    return Expr(:call, :SpinCoupling, orbitals, coeff, tags)
end
function scorbitals(str::AbstractString, n::Int)
    length(str)==0 && return Expr(:kw, :orbitals, Subscript(n))
    @assert str[1:2]=="ob" "scorbitals error: wrong input pattern."
    expr = Meta.parse(str[3:end])
    @assert expr.head∈(:call, :hcat, :vect) "scorbitals error: wrong input pattern for orbitals."
    attrvalue = subscriptexpr(expr)
    return Expr(:kw, :orbitals, attrvalue)
end

"""
    heisenberg"ob[o₁ o₂]" -> Couplings
    heisenberg"xyz ob[o₁ o₂]" -> Couplings
    heisenberg"+-z ob[o₁ o₂]" -> Couplings

The Heisenberg couplings.
"""
macro heisenberg_str(str::String)
    slice = findfirst(r"xyz|\+\-z", str)
    mode = isnothing(slice) ? "+-z" : str[slice]
    isnothing(slice) || (str = strip(str[nextind(str, slice.stop):end]))
    @assert mode=="+-z" || mode=="xyz" "@heisenberg_str error: not supported mode($mode)."
    orbitals = scorbitals(str, 2)
    if mode == "+-z"
        return Expr(:call, :heisenbergpmz, orbitals)
    else
        return Expr(:call, :heisenbergxyz, orbitals)
    end
end
function heisenbergxyz(;orbitals=Subscript(2))
    return Couplings(
        SpinCoupling(1, ('x', 'x'), orbitals=orbitals),
        SpinCoupling(1, ('y', 'y'), orbitals=orbitals),
        SpinCoupling(1, ('z', 'z'), orbitals=orbitals)
    )
end
function heisenbergpmz(;orbitals=Subscript(2))
    return Couplings(
        SpinCoupling(1//2, ('+', '-'), orbitals=orbitals),
        SpinCoupling(1//2, ('-', '+'), orbitals=orbitals),
        SpinCoupling(1//1, ('z', 'z'), orbitals=orbitals)
    )
end

"""
    ising"x ob[o₁ o₂]" -> Couplings
    ising"y ob[o₁ o₂])" -> Couplings
    ising"z ob[o₁ o₂]" -> Couplings

The Ising couplings.
"""
macro ising_str(str::String)
    @assert str[1] ∈ ('x', 'y', 'z') "@ising_str error: wrong input pattern."
    length(str)>1 && @assert str[2]==' ' "@ising_str error: wrong input pattern."
    orbitals = scorbitals(str[3:end], 2)
    return :(Couplings(SpinCoupling(1, ($str[1], $str[1]), $orbitals)))
end

"""
    gamma"x ob[o₁ o₂]" -> Couplings
    gamma"y ob[o₁ o₂]" -> Couplings
    gamma"z ob[o₁ o₂]" -> Couplings

The Gamma couplings.
"""
macro gamma_str(str::String)
    @assert str[1] in ('x', 'y', 'z') "@gamma_str error: wrong input pattern."
    length(str)>1 && @assert str[2]==' ' "@gamma_str error: wrong input pattern."
    t₁, t₂ = str[1]=='x' ? ('y', 'z') : str[1]=='y' ? ('z', 'x') : ('x', 'y')
    orbitals = scorbitals(str[3:end], 2)
    return Expr(:call, :gamma, orbitals, t₁, t₂)
end
function gamma(t₁::Char, t₂::Char; orbitals=Subscript(2))
    return Couplings(
        SpinCoupling(1, (t₁, t₂), orbitals=orbitals),
        SpinCoupling(1, (t₂, t₁), orbitals=orbitals)
        )
end

"""
    dm"x ob[o₁ o₂]" -> Couplings
    dm"y ob[o₁ o₂]" -> Couplings
    dm"z ob[o₁ o₂]" -> Couplings

The DM couplings.
"""
macro dm_str(str::String)
    @assert str[1] in ('x', 'y', 'z') "@dm_str error: wrong input pattern."
    length(str)>1 && @assert str[2]==' ' "@dm_str error: wrong input pattern."
    orbitals = scorbitals(str[3:end], 2)
    t₁, t₂ = str[1]=='x' ? ('y', 'z') : str[1]=='y' ? ('z', 'x') : ('x', 'y')
    return Expr(:call, :dm, orbitals, t₁, t₂)
end
function dm(t₁::Char, t₂::Char; orbitals=Subscript(2))
    return Couplings(
        SpinCoupling(+1, (t₁, t₂), orbitals=orbitals),
        SpinCoupling(-1, (t₂, t₁), orbitals=orbitals)
        )
end

"""
    sˣ"ob[o]" -> Couplings
    sʸ"ob[o]" -> Couplings
    sᶻ"ob[o]" -> Couplings

The single Sˣ/Sʸ/Sᶻ coupling.
"""
macro sˣ_str(str::String) Expr(:call, :Couplings, Expr(:call, :SpinCoupling, scorbitals(str, 1), 1, ('x',))) end
macro sʸ_str(str::String) Expr(:call, :Couplings, Expr(:call, :SpinCoupling, scorbitals(str, 1), 1, ('y',))) end
macro sᶻ_str(str::String) Expr(:call, :Couplings, Expr(:call, :SpinCoupling, scorbitals(str, 1), 1, ('z',))) end

"""
    SpinTerm(id::Symbol, value::Any, bondkind::Any, couplings::Union{Function, Coupling, Couplings},
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false,
        )

Spin term.

Type alias for `Term{:SpinTerm, id, V, B<:Any, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const SpinTerm{id, V, B<:Any, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:SpinTerm, id, V, B, C, A, M}
@inline function SpinTerm(id::Symbol, value::Any, bondkind::Any, couplings::Union{Function, Couplings},
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    Term{:SpinTerm}(id, value, bondkind, couplings=couplings, amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:SpinTerm}) = :sp
@inline ishermitian(::Type{<:SpinTerm}) = true
@inline function couplingcenters(sc::SpinCoupling, ::Bond, ::Val{:SpinTerm})
    @assert rank(sc)%2==0 "couplingcenters error: the rank of the input spin coupling should be even."
    return ntuple(i->2-i%2, Val(rank(sc)))
end

# Phononic systems
"""
    NID{D<:Union{Char, Symbol}} <: SimpleIID

The phonon id.
"""
struct NID{D<:Union{Char, Symbol}} <: SimpleIID
    tag::Char
    dir::D
    function NID(tag::Char, dir::Union{Char, Symbol}=wildcard)
        @assert tag∈('p', 'u') "NID error: wrong tag($tag)."
        isa(dir, Char) && @assert dir∈('x', 'y', 'z') "NID error: wrong direction($dir)."
        new{typeof(dir)}(tag, dir)
    end
end
@inline Base.adjoint(pnid::NID) = pnid
@inline statistics(::Type{<:NID}) = :b
@inline NID(iid::NID, ::CompositeInternal) = iid

"""
    Phonon <: SimpleInternal{NID{Char}}

The phonon internal degrees of freedom.
"""
struct Phonon <: SimpleInternal{NID{Char}}
    ndir::Int
    function Phonon(ndir::Integer)
        @assert ndir∈(1, 2, 3) "Phonon error: wrong number of directions."
        new(ndir)
    end
end
@inline shape(pn::Phonon) = (1:2, 1:pn.ndir)
@inline Base.CartesianIndex(pnid::NID{Char}, ::Phonon) = CartesianIndex(pnid.tag=='u' ? 1 : 2, Int(pnid.dir)-Int('x')+1)
@inline NID(index::CartesianIndex{2}, ::Phonon) = NID(index[1]==1 ? 'u' : 'p', Char(Int('x')+index[2]-1))
@inline shape(iidspace::IIDSpace{NID{Symbol}, Phonon}) = (iidspace.iid.tag=='u' ? (1:1) : (2:2), 1:iidspace.internal.ndir)
@inline function shape(iidspace::IIDSpace{NID{Char}, Phonon})
    dir = Int(iidspace.iid.dir)-Int('x')+1
    @assert 0<dir<iidspace.internal.ndir+1 "shape error: dir out of range."
    return (iidspace.iid.tag=='u' ? (1:1) : (2:2), dir:dir)
end

"""
    script(::Val{:BD}, index::Index{<:AbstractPID, <:NID}, l::LaTeX) -> Char
    script(::Val{:BD}, oid::CompositeOID{<:Index{<:AbstractPID, <:NID}}, l::LaTeX) -> Char
    script(::Val{:site}, index::Index{<:AbstractPID, <:NID}; kwargs...) -> Int
    script(::Val{:dir}, index::Index{<:AbstractPID, <:NID}; kwargs...) -> Char

Get the required script of a phonon oid.
"""
@inline script(::Val{:BD}, index::Index{<:AbstractPID, <:NID}, l::LaTeX) = l.body[index.iid.tag]
@inline script(::Val{:BD}, oid::CompositeOID{<:Index{<:AbstractPID, <:NID}}, l::LaTeX) = l.body[getcontent(oid, :index).iid.tag]
@inline script(::Val{:site}, index::Index{<:AbstractPID, <:NID}; kwargs...) = index.pid.site
@inline script(::Val{:dir}, index::Index{<:AbstractPID, <:NID}; kwargs...) = index.iid.dir

"""
    nlatex

The default LaTeX format for a phonon oid.
"""
const nlatex = LaTeX{(:dir,), (:site,)}(Dict('p'=>'p', 'u'=>'u'), "", "")
@inline latexname(::Type{<:Index{<:AbstractPID, <:NID}}) = Symbol("Index{AbstractPID, NID}")
@inline latexname(::Type{<:CompositeOID{<:Index{<:AbstractPID, <:NID}}}) = Symbol("CompositeOID{Index{AbstractPID, NID}}")
latexformat(Index{<:AbstractPID, <:NID}, nlatex)
latexformat(CompositeOID{<:Index{<:AbstractPID, <:NID}}, nlatex)

"""
    permute(id₁::OID{<:Index{<:AbstractPID, NID{Char}}}, id₂::OID{<:Index{<:AbstractPID, NID{Char}}}) -> Tuple{Vararg{Operator}}

Permute two phonon oids and get the result.
"""
function permute(id₁::OID{<:Index{<:AbstractPID, NID{Char}}}, id₂::OID{<:Index{<:AbstractPID, NID{Char}}})
    if id₁.index.iid.dir==id₂.index.iid.dir && id₁.index.iid.tag≠id₂.index.iid.tag && id₁.rcoord==id₂.rcoord && id₁.icoord==id₂.icoord
        if id₁.index.iid.tag=='u'
            return (Operator(1im), Operator(1, ID(id₂, id₁)))
        else
            return (Operator(-1im), Operator(1, ID(id₂, id₁)))
        end
    else
        return (Operator(1, ID(id₂, id₁)),)
    end
end

"""
    PhononCoupling(value::Number, tags::NTuple{N, Char}, dirs::Subscript{<:Union{NTuple{N, Char}, NTuple{N, Symbol}}}) where N

Phonon coupling.

Type alias for `Coupling{V<:Number, I<:ID{NID}, C<:Subscripts}`.
"""
const PhononCoupling{V<:Number, I<:ID{NID}, C<:Subscripts} = Coupling{V, I, C}
@inline function PhononCoupling(value::Number, tags::NTuple{N, Char}, dirs::Subscript{<:Union{NTuple{N, Char}, NTuple{N, Symbol}}}) where N
    return Coupling(value, ID(NID, tags, dirs.pattern), Subscripts((dir=dirs,)))
end
function Base.show(io::IO, pnc::PhononCoupling)
    @printf io "PhononCoupling(value=%s" decimaltostr(pnc.value)
    @printf io ", tags=[%s]" join(pnc.cid.tags, " ")
    pnc.cid.dirs≠ntuple(i->wildcard, Val(rank(pnc))) && @printf io ", dirs=%s" repr(pnc.subscripts, 1:length(pnc.subscripts), :dir)
    @printf io ")"
end
function Base.repr(pnc::PhononCoupling)
    result = [@sprintf "%s [%s]" decimaltostr(pnc.value) join(pnc.cid.tags, " ")]
    pnc.cid.dirs≠ntuple(i->wildcard, Val(rank(pnc))) && push!(result, @sprintf "dr%s" repr(pnc.subscripts, 1:length(pnc.subscripts), :dir))
    return join(result, " ")
end

"""
    PhononCoupling(value::Number, tags::NTuple{N, Char};
        dirs::Union{NTuple{N, Char}, NTuple{N, Symbol}, Subscript}=Subscript(N)
        ) where N

Construct a phonon coupling.
"""
function PhononCoupling(value::Number, tags::NTuple{N, Char};
        dirs::Union{NTuple{N, Char}, NTuple{N, Symbol}, Subscript}=Subscript(N)
        ) where N
    isa(dirs, Subscript) || (dirs = Subscript(dirs))
    return PhononCoupling(value, tags, dirs)
end

"""
    expand(pnc::PhononCoupling{<:Number, <:ID{NID{Symbol}}}, bond::Bond, hilbert::Hilbert, info::Val{:PhononPotential}) -> PPExpand

Expand the default phonon potential coupling on a given bond.
"""
function expand(pnc::PhononCoupling{<:Number, <:ID{NID{Symbol}}}, bond::Bond, hilbert::Hilbert, info::Val{:PhononPotential})
    R̂ = rcoord(bond)/norm(rcoord(bond))
    @assert pnc.cid.tags==('u', 'u') "expand error: wrong tags of phonon coupling."
    @assert isapprox(pnc.value, 1, atol=atol, rtol=rtol) "expand error: wrong coefficient of phonon coupling."
    pn₁, pn₂ = couplinginternals(pnc, bond, hilbert, info)
    @assert pn₁.ndir==pn₂.ndir==length(R̂) "expand error: mismatched number of directions."
    return PPExpand(R̂, (bond.epoint, bond.spoint))
end
struct PPExpand{N, P<:AbstractPID, D<:Number} <: VectorSpace{Operator{D, ID{OID{Index{P, NID{Char}}, SVector{N, D}}, 2}}}
    direction::SVector{N, D}
    points::NTuple{2, Point{N, P, D}}
end
@inline VectorSpaceStyle(::Type{<:PPExpand}) = VectorSpaceCartesian()
@inline shape(pnce::PPExpand) = (1:length(pnce.direction), 1:length(pnce.direction), 1:4)
function Operator(index::CartesianIndex{3}, pnce::PPExpand)
    dir₁ = Char(Int('x')+index[1]-1)
    dir₂ = Char(Int('x')+index[2]-1)
    coeff = index[3]∈(1, 4) ? 1 : -1
    pos₁, pos₂ = index[3]==1 ? (1, 1) : index[3]==2 ? (1, 2) : index[3]==3 ? (2, 1) : (2, 2)
    oid₁ = OID(Index(pnce.points[pos₁].pid, NID('u', dir₁)), pnce.points[pos₁].rcoord, pnce.points[pos₁].icoord)
    oid₂ = OID(Index(pnce.points[pos₂].pid, NID('u', dir₂)), pnce.points[pos₂].rcoord, pnce.points[pos₂].icoord)
    return Operator(pnce.direction[index[1]]*pnce.direction[index[2]]*coeff, ID(oid₁, oid₂))
end

"""
    kinetic"" -> Couplings

The kinetic energy part of phonon couplings.
"""
macro kinetic_str(::String) Couplings(PhononCoupling(1, ('p', 'p'))) end

"""
    potential"" -> Couplings

The potential energy part of phonon couplings.
"""
macro potential_str(::String) Couplings(PhononCoupling(1, ('u', 'u'))) end

"""
    PhononKinetic(id::Symbol, value::Any;
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )

Kinetic energy part of phonons.

Type alias for `Term{:PhononKinetic, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`.
"""
const PhononKinetic{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:PhononKinetic, id, V, Int, C, A, M}
@inline function PhononKinetic(id::Symbol, value::Any;
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    Term{:PhononKinetic}(id, value, 0, couplings=kinetic"", amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:PhononKinetic}) = :pnk
@inline ishermitian(::Type{<:PhononKinetic}) = true

"""
    PhononPotential(id::Symbol, value::Any, bondkind::Int;
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )

Potential energy part of phonons.

Type alias for `Term{:PhononPotential, id, V, Int, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate}`
"""
const PhononPotential{id, V, C<:TermCouplings, A<:TermAmplitude, M<:TermModulate} = Term{:PhononPotential, id, V, Int, C, A, M}
@inline function PhononPotential(id::Symbol, value::Any, bondkind::Int;
        amplitude::Union{Function, Nothing}=nothing,
        modulate::Union{Function, Bool}=false
        )
    Term{:PhononPotential}(id, value, bondkind, couplings=potential"", amplitude=amplitude, modulate=modulate)
end
@inline abbr(::Type{<:PhononPotential}) = :pnp
@inline ishermitian(::Type{<:PhononPotential}) = true
@inline couplingcenters(::PhononCoupling, ::Bond, ::Val{:PhononPotential}) = (1, 2)

"""
    PhononTerm

Type alias for `Union{PhononKinetic, PhononPotential}`.
"""
const PhononTerm = Union{PhononKinetic, PhononPotential}

end # module


text:
#==
Utility functions
==#

"""
    amplitude(e)

Compute field amplitude in dB.
"""
amplitude(e) = 10log10(abs2(e))  # == 20log10(abs(E))

"""
    amplitudephase(e)

Compute field amplitude in dB and phase in radians and return as (`amplitude`, `phase`).
"""
function amplitudephase(e)
    a = amplitude(e)
    p = angle(e)

    return a, p
end

function amplitudephase(e::AbstractVector)
    a = similar(e, Float64)
    p = similar(a)
    @inbounds for i in eachindex(e)
        a[i] = amplitude(e[i])
        p[i] = angle(e[i])
    end
    return a, p
end

"""
    unwrap!(x)

Unwrap a phase vector `x` in radians in-place.
"""
unwrap!

function unwrap!(x::AbstractVector)
	v = first(x)
	@inbounds for k in eachindex(x)
		x[k] = v = v + rem2pi(x[k]-v, RoundNearest)
	end
	return x
end

unwrap!(x::Number) = x

"""
    pow23(x)

Efficiently compute ``x^(2/3)``.
"""
pow23

pow23(x::Real) = cbrt(x)^2
pow23(z::Complex) = exp(2/3*log(z))

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

module TableDistances

# basic imports
using Tables
using ScientificTypes
using Distances
using Statistics

# custom distances
using StringDistances
using CategoricalArrays
using CoDa

import Distances: pairwise, result_type

include("distances.jl")
include("normalizations.jl")
include("weights.jl")
include("api.jl")

export
  TableDistance,
  pairwise

end


text:
mutable struct PerftTree
    tot::Threads.Atomic{Int64}
    nodes::Array{Int64,1}
    div::Dict{AbstractString,Array{Int64,1}}
end

function perft(b::Board, max_depth::Int64)
    pt = PerftTree(Threads.Atomic{Int64}(0), zeros(max_depth), Dict{UInt64,Int64}())
    explore(pt, b, max_depth, 1)
    return pt
end

function explore(pt::PerftTree, b::Board, max_depth::Int64,
                 depth::Int64, root_move::AbstractString="")

    moves = getMoves(b, b.active)
    Threads.atomic_add!(pt.tot, length(moves.moves))
    pt.nodes[depth] += length(moves.moves)

    if root_move != ""
        pt.div[root_move][depth-1] += length(moves.moves)
    end

    if length(moves.moves) == 0 || depth == max_depth
        return pt
    end

    for m in moves.moves
        if m.type == :none; break end
        if depth == 1
            root_move = UINT2PGN[m.from]*UINT2PGN[m.to]
            if m.promotion != :none
                if m.promotion == :queen
                    root_move *= "q"
                elseif m.promotion == :rook
                    root_move *= "r"
                elseif m.promotion == :bishop
                    root_move *= "b"
                elseif m.promotion == :knight
                    root_move *= "n"
                end
            end
            pt.div[root_move] = zeros(max_depth)
        end
        explore(pt, makeMove(b, m), max_depth, depth+1, root_move)
    end
end


text:
using memo
using Documenter

DocMeta.setdocmeta!(memo, :DocTestSetup, :(using memo); recursive=true)

makedocs(;
    modules=[memo],
    authors="liuxu <liuxu89@mail.ustc.edu.cn> and contributors",
    repo="https://github.com/liuxu89/memo/blob/{commit}{path}#{line}",
    sitename="memo",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://liuxu89.github.io/memo",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/liuxu89/memo",
    devbranch="main",
)


text:
export IndexSet,
       hasindex,
       hasinds,
       hassameinds,
       findindex,
       findinds,
       indexposition,
       swaptags,
       swaptags!,
       swapprime,
       swapprime!,
       mapprime,
       mapprime!,
       commoninds,
       commonindex,
       uniqueinds,
       uniqueindex,
       dims,
       mindim,
       maxdim,
       push,
       permute

struct IndexSet{N}
  inds::MVector{N,Index}
  IndexSet{N}(inds::MVector{N,Index}) where {N} = new{N}(inds)
  IndexSet{0}(::MVector{0}) = new{0}(())
  IndexSet{N}(inds::NTuple{N,Index}) where {N} = new{N}(inds)
  IndexSet{0}() = new{0}(())
  IndexSet{0}(::Tuple{}) = new{0}(())
end
IndexSet(inds::MVector{N,Index}) where {N} = IndexSet{N}(inds)
IndexSet(inds::NTuple{N,Index}) where {N} = IndexSet{N}(inds)

inds(is::IndexSet) = is.inds

# Empty constructor
IndexSet() = IndexSet{0}()
IndexSet(::Tuple{}) = IndexSet()
IndexSet(::MVector{0}) = IndexSet()

# Construct of some size
IndexSet{N}() where {N} = IndexSet{N}(ntuple(_->Index(),Val(N)))
IndexSet(::Val{N}) where {N} = IndexSet{N}()

# Construct from various sets of indices
IndexSet{N}(inds::Vararg{Index,N}) where {N} = IndexSet{N}(NTuple{N,Index}(inds))
IndexSet(inds::Vararg{Index,N}) where {N} = IndexSet{N}(inds...)

IndexSet{N}(ivs::NTuple{N,IndexVal}) where {N} = IndexSet{N}(ntuple(i->ind(ivs[i]),Val(N)))
IndexSet(ivs::NTuple{N,IndexVal}) where {N} = IndexSet{N}(ivs)
IndexSet{N}(ivs::Vararg{IndexVal,N}) where {N} = IndexSet{N}(tuple(ivs...))
IndexSet(ivs::Vararg{IndexVal,N}) where {N} = IndexSet{N}(tuple(ivs...))

# Construct from various sets of IndexSets
IndexSet(inds::IndexSet) = inds
IndexSet(inds::IndexSet,i::Index) = IndexSet(inds...,i)
IndexSet(i::Index,inds::IndexSet) = IndexSet(i,inds...)
IndexSet(is1::IndexSet,is2::IndexSet) = IndexSet(is1...,is2...)

# This is used in type promotion in the Tensor contraction code
Base.promote_rule(::Type{<:IndexSet},::Type{Val{N}}) where {N} = IndexSet{N}

ValLength(::Type{IndexSet{N}}) where {N} = Val{N}
ValLength(::IndexSet{N}) where {N} = Val(N)

# TODO: make a version that accepts an arbitrary set of IndexSets
# as well as mixtures of seperate Indices and Tuples of Indices.
# Look at jointuples in the DenseTensor decomposition logic.
IndexSet(inds::NTuple{2,IndexSet}) = IndexSet(inds...)

# Convert to an Index if there is only one
Index(is::IndexSet) = length(is)==1 ? is[1] : error("Number of Index in IndexSet ≠ 1")

function Base.show(io::IO, is::IndexSet)
  for i in is.inds
    print(io,i)
    print(io," ")
  end
end

Base.getindex(is::IndexSet,n::Integer) = getindex(is.inds,n)
Base.setindex!(is::IndexSet,i::Index,n::Integer) = setindex!(is.inds,i,n)
Base.length(is::IndexSet{N}) where {N} = N
Base.length(::Type{IndexSet{N}}) where {N} = N
order(is::IndexSet) = length(is)
Base.copy(is::IndexSet) = IndexSet(copy(is.inds))
dims(is::IndexSet{N}) where {N} = ntuple(i->dim(is[i]),Val(N))
dim(is::IndexSet) = prod(dim.(is))
dim(is::IndexSet,pos::Integer) = dim(is[pos])

# Optimize this (right own function that extracts dimensions
# with a function)
Base.strides(is::IndexSet) = Base.size_to_strides(1, dims(is)...)
Base.stride(is::IndexSet,k::Integer) = strides(is)[k]

dag(is::IndexSet) = IndexSet(dag.(is.inds))

# Allow iteration
Base.iterate(is::IndexSet{N},state::Int=1) where {N} = state > N ? nothing : (is[state], state+1)

Base.eltype(is::Type{<:IndexSet}) = Index
Base.eltype(is::IndexSet) = eltype(typeof(is))

# Needed for findfirst (I think)
Base.keys(is::IndexSet{N}) where {N} = 1:N

StaticArrays.push(is::IndexSet{N},i::Index) where {N} = IndexSet{N+1}(push(is.inds,i))
StaticArrays.pushfirst(is::IndexSet{N},i::Index) where {N} = IndexSet{N+1}(pushfirst(is.inds,i))

unioninds(is1::IndexSet{N1},is2::IndexSet{N2}) where {N1,N2} = IndexSet{N1+N2}(is1...,is2...)

# This is to help with some generic programming in the Tensor
# code (it helps to construct an IndexSet(::NTuple{N,Index}) where the
# only known thing for dispatch is a concrete type such
# as IndexSet{4})
StaticArrays.similar_type(::Type{<:IndexSet},::Val{N}) where {N} = IndexSet{N}
StaticArrays.similar_type(::Type{<:IndexSet},::Type{Val{N}}) where {N} = IndexSet{N}

sim(is::IndexSet{N}) where {N} = IndexSet{N}(ntuple(i->sim(is[i]),Val(N)))

"""
mindim(is::IndexSet)

Get the minimum dimension of the indices in the index set.

Returns 1 if the IndexSet is empty.
"""
function mindim(is::IndexSet)
  length(is) == 0 && (return 1)
  md = dim(is[1])
  for n ∈ 2:length(is)
    md = min(md,dim(is[n]))
  end
  return md
end

"""
maxdim(is::IndexSet)

Get the maximum dimension of the indices in the index set.

Returns 1 if the IndexSet is empty.
"""
function maxdim(is::IndexSet)
  length(is) == 0 && (return 1)
  md = dim(is[1])
  for n ∈ 2:length(is)
    md = max(md,dim(is[n]))
  end
  return md
end

#
# Set operations
#

# inds has the index i
function hasindex(inds,i::Index)
  is = IndexSet(inds)
  for j ∈ is
    i==j && return true
  end
  return false
end

# Binds is subset of Ainds
function hasinds(Binds,Ainds)
  Ais = IndexSet(Ainds)
  for i ∈ Ais
    !hasindex(Binds,i) && return false
  end
  return true
end
hasinds(Binds,Ainds::Index...) = hasinds(Binds,IndexSet(Ainds...))

# Set equality (order independent)
function hassameinds(Ainds,Binds)
  Ais = IndexSet(Ainds)
  Bis = IndexSet(Binds)
  return hasinds(Ais,Bis) && length(Ais) == length(Bis)
end

"""
==(is1::IndexSet, is2::IndexSet)

IndexSet quality (order dependent)
"""
function ==(Ais::IndexSet,Bis::IndexSet)
  length(Ais) ≠ length(Bis) && return false
  for i ∈ 1:length(Ais)
    Ais[i] ≠ Bis[i] && return false
  end
  return true
end

# Helper function for uniqueinds
# Return true if the Index is not in any
# of the input sets of indices
function _is_unique_index(j::Index,inds::T) where {T<:Tuple}
  for I in inds
    hasindex(I,j) && return false
  end
  return true
end
# Version taking one ITensor or IndexSet
function _is_unique_index(j::Index,inds)
  hasindex(inds,j) && return false
  return true
end

"""
uniqueinds(Ais,Bis...)

Output the IndexSet with Indices in Ais but not in
the IndexSets Bis.
"""
function uniqueinds(Ainds,Binds)
  Ais = IndexSet(Ainds)
  Cis = IndexSet()
  for j ∈ Ais
    _is_unique_index(j,Binds) && (Cis = push(Cis,j))
  end
  return Cis
end

"""
uniqueindex(Ais,Bis)

Output the Index in Ais but not in the IndexSets Bis.
Otherwise, return a default constructed Index.

In the future, this may throw an error if more than
one Index is found.
"""
function uniqueindex(Ainds,Binds)
  Ais = IndexSet(Ainds)
  for j ∈ Ais
    _is_unique_index(j,Binds) && return j
  end
  return nothing
end
# This version can check for repeats, but is a bit
# slower because of IndexSet allocation
#uniqueindex(Ais,Bis) = Index(uniqueinds(Ais,Bis))

setdiff(Ais::IndexSet,Bis) = uniqueinds(Ais,Bis)

"""
commoninds(Ais,Bis)

Output the IndexSet in the intersection of Ais and Bis
"""
function commoninds(Ainds,Binds)
  Ais = IndexSet(Ainds)
  Cis = IndexSet()
  for i ∈ Ais
    hasindex(Binds,i) && (Cis = push(Cis,i))
  end
  return Cis
end

"""
commonindex(Ais,Bis)

Output the Index common to Ais and Bis.
If more than one Index is found, throw an error.
Otherwise, return a default constructed Index.
"""
function commonindex(Ainds,Binds)
  Ais = IndexSet(Ainds)
  for i ∈ Ais
    hasindex(Binds,i) && return i
  end
  return nothing
end
# This version checks if there are more than one indices
#commonindex(Ais,Bis) = Index(commoninds(Ais,Bis))

"""
findinds(inds,tags)

Output the IndexSet containing the subset of indices
of inds containing the tags in the input tagset.
"""
function findinds(inds,tags)
  is = IndexSet(inds)
  ts = TagSet(tags)
  found_inds = IndexSet()
  for i ∈ is
    if hastags(i,ts)
      found_inds = push(found_inds,i)
    end
  end
  return found_inds
end
"""
findindex(inds,tags)

Output the Index containing the tags in the input tagset.
If more than one Index is found, throw an error.
Otherwise, return a default constructed Index.
"""
function findindex(inds,tags)
  is = IndexSet(inds)
  ts = TagSet(tags)
  for i ∈ is
    if hastags(i,ts)
      return i
    end
  end
  # TODO: should this return `nothing` if no Index is found?
  return nothing
end
# This version checks if there are more than one indices
#findindex(inds, tags) = Index(findinds(inds,tags))

# TODO: Should this return `nothing` like `findfirst`?
# Should this just use `findfirst`?
function indexposition(is::IndexSet,
                       i::Index)
  for (n,j) in enumerate(is)
    if i==j
      return n
    end
  end
  return nothing
end

# From a tag set or index set, find the positions
# of the matching indices as a vector of integers
indexpositions(inds) = collect(1:length(inds))
indexpositions(inds, match::Nothing) = collect(1:length(inds))
#indexpositions(inds, match::Tuple{}) = collect(1:length(inds))
# Version for matching a tag set
function indexpositions(inds, match::T) where {T<:Union{AbstractString,
                                                        Tuple{<:AbstractString,<:Integer},
                                                        TagSet}}
  is = IndexSet(inds)
  tsmatch = TagSet(match)
  pos = Int[]
  for (j,I) ∈ enumerate(is)
    hastags(I,tsmatch) && push!(pos,j)
  end
  return pos
end
# Version for matching a collection of indices
function indexpositions(inds, match)
  is = IndexSet(inds)
  ismatch = IndexSet(match)
  pos = Int[]
  for (j,I) ∈ enumerate(is)
    hasindex(ismatch,I) && push!(pos,j)
  end
  return pos
end
# Version for matching a list of indices
indexpositions(inds, match_inds::Index...) = indexpositions(inds, IndexSet(match_inds...))

#
# Tagging functions
#

function prime!(is::IndexSet, plinc::Integer, match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = prime(is[jj],plinc)
  end
  return is
end
prime!(is::IndexSet,match=nothing) = prime!(is,1,match)
prime(is::IndexSet, vargs...) = prime!(copy(is), vargs...)
# For is' notation
adjoint(is::IndexSet) = prime(is)

function setprime!(is::IndexSet, plev::Integer, match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = setprime(is[jj],plev)
  end
  return is
end
setprime(is::IndexSet, vargs...) = setprime!(copy(is), vargs...)

noprime!(is::IndexSet, match = nothing) = setprime!(is, 0, match)
noprime(is::IndexSet, vargs...) = noprime!(copy(is), vargs...)

function swapprime!(is::IndexSet,
                    pl1::Int,
                    pl2::Int,
                    vargs...)
  pos = indexpositions(is,vargs...)
  for n in pos
    if plev(is[n])==pl1
      is[n] = setprime(is[n],pl2)
    elseif plev(is[n])==pl2
      is[n] = setprime(is[n],pl1)
    end
  end
  return is
end

swapprime(is::IndexSet,pl1::Int,pl2::Int,vargs...) = swapprime!(copy(is),pl1,pl2,vargs...)

function mapprime!(is::IndexSet,
                   plold::Integer,
                   plnew::Integer,
                   match = nothing)
  pos = indexpositions(is,match)
  for n in pos
    if plev(is[n])==plold
      is[n] = setprime(is[n],plnew)
    end
  end
  return is
end

function mapprime(is::IndexSet,
                  plold::Integer,
                  plnew::Integer,
                  match=nothing)
  return mapprime!(copy(is),plold,plnew,match)
end

function addtags!(is::IndexSet,
                  tags,
                  match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = addtags(is[jj],tags)
  end
  return is
end
addtags(is, vargs...) = addtags!(copy(is), vargs...)

function settags!(is::IndexSet,
                  ts,
                  match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = settags(is[jj],ts)
  end
  return is
end
settags(is, vargs...) = settags!(copy(is), vargs...)

function removetags!(is::IndexSet,
                     tags,
                     match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = removetags(is[jj],tags)
  end
  return is
end
removetags(is, vargs...) = removetags!(copy(is), vargs...)

function replacetags!(is::IndexSet,
                      tags_old, tags_new,
                      match = nothing)
  pos = indexpositions(is, match)
  for jj ∈ pos
    is[jj] = replacetags(is[jj],tags_old,tags_new)
  end
  return is
end
replacetags(is, vargs...) = replacetags!(copy(is), vargs...)

# TODO: write more efficient version in terms
# of indexpositions like swapprime!
function swaptags!(is::IndexSet,
                   tags1, tags2,
                   match = nothing)
  ts1 = TagSet(tags1)
  ts2 = TagSet(tags2)
  # TODO: add debug check that this "random" tag
  # doesn't clash with ts1 or ts2
  tstemp = TagSet("e43efds")
  plev(ts1) ≥ 0 && (tstemp = setprime(tstemp,431534))
  replacetags!(is, ts1, tstemp, match)
  replacetags!(is, ts2, ts1, match)
  replacetags!(is, tstemp, ts2, match)
  return is
end
swaptags(is, vargs...) = swaptags!(copy(is), vargs...)

#
# Helper functions for contracting ITensors
#

function compute_contraction_labels(Ai::IndexSet{N1},
                                    Bi::IndexSet{N2}) where {N1,N2}
  rA = length(Ai)
  rB = length(Bi)
  Aind = MVector{N1,Int}(ntuple(_->0,Val(N1)))
  Bind = MVector{N2,Int}(ntuple(_->0,Val(N2)))

  ncont = 0
  for i = 1:rA, j = 1:rB
    if Ai[i]==Bi[j]
      Aind[i] = Bind[j] = -(1+ncont)
      ncont += 1
    end
  end

  u = ncont
  for i = 1:rA
    if(Aind[i]==0) Aind[i] = (u+=1) end
  end
  for j = 1:rB
    if(Bind[j]==0) Bind[j] = (u+=1) end
  end

  return (NTuple{N1,Int}(Aind),NTuple{N2,Int}(Bind))
end

# Move this to tensor, since this logic is different
# for contracting different kinds of storage
# Also, generalize this to not just use IndexSet
#function contract_inds(Ais::IndexSet{N1},
#                       Alabels::Labels{N1},
#                       Bis::IndexSet{N2},
#                       Blabels::Labels{N2}) where {N1,N2}
#  ncont = 0
#  for i in Alabel
#    i < 0 && (ncont += 1)
#  end
#  NR = N1+N2-2*ncont
#  Clabels = Vector{Int}(undef,NR)
#  Cis = Vector{Index}(undef,NR)
#  u = 1
#  @inbounds for i ∈ 1:N1
#    if(Alabel[i] > 0)
#      Clabels[u] = Alabel[i];
#      Cis[u] = Ais[i];
#      u += 1
#    end
#  end
#  @inbounds for i ∈ 1:N2
#    if(Blabel[i] > 0)
#      Clabels[u] = Blabel[i];
#      Cis[u] = Bis[i];
#      u += 1
#    end
#  end
#  return IndexSet{NR}(Cis...),NTuple{NR,Int}(Clabels)
#end

# TODO: implement this in terms of a tuple,
# overload Base.strides and implement strides(inds,j)
# to get the jth stride
# TODO: should the IndexSet store the strides?
function compute_strides(inds::IndexSet)
  r = order(inds)
  stride = zeros(Int, r)
  s = 1
  for i = 1:r
    stride[i] = s
    s *= dim(inds[i])
  end
  return stride
end

#
# More general set functions
#

function permute(is::IndexSet{N},perm) where {N}
  indsp = ntuple(i->is[perm[i]], Val(N))
  return IndexSet(indsp)
end

# Permute some other type by perm
# (for example, tuple, MVector, etc.)
# as long as the constructor accepts a tuple
function permute(is::T,perm) where {T}
  indsp = ntuple(i->is[perm[i]], Val(length(is)))
  return T(indsp)
end

"""
getperm(col1,col2)

Get the permutation that takes collection 2 to collection 1,
such that col2[p].==col1
"""
function getperm(s1::Union{IndexSet{N},NTuple{N}}, s2::Union{IndexSet{N},NTuple{N}}) where {N}
  return ntuple(i->findfirst(==(s1[i]),s2),Val(N))
end

"""
getperm(col1,col2,col3)

Get the permutations that takes collections 2 and 3 to collection 1.
"""
function getperms(s::IndexSet{N},s1::IndexSet{N1},s2::IndexSet{N2}) where {N1,N2,N}
  N1+N2≠N && error("Size of partial sets don't match with total set")
  perm1 = ntuple(i->findfirst(==(s1[i]),s),Val(N1))
  perm2 = ntuple(i->findfirst(==(s2[i]),s),Val(N2))
  isperm((perm1...,perm2...)) || error("Combined permutations are $((perm1...,perm2...)), not a valid permutation")
  return perm1,perm2
end

"""
Determine if s1 and s2 have no overlapping elements.
"""
function isdisjoint(s1,s2)
  for i1 ∈ 1:length(s1)
    for i2 ∈ 1:length(s2)
      s1[i1] == s2[i2] && return false
    end
  end
  return true
end

"""
Determine if P is a trivial permutation. Errors if P is not a valid
permutation.
"""
function is_trivial_permutation(P)
  isperm(P) || error("Input is not a permutation")
  # TODO: use `all(n->P[n]==n,1:length(P))`?
  N = length(P)
  for n = 1:N
    P[n]!=n && return false
  end
  return true
end

function count_unique(labelsT1,labelsT2)
  count = 0
  for l1 ∈ labelsT1
    l1 ∉ labelsT2 && (count += 1)
  end
  return count
end

function count_common(labelsT1,labelsT2)
  count = 0
  for l1 ∈ labelsT1
    l1 ∈ labelsT2 && (count += 1)
  end
  return count
end

function intersect_positions(labelsT1,labelsT2)
  for i1 = 1:length(labelsT1)
    for i2 = 1:length(labelsT2)
      if labelsT1[i1] == labelsT2[i2]
        return i1,i2
      end
    end
  end
  return nothing
end

function is_replacement(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) == 1 &&
         count_common(labelsT1,labelsT2) == 1
end

function is_combiner(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) == 1 &&
         count_common(labelsT1,labelsT2) > 1
end

function is_uncombiner(labelsT1,labelsT2)
  return count_unique(labelsT1,labelsT2) > 1 &&
         count_common(labelsT1,labelsT2) == 1
end

function readcpp(io::IO,::Type{IndexSet};kwargs...)
  format = get(kwargs,:format,"v3")
  is = IndexSet()
  if format=="v3"
    size = read(io,Int)
    function readind(io,n)
      i = readcpp(io,Index;kwargs...)
      stride = read(io,UInt64)
      return i
    end
    is = IndexSet(ntuple(n->readind(io,n),size))
  else
    throw(ArgumentError("read IndexSet: format=$format not supported"))
  end
  return is
end


text:
#GObject not defined
GtkColorButtonLeaf() = GtkColorButtonLeaf(ccall((:gtk_color_button_new, libgtk), Ptr{GObject}, ()))
GtkColorButtonLeaf(color::GdkRGBA) = GtkColorButtonLeaf(ccall((:gtk_color_button_new_with_rgba, libgtk), Ptr{GObject},
    (Ref{GdkRGBA},), color))

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
"""
    Avro.write([filename|io,] x::T; kw...)

Write an object `x` of avro-supported type `T` in the avro format.
If a file name is provided as a `String` as 1st argument, the avro
data will be written out to disk. Similarly, an `IO` argument can be
provided as 1st argument. If no destination 1st argument is provided,
a byte buffer `Vector{UInt8}` will be returned with avro data written
to it. Supported keyword arguments include:
  * `schema`: the type that should be used when encoding the object in
  the avro format; most common is providing a `Union{...}` type to write
  the data out specifically as a "union type" instead of only the type of the object;
  alternatively, a valid `Avro.Schema` can be passed, like the result of
  `Avro.parseschema(src)`
"""
function write end

function write(obj::T;
        schema=schematype(T),
        jsonencoding::Bool=false,
        kw...) where {T}
    E = jsonencoding ? JSON() : Binary()
    pos = 1
    sch = schematype(schema)
    buf = Vector{UInt8}(undef, nbytes(sch, obj))
    pos = writevalue(E, sch, obj, buf, pos, length(buf), kw)
    return buf
end

function write(io::IO, obj::T; kw...) where {T}
    buf = write(obj; kw...)
    return Base.write(io, buf)
end

function write(fname::String, obj::T; kw...) where {T}
    open(fname, "w") do io
        write(io, obj; kw...)
    end
    fname
end

"""
    Avro.read(source, T_or_sch) => T

Read an avro-encoded object of type `T` or avro schema `sch` from `source`,
which can be a byte buffer `AbstractVector{UInt8}`, file name `String`, or `IO`.

The data in `source` must be avro-formatted data, as no schema verification
can be done. Note that "avro object container files" should be processed
using [`Avro.readtable`](@ref) instead, where the data schema is encoded
in the file itself. Also note that the 2nd argument can be a Julia type like
`Vector{String}`, or a valid `Avro.Schema` type object, like is returned from
`Avro.parseschema(src)`.
"""
function read end

function read(buf::AbstractVector{UInt8}, ::Type{T}; schema=schematype(T), jsonencoding::Bool=false, kw...) where {T}
    x, pos = readvalue(jsonencoding ? JSON() : Binary(), schematype(schema), T, buf, 1, length(buf), kw)
    return x
end

function read(buf::AbstractVector{UInt8}, sch::Schema; T=juliatype(sch), jsonencoding::Bool=false, kw...)
    x, pos = readvalue(jsonencoding ? JSON() : Binary(), sch, T, buf, 1, length(buf), kw)
    return x
end

read(fname::String, sch; kw...) = read(Mmap.mmap(fname), sch; kw...)
read(io::IO, sch; kw...) = read(Base.read(io), sch; kw...)
read(io::IOBuffer, sch; kw...) = read(io.data, sch; kw...)

# fallbacks
writevalue(E::Encoding, x::T, buf, pos, len, opts) where {T} =
    writevalue(E, schematype(T), x, buf, pos, len, opts)
readvalue(E::Encoding, ::Type{T}, buf, pos, len, opts) where {T} =
    readvalue(E, schematype(T), T, buf, pos, len, opts)
readvalue(B::Binary, sch::String, ::Type{T}, buf, pos, len, opts) where {T} =
    readvalue(B, PrimitiveType(sch), T, buf, pos, len, opts)
skipvalue(E::Encoding, ::Type{T}, buf, pos, len, opts) where {T} =
    skipvalue(E, schematype(T), T, buf, pos, len, opts)
skipvalue(B::Binary, sch::String, ::Type{T}, buf, pos, len, opts) where {T} =
    skipvalue(B, PrimitiveType(sch), T, buf, pos, len, opts)

nbytes(x::T) where {T} = nbytes(schematype(T), x)

# null
writevalue(::Binary, ::NullType, x, buf, pos, len, opts) = pos
readvalue(::Binary, ::NullType, ::Type{T}, buf, pos, len, opts) where {T} = (StructTypes.construct(T), pos)
skipvalue(::Binary, ::NullType, ::Type{T}, buf, pos, len, opts) where {T} = pos
nbytes(::NullType, T) = 0

# boolean
function writevalue(::Binary, ::BooleanType, x, buf, pos, len, opts)
    @check 1
    @inbounds buf[pos] = UInt8(Bool(x))
    return pos + 1
end

function readvalue(::Binary, ::BooleanType, ::Type{T}, buf, pos, len, opts) where {T}
    @readcheck 1
    @inbounds x = buf[pos]
    return StructTypes.construct(T, Bool(x)), pos + 1
end

skipvalue(::Binary, ::BooleanType, T, buf, pos, len, opts) = pos + 1
nbytes(::BooleanType, T) = 1

# generic fallbacks to convert to supported number type
const NumberType = Union{IntType, LongType}
const FloatTypes = Union{FloatType, DoubleType}

writevalue(B::Binary, S::NumberType, y::T, buf, pos, len, opts) where {T} =
    writevalue(B, S, StructTypes.construct(StructTypes.numbertype(T), y), buf, pos, len, opts)

function readvalue(B::Binary, S::Union{NumberType, FloatTypes}, ::Type{T}, buf, pos, len, opts) where {T}
    x, pos = _readvalue(B, S, StructTypes.numbertype(T), buf, pos, len, opts)
    return StructTypes.construct(T, x), pos
end

# read/write integers in vint zigzag format: https://lucene.apache.org/core/3_5_0/fileformats.html#VInt
writevalue(B::Binary, N::NumberType, y::T, buf, pos, len, opts) where {T <: Unsigned} =
    writevalue(B, N, signed(widen(y)), buf, pos, len, opts)

function writevalue(::Binary, ::NumberType, y::T, buf, pos, len, opts) where {T <: Signed}
    x = tozigzag(y)
    while true
        # are any bits higher than least 7 set? if not, we're done
        if (x & ((~T(0)) << 7)) == 0
            @inbounds buf[pos] = x % UInt8
            return pos + 1
        end
        @inbounds buf[pos] = ((x & 0b1111111) | 0b10000000) % UInt8
        pos += 1
        x >>>= 7
    end
end

function nbytes(::NumberType, y::T) where {T <: Integer}
    x = tozigzag(y)
    N = 1
    while true
        # are any bits higher than least 7 set? if not, we're done
        if (x & ((~T(0)) << 7)) == 0
            return N
        end
        x >>>= 7
        N += 1
    end
end

function _readvalue(B::Binary, N::NumberType, ::Type{T}, buf, pos, len, opts) where {T <: Unsigned}
    x, pos = _readvalue(B, N, signed(T), buf, pos, len, opts)
    return Core.bitcast(T, x), pos
end

function _readvalue(::Binary, ::NumberType, ::Type{T}, buf, pos, len, opts) where {T <: Signed}
    x = T(0)
    shift = 0
    len = length(buf)
    while pos <= len
        @inbounds b = buf[pos]
        pos += 1
        x |= T(b & 0x7F) << shift
        shift += 7
        (b & 0x80) == 0 && return fromzigzag(x), pos
    end
    return x, pos
end

function skipvalue(::Binary, ::NumberType, ::Type{T}, buf, pos, len, opts) where {T <: Integer}
    len = length(buf)
    @inbounds while pos <= len && (buf[pos] & 0x80) > 0
        pos += 1
    end
    return pos + 1
end

# float/double

function writevalue(::Binary, ::FloatTypes, x::T, buf, pos, len, opts) where {T <: Base.IEEEFloat}
    N = sizeof(T)
    @check N
    ref = Ref(x)
    GC.@preserve ref unsafe_copyto!(pointer(buf, pos), convert(Ptr{UInt8}, Base.unsafe_convert(Ref{T}, ref)), N)
    return pos + N
end

nbytes(::FloatTypes, x::T) where {T <: Base.IEEEFloat} = sizeof(T)

function _readvalue(::Binary, ::FloatTypes, ::Type{T}, buf, pos, len, opts) where {T <: Base.IEEEFloat}
    @readcheck sizeof(T)
    GC.@preserve buf begin
        ptr::Ptr{T} = pointer(buf, pos)
        x = unsafe_load(ptr)
    end
    return x, pos + sizeof(T)
end

skipvalue(::Binary, ::FloatTypes, ::Type{T}, buf, pos, len, opts) where {T <: Base.IEEEFloat} = pos + sizeof(T)

# bytes/strings
_codeunits(x::AbstractVector{UInt8}) = x
_codeunits(x) = codeunits(Base.string(x))

# generic fallbacks to convert to string type
writevalue(B::Binary, S::StringType, x, buf, pos, len, opts) =
    _writevalue(B, S, Base.string(x), buf, pos, len, opts)
writevalue(B::Binary, S::BytesType, x, buf, pos, len, opts) =
    _writevalue(B, S, x, buf, pos, len, opts)

const BytesOrString = Union{BytesType, StringType}

function _writevalue(B::Binary, ::BytesOrString, x, buf, pos, len, opts)
    N = sizeof(x)
    pos = writevalue(B, long, N, buf, pos, len, opts)
    @check N
    for b in _codeunits(x)
        @inbounds buf[pos] = b
        pos += 1
    end
    return pos
end

function _writevalue(B::Binary, ::BytesType, x::AbstractVector{UInt8}, buf, pos, len, opts)
    N = sizeof(x)
    pos = writevalue(B, long, N, buf, pos, len, opts)
    @check N
    copyto!(buf, pos, x, 1, N)
    return pos + N
end

nbytes(::BytesOrString, x) = nbytes(long, sizeof(_codeunits(x))) + sizeof(_codeunits(x))

function readvalue(B::Binary, ::BytesType, ::Type{T}, buf, pos, len, opts) where {T}
    N, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    @readcheck N
    return StructTypes.construct(T, view(buf, pos:(pos + N - 1))), N + pos
end

function skipvalue(B::Binary, ::BytesOrString, ::Type{T}, buf, pos, len, opts) where {T}
    N, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    return N + pos
end

function readvalue(B::Binary, ::StringType, ::Type{T}, buf, pos, len, opts) where {T}
    N, pos = readvalue(B, long, Int64, buf, pos, len, opts)
    @readcheck N
    return StructTypes.construct(T, pointer(buf, pos), N), N + pos
end


text:
"""
    data(dict::Dict, keys::Vector{Symbol})
Construct a tuple with values filled according to `dict` and keys
according to `keys`.
"""
function data(dict::Dict, keys::Vector{Symbol})

    @assert mapreduce(k -> haskey(dict, k), &, keys)

    r = Expr(:tuple)
    for k in keys
        push!(r.args, Expr(:(=), k, dict[k]))
    end
    return Main.eval(r)
end


text:
export cartesian_product

"""
    cartesian_product(X::LazySet, Y::LazySet; [backend]=nothing, [algorithm]::String="vrep")

Compute the Cartesian product of two sets.

### Input

- `X`         -- set
- `Y`         -- another set
- `backend`   -- (optional, default: `nothing`) the polyhedral computations backend
- `algorithm` -- (optional, default: "hrep") the method used to transform each set
                 `X` and `Y` before taking the Cartesian product; choose between
                 "vrep" (use the vertex representation) and "hrep" (use the constraint representation)

### Output

The `VPolytope` (if "vrep" was used) or `HPolytope` (if "hrep" was used)
obtained by the concrete Cartesian product of `X` and `Y`.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).

If `X` can be converted to a one-dimensional interval and the vertices of `Y`
are available use `algorithm="vrep"`.
"""
function cartesian_product(X::LazySet, Y::LazySet; backend=nothing, algorithm::String="hrep")

    if algorithm == "vrep"
        Yv = VPolytope(vertices_list(Y))
        if dim(X) == 1 # special case
            Xv = convert(Interval, X)  # works if X is convex
            return cartesian_product(Xv, Yv)
        end
        Xv = VPolytope(vertices_list(X))
        Pout = cartesian_product(Xv, Yv, backend=backend)

    elseif algorithm == "hrep"
        Xp = HPolyhedron(constraints_list(X))
        if isempty(Xp.constraints) && isuniversal(X)
            Xp = X isa Universe ? X : Universe(dim(X))
        end
        Yp = HPolyhedron(constraints_list(Y))
        if isempty(Yp.constraints) && isuniversal(Y)
            Yp = Y isa Universe ? Y : Universe(dim(Y))
        end
        Pout = cartesian_product(Xp, Yp, backend=backend)

    else
        throw(ArgumentError("expected algorithm `vrep` or `hrep`, got $algorithm"))
    end
    return Pout
end

"""
    cartesian_product(P1::VPolytope, P2::VPolytope; [backend]=nothing)

Compute the Cartesian product of two polytopes in V-representation.

### Input

- `P1`      -- polytope
- `P2`      -- another polytope
- `backend` -- (optional, default: `nothing`) the polyhedral computations backend

### Output

The `VPolytope` obtained by the concrete Cartesian product of `P1` and `P2`.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function cartesian_product(P1::VPolytope, P2::VPolytope; backend=nothing)
    require(:Polyhedra; fun_name="cartesian_product")

    return _cartesian_product_vrep(P1, P2, backend1=backend, backend2=backend)
end

function _cartesian_product_vrep(P1, P2;
                                 backend1=default_polyhedra_backend(P1),
                                 backend2=default_polyhedra_backend(P2))
    if isnothing(backend1)
        backend1 = default_polyhedra_backend(P1)
    end
    if isnothing(backend2)
        backend2 = default_polyhedra_backend(P2)
    end

    P1′ = polyhedron(P1; backend=backend1)
    P2′ = polyhedron(P2; backend=backend2)
    Pout = Polyhedra.vcartesianproduct(P1′, P2′)
    return VPolytope(Pout)
end

"""
    cartesian_product(P1::HPoly, P2::HPoly; [backend]=nothing)

Compute the Cartesian product of two polyhedra in H-representaion.

### Input

- `P1`         -- polyhedron
- `P2`         -- another polyhedron
- `backend`    -- (optional, default: `nothing`) the polyhedral computations backend

### Output

The polyhedron obtained by the concrete cartesian product of `P1` and `P2`.

### Notes

For further information on the supported backends see
[Polyhedra's documentation](https://juliapolyhedra.github.io/).
"""
function cartesian_product(P1::HPoly, P2::HPoly; backend=nothing)
    require(:Polyhedra; fun_name="`cartesian_product")

    return _cartesian_product_hrep(P1, P2, backend1=backend, backend2=backend)
end

function cartesian_product(U1::Universe, U2::Universe)
    return Universe(dim(U1) + dim(U2))
end

function _cartesian_product_hrep(P1::PT1, P2::PT2; backend1, backend2) where {PT1, PT2}
    if isnothing(backend1)
        backend1 = default_polyhedra_backend(P1)
    end
    if isnothing(backend2)
        backend2 = default_polyhedra_backend(P2)
    end

    P1′ = polyhedron(P1; backend=backend1)
    P2′ = polyhedron(P2; backend=backend2)
    Pout = Polyhedra.hcartesianproduct(P1′, P2′)

    PT = isboundedtype(PT1) && isboundedtype(PT2) ? HPolytope : HPolyhedron
    return convert(PT, Pout)
end

# if the first set is an interval => the result is always a polytope
function cartesian_product(I::Interval, P::Union{VPolygon, VPolytope})
    a = min(I)
    b = max(I)
    vlist = vertices_list(P)
    m = length(vlist)

    vout = Vector{eltype(vlist)}(undef, 2m)
    @inbounds for i in 1:m
        vout[i] = vcat(a, vlist[i])
        vout[i + m] = vcat(b, vlist[i])
    end
    return VPolytope(vout)
end

function cartesian_product(S1::AbstractSingleton, S2::AbstractSingleton)
    v1 = element(S1)
    v2 = element(S2)
    return Singleton(vcat(v1, v2))
end

function cartesian_product(Z1::AbstractZonotope, Z2::AbstractZonotope)
    N = promote_type(eltype(Z1), eltype(Z2))
    z1 = zeros(N, dim(Z1), ngens(Z2))
    z2 = zeros(N, dim(Z2), ngens(Z1))
    g1 = genmat(Z1)
    g2 = genmat(Z2)
    c = vcat(center(Z1), center(Z2))
    G = [g1 z1; z2 g2]
    return Zonotope(c, G)
end

function cartesian_product(H1::AbstractHyperrectangle, H2::AbstractHyperrectangle)
    c = vcat(center(H1), center(H2))
    r = vcat(radius_hyperrectangle(H1), radius_hyperrectangle(H2))
    return Hyperrectangle(c, r)
end

function cartesian_product(P::AbstractPolyhedron, U::Universe)
    clist = [cartesian_product(H, U) for H in constraints_list(P)]
    return HPolyhedron(clist)
end

function cartesian_product(U::Universe, P::AbstractPolyhedron)
    clist = [cartesian_product(U, H) for H in constraints_list(P)]
    return HPolyhedron(clist)
end

function cartesian_product(H::HalfSpace, U::Universe)
    return HalfSpace(append_zeros(H.a, dim(U)), H.b)
end

function cartesian_product(U::Universe, H::HalfSpace)
    return HalfSpace(prepend_zeros(H.a, dim(U)), H.b)
end


text:
module TopicTiling
include("Utls.jl")
using Statistics
using .Utls

"""
    TopicTiling.SegmentObject(window_size, smooth_window_size, lda_model, dictionary)
TopicTiling is an extension of TextTiling that uses the topic IDs of words in a sentence to calculate the similarity between blocks.
# Arguments
- `window_size`: Sliding window size.
- `smooth_window_size`: Window size for smoothing depth scores.
- `lda_model`: Trained LDA topic model.
- `dictionary`: A dictionary showing word-id mappings.
"""
mutable struct SegmentObject
    window_size::Int
    smooth_window_size::Int
    lda_model::Any
    dictionary::Any
end

function preprocessing(seg::SegmentObject, tokenized_document)
    n = length(tokenized_document)
    seg.window_size = maximum([minimum([seg.window_size, n / 3]), 1])
    return [
        Utls.count_elements(
            aggregate_topic_id(
                seg.lda_model.get_document_topics(
                    seg.dictionary.doc2bow(tokenized_document[i]),
                    per_word_topics = true,
                )[2],
            ),
        ) for i = 1:n
    ]
end

function aggregate_topic_id(ids)
    sentence_ids = []
    for i in ids
        append!(sentence_ids, string.(i[2]))
    end
    return sentence_ids
end

function calculate_gap_score(seg::SegmentObject, preprocessed_document)
    n = length(preprocessed_document)
    gap_score = [0.0 for _ = 1:n]

    for i = 1:n
        sz = minimum([minimum([i, n - i]), seg.window_size])
        left_side, right_side = Utls.SentenceElements(Dict{String,Int}()),
        Utls.SentenceElements(Dict{String,Int}())

        for j = Int(i - sz + 1):Int(i)
            Utls.merge_elements(left_side, preprocessed_document[j])
        end

        for j = Int(i + 1):Int(i + sz)
            Utls.merge_elements(right_side, preprocessed_document[j])
        end

        gap_score[i] =
            Utls.calculate_cosin_similarity(left_side.elements_dct, right_side.elements_dct)
    end
    return gap_score
end

function calculate_depth_score(seg::SegmentObject, gap_score)
    n = length(gap_score)
    depth_score = [0.0 for _ = 1:n]

    for i = 1:n
        if i <= seg.window_size || i + seg.window_size > n
            continue
        end

        ptr = i - 1
        while ptr >= 1 && gap_score[ptr] >= gap_score[ptr+1]
            ptr -= 1
        end
        lval = gap_score[ptr+1]

        ptr = i + 1
        while ptr < n && gap_score[ptr] >= gap_score[ptr-1]
            ptr += 1
        end
        rval = gap_score[ptr-1]

        depth_score[i] = 1 / 2 * (lval - gap_score[i] + rval - gap_score[i])
    end
    return depth_score
end

function smoothing(seg::SegmentObject, depth_score)
    n = length(depth_score)
    smooth_depth_score = [0.0 for _ = 1:n]

    for i = 1:n
        if i - seg.smooth_window_size < 1 || i + seg.smooth_window_size >= n
            smooth_depth_score[i] = depth_score[i]
        else
            smooth_depth_score[i] =
                mean(depth_score[(i-seg.smooth_window_size):(i+seg.smooth_window_size)])
        end
    end
    return smooth_depth_score
end

function determine_boundaries(seg::SegmentObject, smooth_depth_score, num_topics)
    n = length(smooth_depth_score)
    boundaries = ["0" for _ = 1:n]
    cutoff_threshold = mean(smooth_depth_score) - std(smooth_depth_score) / 2.0

    depth_tuples = [(x, y) for (x, y) in zip(smooth_depth_score, [i for i = 1:n])]
    depth_tuples = reverse(sort(depth_tuples))

    if num_topics == Nothing
        for x in depth_tuples
            if x[1] > cutoff_threshold
                boundaries[x[2]] = "1"
            end
        end
    else
        for (i, x) in enumerate(depth_tuples)
            if x[1] > cutoff_threshold && i <= num_topics - 1
                boundaries[x[2]] = "1"
            end
        end
    end
    return join(boundaries[1:end-1])
end

"""
    TopicTiling.segment(seg, document, [num_topics]) -> String
Performs the splitting of the document entered in the `document` argument.
# Arguments
- `seg`: Segment object.
- `document`: The document to be text segmented.
- `num_topics`: num_topics is the number of topics in the document. If this value is specified, segment boundaries are determined by the number of num_topics, starting with the highest depth score.

# Examples
```jldoctest
using TextSegmentation

# LDA Topic Model
pygensim = pyimport("gensim")
# train_document
# Data to be used when training the LDA topic model.
# Data from the same domain as the text to be segmented is preferred.
tokenized_train_document = [Utls.tokenize(i) for i in train_document]
dictionary = pygensim.corpora.Dictionary(tokenized_train_document)
corpus = [dictionary.doc2bow(text) for text in tokenized_train_document]
lda_model = pygensim.models.ldamodel.LdaModel(
    corpus = corpus,
    id2word = dictionary,
    minimum_probability = 0.0001,
    num_topics = 3,
    random_state=1234,
)

# TopicTiling
window_size = 2
smooth_window_size = 1
num_topics = 3
to = TopicTiling.SegmentObject(window_size, smooth_window_size, lda_model, dictionary)
result = TopicTiling.segment(to, document, num_topics)
println(result)
000110000000
```
"""
function segment(seg, test_document, num_topics=Nothing)
    tokenized_test_document = [Utls.tokenize(i) for i in test_document]
    preprocessed_document = preprocessing(seg, tokenized_test_document)
    gap_score = calculate_gap_score(seg, preprocessed_document)
    depth_score = calculate_depth_score(seg, gap_score)
    smooth_depth_score = smoothing(seg, depth_score)
    return determine_boundaries(seg, smooth_depth_score, num_topics)
end
end


text:
import AutoGrad: ungetindex
using Knet.KnetArrays: DevArray, KnetArray
using AutoGrad: Value

function ungetindex(x::DevArray{T},dxi,i) where T
    if isbitstype(T)
        if dxi isa Value
            forw(addto!, zeroslike(x), forw(ungetindex, x, dxi, i))
        elseif recording()
            addtoindex!(zero(x), dxi, i...)
        else
            Sparse(x,Any[dxi],Any[i])
        end
    else
        # Using addtoindex! instead of setindex! to handle repeated indices
        addtoindex!(Array{Union{T,Nothing}}(nothing, size(x)), dxi, i...)
    end
end